src/float/dfloat/conv/cl_I_to_double.cc

   1 // double_approx().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/integer.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_DF.h"
  13 #include "cl_I.h"
  14 #include "cl_DS.h"
  15 #include "cl_F.h"
  16
  17 namespace cln {
  18
  19 double double_approx (const cl_I& x)
  20 {
  21 // Method: same as cl_I_to_DF().
  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 uintL 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=,,cl_false,);
  30       // MSDptr/len/LSDptr ist die NDS zu x, len>0.
  31       // Führende Digits holen: Brauche DF_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 (cl_word_size==64)
  36       var uint64 msdd = 0; // weitere min(len-1,64/intDsize) Digits
  37       #define NEXT_DIGIT(i)  \
  38         { if (--len == 0) goto ok;                                   \
  39           msdd |= (uint64)msprefnext(MSDptr) << (64-(i+1)*intDsize); \
  40         }
  41       DOCONSTTIMES(64/intDsize,NEXT_DIGIT);
  42       #undef NEXT_DIGIT
  43       #else
  44       var uint32 msdd = 0; // weitere min(len-1,32/intDsize) Digits
  45       var uint32 msddf = 0; // weitere maximal 32/intDsize Digits
  46       #define NEXT_DIGIT(i)  \
  47         { if (--len == 0) goto ok;                                   \
  48           msdd |= (uint32)msprefnext(MSDptr) << (32-(i+1)*intDsize); \
  49         }
  50       DOCONSTTIMES(32/intDsize,NEXT_DIGIT);
  51       #undef NEXT_DIGIT
  52       #define NEXT_DIGIT(i)  \
  53         { if (--len == 0) goto ok;                                    \
  54           msddf |= (uint32)msprefnext(MSDptr) << (32-(i+1)*intDsize); \
  55         }
  56       DOCONSTTIMES(32/intDsize,NEXT_DIGIT);
  57       #undef NEXT_DIGIT
  58       #endif
  59       --len; ok:
  60       #if (cl_word_size==64)
  61       // Die NDS besteht aus msd, msdd und len weiteren Digits.
  62       // Das höchste in 2^64*msd+msdd gesetzte Bit ist Bit Nummer
  63       // 63 + (exp mod intDsize).
  64       var uintL shiftcount = exp % intDsize;
  65       var uint64 mant = // führende 64 Bits
  66         (shiftcount==0
  67           ? msdd
  68           : (((uint64)msd << (64-shiftcount)) | (msdd >> shiftcount))
  69         );
  70       // Das höchste in mant gesetzte Bit ist Bit Nummer 63.
  71       if ( ((mant & bit(62-DF_mant_len)) ==0) // Bit 10 =0 -> abrunden
  72            || ( ((mant & (bit(62-DF_mant_len)-1)) ==0) // Bit 10 =1 und Bits 9..0 =0
  73                 && ((msdd & (bit(shiftcount)-1)) ==0) // und weitere Bits aus msddf =0
  74                 && (!test_loop_msp(MSDptr,len)) // und alle weiteren Digits =0
  75                 // round-to-even, je nach Bit 11 :
  76                 && ((mant & bit(63-DF_mant_len)) ==0)
  77          )    )
  78         // abrunden
  79         { mant = mant >> (63-DF_mant_len); }
  80         else
  81         // aufrunden
  82         { mant = mant >> (63-DF_mant_len);
  83           mant += 1;
  84           if (mant >= bit(DF_mant_len+1)) // rounding overflow?
  85             { mant = mant>>1; exp = exp+1; }
  86         }
  87       union { dfloat eksplicit; double machine_double; } u;
  88       if ((sintL)exp > (sintL)(DF_exp_high-DF_exp_mid))
  89         { u.eksplicit =
  90               ((sint64)sign & bit(63))
  91             | ((uint64)(bit(DF_exp_len)-1) << DF_mant_len); // Infinity
  92         }
  93       else
  94         { u.eksplicit =
  95               ((sint64)sign & bit(63))                         /* Vorzeichen */
  96             | ((uint64)((sintL)exp+DF_exp_mid) << DF_mant_len) /* Exponent   */
  97             | ((uint64)mant & (bit(DF_mant_len)-1));           /* Mantisse   */
  98         }
  99       return u.machine_double;
 100       #else
 101       // Die NDS besteht aus msd, msdd, msddf und len weiteren Digits.
 102       // Das höchste in 2^64*msd+2^32*msdd+msddf gesetzte Bit ist Bit Nummer
 103       // 63 + (exp mod intDsize).
 104       var uintL shiftcount = exp % intDsize;
 105       var uint32 manthi; // führende 32 Bits
 106       var uint32 mantlo; // nächste 32 Bits
 107       if (shiftcount==0)
 108         { manthi = msdd; mantlo = msddf; }
 109         else
 110         { manthi = ((uint32)msd << (32-shiftcount)) | (msdd >> shiftcount);
 111           mantlo = (msdd << (32-shiftcount)) | (msddf >> shiftcount);
 112         }
 113       // Das höchste in mant gesetzte Bit ist Bit Nummer 63.
 114       if ( ((mantlo & bit(62-DF_mant_len)) ==0) // Bit 10 =0 -> abrunden
 115            || ( ((mantlo & (bit(62-DF_mant_len)-1)) ==0) // Bit 10 =1 und Bits 9..0 =0
 116                 && ((msddf & (bit(shiftcount)-1)) ==0) // und weitere Bits aus msddf =0
 117                 && (!test_loop_msp(MSDptr,len)) // und alle weiteren Digits =0
 118                 // round-to-even, je nach Bit 11 :
 119                 && ((mantlo & bit(63-DF_mant_len)) ==0)
 120          )    )
 121         // abrunden
 122         { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
 123           manthi = manthi >> (63-DF_mant_len);
 124         }
 125         else
 126         // aufrunden
 127         { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
 128           manthi = manthi >> (63-DF_mant_len);
 129           mantlo += 1;
 130           if (mantlo==0)
 131             { manthi += 1;
 132               if (manthi >= bit(DF_mant_len-32+1)) // rounding overflow?
 133                 { manthi = manthi>>1; exp = exp+1; }
 134         }   }
 135       union { dfloat eksplicit; double machine_double; } u;
 136       if ((sintL)exp > (sintL)(DF_exp_high-DF_exp_mid))
 137         { u.eksplicit.semhi =
 138               ((sint32)sign & bit(31))
 139             | ((uint32)(bit(DF_exp_len)-1) << (DF_mant_len-32)); // Infinity
 140           u.eksplicit.mlo = 0;
 141         }
 142       else
 143         { u.eksplicit.semhi =
 144               ((sint32)sign & bit(31))                              /* Vorzeichen */
 145             | ((uint32)((sintL)exp+DF_exp_mid) << (DF_mant_len-32)) /* Exponent   */
 146             | ((uint32)manthi & (bit(DF_mant_len-32)-1));           /* Mantisse   */
 147           u.eksplicit.mlo = mantlo;
 148         }
 149       return u.machine_double;
 150       #endif
 151 }
 152
 153 }  // namespace cln