src/float/dfloat/conv/cl_RA_to_double.cc

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