src/float/dfloat/conv/cl_RA_to_double.cc

   1 // double_approx().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/rational.h"
   8
   9
  10 // Implementation.
  11
  12 #include "float/dfloat/cl_DF.h"
  13 #include "rational/cl_RA.h"
  14 #include "cln/integer.h"
  15 #include "integer/cl_I.h"
  16 #include "float/cl_F.h"
  17
  18 namespace cln {
  19
  20 double double_approx (const cl_RA& x)
  21 {
  22 // Method: same as cl_RA_to_DF().
  23       if (integerp(x)) {
  24         DeclareType(cl_I,x);
  25         return double_approx(x);
  26       }
  27  {    // x Ratio
  28       DeclareType(cl_RT,x);
  29       union { dfloat eksplicit; double machine_double; } u;
  30       var cl_I a = numerator(x); // +/- a
  31       var const cl_I& b = denominator(x); // b
  32       var cl_signean sign = -(cl_signean)minusp(a); // Vorzeichen
  33       if (!(sign==0)) { a = -a; } // Betrag nehmen, liefert a
  34       var sintC lendiff = (sintC)integer_length(a) // (integer-length a)
  35                           - (sintC)integer_length(b); // (integer-length b)
  36       if (lendiff > DF_exp_high-DF_exp_mid) // Exponent >= n-m > Obergrenze ?
  37         {
  38           #if (cl_word_size==64)
  39           u.eksplicit =
  40               ((sint64)sign & bit(63))
  41             | ((uint64)(bit(DF_exp_len)-1) << DF_mant_len); // Infinity
  42           #else
  43           u.eksplicit.semhi =
  44               ((sint32)sign & bit(31))
  45             | ((uint32)(bit(DF_exp_len)-1) << (DF_mant_len-32)); // Infinity
  46           u.eksplicit.mlo = 0;
  47           #endif
  48           return u.machine_double;
  49         }
  50       if (lendiff < DF_exp_low-DF_exp_mid-2) // Exponent <= n-m+2 < Untergrenze ?
  51         {
  52           #if (cl_word_size==64)
  53           u.eksplicit = ((sint64)sign & bit(63)); // 0.0
  54           #else
  55           u.eksplicit.semhi = ((sint32)sign & bit(31)); // 0.0
  56           u.eksplicit.mlo = 0;
  57           #endif
  58           return u.machine_double;
  59         }
  60       var cl_I zaehler;
  61       var cl_I nenner;
  62       if (lendiff >= DF_mant_len+2)
  63         // n-m-54>=0
  64         { nenner = ash(b,lendiff - (DF_mant_len+2)); // (ash b n-m-54)
  65           zaehler = a; // a
  66         }
  67         else
  68         { zaehler = ash(a,(DF_mant_len+2) - lendiff); // (ash a -n+m+54)
  69           nenner = b; // b
  70         }
  71       // Division zaehler/nenner durchführen:
  72       var cl_I_div_t q_r = cl_divide(zaehler,nenner);
  73       var cl_I& q = q_r.quotient;
  74       var cl_I& r = q_r.remainder;
  75       // 2^53 <= q < 2^55, also ist q Bignum mit ceiling(55/intDsize) Digits.
  76       var const uintD* ptr = BN_MSDptr(q);
  77       #if (cl_word_size==64)
  78       var uint64 mant = get_max64_Dptr(55,ptr);
  79       if (mant >= bit(DF_mant_len+2))
  80         // 2^54 <= q < 2^55, schiebe um 2 Bits nach rechts
  81         { var uint64 rounding_bits = mant & (bit(2)-1);
  82           lendiff = lendiff+1; // Exponent := n-m+1
  83           mant = mant >> 2;
  84           if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
  85                || ( (rounding_bits == bit(1)) // 10
  86                     && (eq(r,0)) // und genau halbzahlig (r=0)
  87                     && ((mant & bit(0)) ==0) // -> round-to-even
  88              )    )
  89             // abrunden
  90             goto ab;
  91             else
  92             // aufrunden
  93             goto auf;
  94         }
  95         else
  96         { var uint64 rounding_bit = mant & bit(0);
  97           mant = mant >> 1;
  98           if ( (rounding_bit == 0) // 0 wird abgerundet
  99                || ( (eq(r,0)) // genau halbzahlig (r=0)
 100                     && ((mant & bit(0)) ==0) // -> round-to-even
 101              )    )
 102             // abrunden
 103             goto ab;
 104             else
 105             // aufrunden
 106             goto auf;
 107         }
 108       auf:
 109       mant += 1;
 110       if (mant >= bit(DF_mant_len+1)) // rounding overflow?
 111         { mant = mant>>1; lendiff = lendiff+1; }
 112       ab:
 113       // Fertig.
 114       if (lendiff < (sintL)(DF_exp_low-DF_exp_mid))
 115         { u.eksplicit = ((sint64)sign & bit(63)); } // 0.0
 116       else if (lendiff > (sintL)(DF_exp_high-DF_exp_mid))
 117         { u.eksplicit =
 118               ((sint64)sign & bit(63))
 119             | ((uint64)(bit(DF_exp_len)-1) << DF_mant_len); // Infinity
 120         }
 121       else
 122         { u.eksplicit =
 123               ((sint64)sign & bit(63))                      /* Vorzeichen */
 124             | ((uint64)(lendiff+DF_exp_mid) << DF_mant_len) /* Exponent   */
 125             | ((uint64)mant & (bit(DF_mant_len)-1));        /* Mantisse   */
 126         }
 127       return u.machine_double;
 128       #else
 129       var uint32 manthi = get_max32_Dptr(23,ptr);
 130       var uint32 mantlo = get_32_Dptr(ptr mspop ceiling(23,intDsize));
 131       if (manthi >= bit(DF_mant_len-32+2))
 132         // 2^54 <= q < 2^55, schiebe um 2 Bits nach rechts
 133         { var uintL rounding_bits = mantlo & (bit(2)-1);
 134           lendiff = lendiff+1; // Exponent := n-m+1
 135           mantlo = (mantlo >> 2) | (manthi << 30); manthi = manthi >> 2;
 136           if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
 137                || ( (rounding_bits == bit(1)) // 10
 138                     && (eq(r,0)) // und genau halbzahlig (r=0)
 139                     && ((mantlo & bit(0)) ==0) // -> round-to-even
 140              )    )
 141             // abrunden
 142             goto ab;
 143             else
 144             // aufrunden
 145             goto auf;
 146         }
 147         else
 148         { var uintL rounding_bit = mantlo & bit(0);
 149           mantlo = (mantlo >> 1) | (manthi << 31); manthi = manthi >> 1;
 150           if ( (rounding_bit == 0) // 0 wird abgerundet
 151                || ( (eq(r,0)) // genau halbzahlig (r=0)
 152                     && ((mantlo & bit(0)) ==0) // -> round-to-even
 153              )    )
 154             // abrunden
 155             goto ab;
 156             else
 157             // aufrunden
 158             goto auf;
 159         }
 160       auf:
 161       mantlo += 1;
 162       if (mantlo==0)
 163         { manthi += 1;
 164           if (manthi >= bit(DF_mant_len-32+1)) // rounding overflow?
 165             { manthi = manthi>>1; lendiff = lendiff+1; }
 166         }
 167       ab:
 168       // Fertig.
 169       if (lendiff < (sintL)(DF_exp_low-DF_exp_mid))
 170         { u.eksplicit.semhi = ((sint32)sign & bit(31)); // 0.0
 171           u.eksplicit.mlo = 0;
 172         }
 173       else if (lendiff > (sintL)(DF_exp_high-DF_exp_mid))
 174         { u.eksplicit.semhi =
 175               ((sint32)sign & bit(31))
 176             | ((uint32)(bit(DF_exp_len)-1) << (DF_mant_len-32)); // Infinity
 177           u.eksplicit.mlo = 0;
 178         }
 179       else
 180         { u.eksplicit.semhi =
 181               ((sint32)sign & bit(31))                           /* Vorzeichen */
 182             | ((uint32)(lendiff+DF_exp_mid) << (DF_mant_len-32)) /* Exponent   */
 183             | ((uint32)manthi & (bit(DF_mant_len-32)-1));        /* Mantisse   */
 184           u.eksplicit.mlo = mantlo;
 185         }
 186       return u.machine_double;
 187       #endif
 188 }}
 189
 190 }  // namespace cln