src/float/dfloat/elem/cl_DF_from_I.cc

   1 // cl_I_to_DF().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_DF.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_integer.h"
  13 #include "cl_I.h"
  14 #include "cl_DS.h"
  15 #include "cl_F.h"
  16
  17 const cl_DF cl_I_to_DF (const cl_I& x)
  18 {
  19 // Methode:
  20 // x=0 -> Ergebnis 0.0
  21 // Merke Vorzeichen von x.
  22 // x:=(abs x)
  23 // Exponent:=(integer-length x)
  24 //   Greife die 54 höchstwertigen Bits heraus (angeführt von einer 1).
  25 //   Runde das letzte Bit weg:
  26 //     Bit 0 = 0 -> abrunden,
  27 //     Bit 0 = 1 und Rest =0 -> round-to-even,
  28 //     Bit 0 = 1 und Rest >0 -> aufrunden.
  29 //   Dabei um ein Bit nach rechts schieben.
  30 //   Bei Aufrundung auf 2^53 (rounding overflow) Mantisse um 1 Bit nach rechts
  31 //     schieben und Exponent incrementieren.
  32       if (eq(x,0)) { return cl_DF_0; }
  33       var cl_signean sign = -(cl_signean)minusp(x); // Vorzeichen
  34       var cl_I abs_x = (sign==0 ? x : -x);
  35       var uintL exp = integer_length(abs_x); // (integer-length x)
  36       // NDS zu |x|>0 bilden:
  37       var const uintD* MSDptr;
  38       var uintC len;
  39       I_to_NDS_nocopy(abs_x, MSDptr=,len=,,cl_false,);
  40       // MSDptr/len/LSDptr ist die NDS zu x, len>0.
  41       // Führende Digits holen: Brauche DF_mant_len+1 Bits, dazu intDsize
  42       // Bits (die NDS kann mit bis zu intDsize Nullbits anfangen).
  43       // Dann werden diese Bits um (exp mod intDsize) nach rechts geschoben.
  44       var uintD msd = msprefnext(MSDptr); // erstes Digit
  45       #if (cl_word_size==64)
  46       var uint64 msdd = 0; // weitere min(len-1,64/intDsize) Digits
  47       #define NEXT_DIGIT(i)  \
  48         { if (--len == 0) goto ok;                                   \
  49           msdd |= (uint64)msprefnext(MSDptr) << (64-(i+1)*intDsize); \
  50         }
  51       DOCONSTTIMES(64/intDsize,NEXT_DIGIT);
  52       #undef NEXT_DIGIT
  53       #else
  54       var uint32 msdd = 0; // weitere min(len-1,32/intDsize) Digits
  55       var uint32 msddf = 0; // weitere maximal 32/intDsize Digits
  56       #define NEXT_DIGIT(i)  \
  57         { if (--len == 0) goto ok;                                   \
  58           msdd |= (uint32)msprefnext(MSDptr) << (32-(i+1)*intDsize); \
  59         }
  60       DOCONSTTIMES(32/intDsize,NEXT_DIGIT);
  61       #undef NEXT_DIGIT
  62       #define NEXT_DIGIT(i)  \
  63         { if (--len == 0) goto ok;                                    \
  64           msddf |= (uint32)msprefnext(MSDptr) << (32-(i+1)*intDsize); \
  65         }
  66       DOCONSTTIMES(32/intDsize,NEXT_DIGIT);
  67       #undef NEXT_DIGIT
  68       #endif
  69       --len; ok:
  70       #if (cl_word_size==64)
  71       // Die NDS besteht aus msd, msdd und len weiteren Digits.
  72       // Das höchste in 2^64*msd+msdd gesetzte Bit ist Bit Nummer
  73       // 63 + (exp mod intDsize).
  74       var uintL shiftcount = exp % intDsize;
  75       var uint64 mant = // führende 64 Bits
  76         (shiftcount==0
  77           ? msdd
  78           : (((uint64)msd << (64-shiftcount)) | (msdd >> shiftcount))
  79         );
  80       // Das höchste in mant gesetzte Bit ist Bit Nummer 63.
  81       if ( ((mant & bit(62-DF_mant_len)) ==0) // Bit 10 =0 -> abrunden
  82            || ( ((mant & (bit(62-DF_mant_len)-1)) ==0) // Bit 10 =1 und Bits 9..0 =0
  83                 && ((msdd & (bit(shiftcount)-1)) ==0) // und weitere Bits aus msddf =0
  84                 && (!test_loop_msp(MSDptr,len)) // und alle weiteren Digits =0
  85                 // round-to-even, je nach Bit 11 :
  86                 && ((mant & bit(63-DF_mant_len)) ==0)
  87          )    )
  88         // abrunden
  89         { mant = mant >> (63-DF_mant_len); }
  90         else
  91         // aufrunden
  92         { mant = mant >> (63-DF_mant_len);
  93           mant += 1;
  94           if (mant >= bit(DF_mant_len+1)) // rounding overflow?
  95             { mant = mant>>1; exp = exp+1; }
  96         }
  97       return encode_DF(sign,(sintL)exp,mant);
  98       #else
  99       // Die NDS besteht aus msd, msdd, msddf und len weiteren Digits.
 100       // Das höchste in 2^64*msd+2^32*msdd+msddf gesetzte Bit ist Bit Nummer
 101       // 63 + (exp mod intDsize).
 102       var uintL shiftcount = exp % intDsize;
 103       var uint32 manthi; // führende 32 Bits
 104       var uint32 mantlo; // nächste 32 Bits
 105       if (shiftcount==0)
 106         { manthi = msdd; mantlo = msddf; }
 107         else
 108         { manthi = ((uint32)msd << (32-shiftcount)) | (msdd >> shiftcount);
 109           mantlo = (msdd << (32-shiftcount)) | (msddf >> shiftcount);
 110         }
 111       // Das höchste in mant gesetzte Bit ist Bit Nummer 63.
 112       if ( ((mantlo & bit(62-DF_mant_len)) ==0) // Bit 10 =0 -> abrunden
 113            || ( ((mantlo & (bit(62-DF_mant_len)-1)) ==0) // Bit 10 =1 und Bits 9..0 =0
 114                 && ((msddf & (bit(shiftcount)-1)) ==0) // und weitere Bits aus msddf =0
 115                 && (!test_loop_msp(MSDptr,len)) // und alle weiteren Digits =0
 116                 // round-to-even, je nach Bit 11 :
 117                 && ((mantlo & bit(63-DF_mant_len)) ==0)
 118          )    )
 119         // abrunden
 120         { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
 121           manthi = manthi >> (63-DF_mant_len);
 122         }
 123         else
 124         // aufrunden
 125         { mantlo = (mantlo >> (63-DF_mant_len)) | (manthi << (DF_mant_len-32+1));
 126           manthi = manthi >> (63-DF_mant_len);
 127           mantlo += 1;
 128           if (mantlo==0)
 129             { manthi += 1;
 130               if (manthi >= bit(DF_mant_len-32+1)) // rounding overflow?
 131                 { manthi = manthi>>1; exp = exp+1; }
 132         }   }
 133       return encode_DF(sign,(sintL)exp,manthi,mantlo);
 134       #endif
 135 }