src/float/dfloat/elem/cl_DF_from_RA.cc

   1 // cl_RA_to_DF().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "float/dfloat/cl_DF.h"
   8
   9
  10 // Implementation.
  11
  12 #include "rational/cl_RA.h"
  13 #include "cln/integer.h"
  14 #include "integer/cl_I.h"
  15 #include "float/cl_F.h"
  16
  17 namespace cln {
  18
  19 const cl_DF cl_RA_to_DF (const cl_RA& x)
  20 {
  21 // Methode:
  22 // x ganz -> klar.
  23 // x = +/- a/b mit Integers a,b>0:
  24 //   Seien n,m so gewählt, daß
  25 //     2^(n-1) <= a < 2^n, 2^(m-1) <= b < 2^m.
  26 //   Dann ist 2^(n-m-1) < a/b < 2^(n-m+1).
  27 //   Berechne n=(integer-length a) und m=(integer-length b) und
  28 //   floor(2^(-n+m+54)*a/b) :
  29 //   Bei n-m>=54 dividiere a durch (ash b (n-m-54)),
  30 //   bei n-m<54 dividiere (ash a (-n+m+54)) durch b.
  31 //   Der erste Wert ist >=2^53, <2^55.
  32 //   Falls er >=2^54 ist, runde 2 Bits weg,
  33 //   falls er <2^54 ist, runde 1 Bit weg.
  34       if (integerp(x)) {
  35         DeclareType(cl_I,x);
  36         return cl_I_to_DF(x);
  37       }
  38  {    // x Ratio
  39       DeclareType(cl_RT,x);
  40       var cl_I a = numerator(x); // +/- a
  41       var const cl_I& b = denominator(x); // b
  42       var cl_signean sign = -(cl_signean)minusp(a); // Vorzeichen
  43       if (!(sign==0)) { a = -a; } // Betrag nehmen, liefert a
  44       var sintC lendiff = (sintC)integer_length(a) // (integer-length a)
  45                           - (sintC)integer_length(b); // (integer-length b)
  46       if (lendiff > DF_exp_high-DF_exp_mid) // Exponent >= n-m > Obergrenze ?
  47         { throw floating_point_overflow_exception(); } // -> Overflow
  48       if (lendiff < DF_exp_low-DF_exp_mid-2) // Exponent <= n-m+2 < Untergrenze ?
  49         { if (underflow_allowed())
  50             { throw floating_point_underflow_exception(); } // -> Underflow
  51             else
  52             { return cl_DF_0; }
  53         }
  54       var cl_I zaehler;
  55       var cl_I nenner;
  56       if (lendiff >= DF_mant_len+2)
  57         // n-m-54>=0
  58         { nenner = ash(b,lendiff - (DF_mant_len+2)); // (ash b n-m-54)
  59           zaehler = a; // a
  60         }
  61         else
  62         { zaehler = ash(a,(DF_mant_len+2) - lendiff); // (ash a -n+m+54)
  63           nenner = b; // b
  64         }
  65       // Division zaehler/nenner durchführen:
  66       var cl_I_div_t q_r = cl_divide(zaehler,nenner);
  67       var cl_I& q = q_r.quotient;
  68       var cl_I& r = q_r.remainder;
  69       #if (cl_word_size==64)
  70       #if (cl_value_len>=55)
  71       // 2^53 <= q < 2^55: q is fixnum.
  72       var uint64 mant = FN_to_UV(q);
  73       #else
  74       // 2^53 <= q < 2^55: q is bignum with one Digit.
  75       var const uintD* ptr = BN_MSDptr(q);
  76       var uint64 mant = get_max64_Dptr(55,ptr);
  77       #endif
  78       if (mant >= bit(DF_mant_len+2))
  79         // 2^54 <= q < 2^55, schiebe um 2 Bits nach rechts
  80         { var uint64 rounding_bits = mant & (bit(2)-1);
  81           lendiff = lendiff+1; // Exponent := n-m+1
  82           mant = mant >> 2;
  83           if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
  84                || ( (rounding_bits == bit(1)) // 10
  85                     && (eq(r,0)) // und genau halbzahlig (r=0)
  86                     && ((mant & bit(0)) ==0) // -> round-to-even
  87              )    )
  88             // abrunden
  89             goto ab;
  90             else
  91             // aufrunden
  92             goto auf;
  93         }
  94         else
  95         { var uint64 rounding_bit = mant & bit(0);
  96           mant = mant >> 1;
  97           if ( (rounding_bit == 0) // 0 wird abgerundet
  98                || ( (eq(r,0)) // genau halbzahlig (r=0)
  99                     && ((mant & bit(0)) ==0) // -> round-to-even
 100              )    )
 101             // abrunden
 102             goto ab;
 103             else
 104             // aufrunden
 105             goto auf;
 106         }
 107       auf:
 108       mant += 1;
 109       if (mant >= bit(DF_mant_len+1)) // rounding overflow?
 110         { mant = mant>>1; lendiff = lendiff+1; }
 111       ab:
 112       // Fertig.
 113       return encode_DF(sign,lendiff,mant);
 114       #else  // (cl_word_size<64)
 115       // 2^53 <= q < 2^55: q is bignum with two Digits.
 116       var const uintD* ptr = BN_MSDptr(q);
 117       var uint32 manthi = get_max32_Dptr(23,ptr);
 118       var uint32 mantlo = get_32_Dptr(ptr mspop ceiling(23,intDsize));
 119       if (manthi >= bit(DF_mant_len-32+2))
 120         // 2^54 <= q < 2^55, schiebe um 2 Bits nach rechts
 121         { var uintL rounding_bits = mantlo & (bit(2)-1);
 122           lendiff = lendiff+1; // Exponent := n-m+1
 123           mantlo = (mantlo >> 2) | (manthi << 30); manthi = manthi >> 2;
 124           if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
 125                || ( (rounding_bits == bit(1)) // 10
 126                     && (eq(r,0)) // und genau halbzahlig (r=0)
 127                     && ((mantlo & bit(0)) ==0) // -> round-to-even
 128              )    )
 129             // abrunden
 130             goto ab;
 131             else
 132             // aufrunden
 133             goto auf;
 134         }
 135         else
 136         { var uintL rounding_bit = mantlo & bit(0);
 137           mantlo = (mantlo >> 1) | (manthi << 31); manthi = manthi >> 1;
 138           if ( (rounding_bit == 0) // 0 wird abgerundet
 139                || ( (eq(r,0)) // genau halbzahlig (r=0)
 140                     && ((mantlo & bit(0)) ==0) // -> round-to-even
 141              )    )
 142             // abrunden
 143             goto ab;
 144             else
 145             // aufrunden
 146             goto auf;
 147         }
 148       auf:
 149       mantlo += 1;
 150       if (mantlo==0)
 151         { manthi += 1;
 152           if (manthi >= bit(DF_mant_len-32+1)) // rounding overflow?
 153             { manthi = manthi>>1; lendiff = lendiff+1; }
 154         }
 155       ab:
 156       // Fertig.
 157       return encode_DF(sign,lendiff,manthi,mantlo);
 158       #endif
 159 }}
 160
 161 }  // namespace cln