src/float/ffloat/conv/cl_I_to_float.cc

   1 // float_approx().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/integer.h"
   8
   9
  10 // Implementation.
  11
  12 #include "float/ffloat/cl_FF.h"
  13 #include "integer/cl_I.h"
  14 #include "base/digitseq/cl_DS.h"
  15 #include "float/cl_F.h"
  16
  17 namespace cln {
  18
  19 float float_approx (const cl_I& x)
  20 {
  21 // Method: same as cl_I_to_FF().
  22       if (eq(x,0)) { return 0.0; }
  23       var cl_signean sign = -(cl_signean)minusp(x); // Vorzeichen
  24       var cl_I abs_x = (sign==0 ? x : -x);
  25       var uintC exp = integer_length(abs_x); // (integer-length x)
  26       // NDS zu |x|>0 bilden:
  27       var const uintD* MSDptr;
  28       var uintC len;
  29       I_to_NDS_nocopy(abs_x, MSDptr=,len=,,false,);
  30       // MSDptr/len/LSDptr ist die NDS zu x, len>0.
  31       // Führende Digits holen: Brauche FF_mant_len+1 Bits, dazu intDsize
  32       // Bits (die NDS kann mit bis zu intDsize Nullbits anfangen).
  33       // Dann werden diese Bits um (exp mod intDsize) nach rechts geschoben.
  34       var uintD msd = msprefnext(MSDptr); // erstes Digit
  35       #if (intDsize==64)
  36       var uintD msdd = 0; // weiteres Digit
  37       if (--len == 0) goto ok;
  38       msdd = msprefnext(MSDptr);
  39       #else // (intDsize<=32)
  40       var uint32 msdd = 0; // weitere min(len-1,32/intDsize) Digits
  41       #define NEXT_DIGIT(i)  \
  42         { if (--len == 0) goto ok;                                   \
  43           msdd |= (uint32)msprefnext(MSDptr) << (32-(i+1)*intDsize); \
  44         }
  45       DOCONSTTIMES(32/intDsize,NEXT_DIGIT);
  46       #undef NEXT_DIGIT
  47       #endif
  48       --len; ok:
  49       #if (intDsize==64)
  50       // Die NDS besteht aus msd, msdd, und len weiteren Digits.
  51       // Das höchste in 2^intDsize*msd+msdd gesetzte Bit ist Bit Nummer
  52       // intDsize-1 + (exp mod intDsize).
  53       var uintL shiftcount = exp % intDsize;
  54       var uint64 mant = // führende 64 Bits
  55         (shiftcount==0
  56          ? msdd
  57          : ((msd << (64-shiftcount)) | (msdd >> shiftcount))
  58         );
  59       // Das höchste in mant gesetzte Bit ist Bit Nummer 63.
  60       if ( ((mant & bit(62-FF_mant_len)) ==0) // Bit 39 =0 -> abrunden
  61            || ( ((mant & (bit(62-FF_mant_len)-1)) ==0) // Bit 39 =1 und Bits 38..0 =0
  62                 && ((msdd & (bit(shiftcount)-1)) ==0) // und weitere Bits aus msdd =0
  63                 && (!test_loop_msp(MSDptr,len)) // und alle weiteren Digits =0
  64                 // round-to-even, je nach Bit 40 :
  65                 && ((mant & bit(63-FF_mant_len)) ==0)
  66          )    )
  67         // abrunden
  68         { mant = mant >> (63-FF_mant_len); }
  69         else
  70         // aufrunden
  71         { mant = mant >> (63-FF_mant_len);
  72           mant += 1;
  73           if (mant >= bit(FF_mant_len+1)) // rounding overflow?
  74             { mant = mant>>1; exp = exp+1; }
  75         }
  76       #else
  77       // Die NDS besteht aus msd, msdd, und len weiteren Digits.
  78       // Das höchste in 2^32*msd+msdd gesetzte Bit ist Bit Nummer
  79       // 31 + (exp mod intDsize).
  80       var uintL shiftcount = exp % intDsize;
  81       var uint32 mant = // führende 32 Bits
  82         (shiftcount==0
  83          ? msdd
  84          : (((uint32)msd << (32-shiftcount)) | (msdd >> shiftcount))
  85         );
  86       // Das höchste in mant gesetzte Bit ist Bit Nummer 31.
  87       if ( ((mant & bit(30-FF_mant_len)) ==0) // Bit 7 =0 -> abrunden
  88            || ( ((mant & (bit(30-FF_mant_len)-1)) ==0) // Bit 7 =1 und Bits 6..0 =0
  89                 && ((msdd & (bit(shiftcount)-1)) ==0) // und weitere Bits aus msdd =0
  90                 && (!test_loop_msp(MSDptr,len)) // und alle weiteren Digits =0
  91                 // round-to-even, je nach Bit 8 :
  92                 && ((mant & bit(31-FF_mant_len)) ==0)
  93          )    )
  94         // abrunden
  95         { mant = mant >> (31-FF_mant_len); }
  96         else
  97         // aufrunden
  98         { mant = mant >> (31-FF_mant_len);
  99           mant += 1;
 100           if (mant >= bit(FF_mant_len+1)) // rounding overflow?
 101             { mant = mant>>1; exp = exp+1; }
 102         }
 103       #endif
 104       union { ffloat eksplicit; float machine_float; } u;
 105       if ((sintL)exp > (sintL)(FF_exp_high-FF_exp_mid))
 106         { u.eksplicit = make_FF_word(sign,bit(FF_exp_len)-1,0); } // Infinity
 107       else
 108         { u.eksplicit = make_FF_word(sign,(sintL)exp+FF_exp_mid,mant); }
 109       return u.machine_float;
 110 }
 111
 112 }  // namespace cln