src/float/ffloat/conv/cl_RA_to_float.cc

   1 // float_approx().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/rational.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_FF.h"
  13 #include "cl_RA.h"
  14 #include "cln/integer.h"
  15 #include "cl_I.h"
  16 #include "cl_F.h"
  17
  18 namespace cln {
  19
  20 float float_approx (const cl_RA& x)
  21 {
  22 // Method: same as cl_RA_to_FF().
  23       if (integerp(x)) {
  24         DeclareType(cl_I,x);
  25         return float_approx(x);
  26       }
  27  {    // x Ratio
  28       DeclareType(cl_RT,x);
  29       union { ffloat eksplicit; float machine_float; } 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 > FF_exp_high-FF_exp_mid) // Exponent >= n-m > Obergrenze ?
  37         { u.eksplicit = make_FF_word(sign,bit(FF_exp_len)-1,0); // Infinity
  38           return u.machine_float;
  39         }
  40       if (lendiff < FF_exp_low-FF_exp_mid-2) // Exponent <= n-m+2 < Untergrenze ?
  41         { u.eksplicit = make_FF_word(sign,0,0); // 0.0
  42           return u.machine_float;
  43         }
  44       var cl_I zaehler;
  45       var cl_I nenner;
  46       if (lendiff >= FF_mant_len+2)
  47         // n-m-25>=0
  48         { nenner = ash(b,lendiff - (FF_mant_len+2)); // (ash b n-m-25)
  49           zaehler = a; // a
  50         }
  51         else
  52         { zaehler = ash(a,(FF_mant_len+2) - lendiff); // (ash a -n+m+25)
  53           nenner = b; // b
  54         }
  55       // Division zaehler/nenner durchführen:
  56       var cl_I_div_t q_r = cl_divide(zaehler,nenner);
  57       var cl_I& q = q_r.quotient;
  58       var cl_I& r = q_r.remainder;
  59       // 2^24 <= q < 2^26, also ist q Fixnum oder Bignum mit bn_minlength Digits.
  60       var uint32 mant = ((FF_mant_len+3 < cl_value_len)
  61                           ? FN_to_UV(q)
  62                           : cl_I_to_UL(q)
  63                         );
  64       if (mant >= bit(FF_mant_len+2))
  65         // 2^25 <= q < 2^26, schiebe um 2 Bits nach rechts
  66         { var uintL rounding_bits = mant & (bit(2)-1);
  67           lendiff = lendiff+1; // Exponent := n-m+1
  68           mant = mant >> 2;
  69           if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
  70                || ( (rounding_bits == bit(1)) // 10
  71                     && (eq(r,0)) // und genau halbzahlig (r=0)
  72                     && ((mant & bit(0)) ==0) // -> round-to-even
  73              )    )
  74             // abrunden
  75             goto ab;
  76             else
  77             // aufrunden
  78             goto auf;
  79         }
  80         else
  81         { var uintL rounding_bit = mant & bit(0);
  82           mant = mant >> 1;
  83           if ( (rounding_bit == 0) // 0 wird abgerundet
  84                || ( (eq(r,0)) // 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       auf:
  94       mant += 1;
  95       if (mant >= bit(FF_mant_len+1)) // rounding overflow?
  96         { mant = mant>>1; lendiff = lendiff+1; }
  97       ab:
  98       // Fertig.
  99       if (lendiff < (sintL)(FF_exp_low-FF_exp_mid))
 100         { u.eksplicit = make_FF_word(sign,0,0); }
 101       else if (lendiff > (sintL)(FF_exp_high-FF_exp_mid))
 102         { u.eksplicit = make_FF_word(sign,bit(FF_exp_len)-1,0); } // Infinity
 103       else
 104         { u.eksplicit = make_FF_word(sign,lendiff+FF_exp_mid,mant); }
 105       return u.machine_float;
 106 }}
 107
 108 }  // namespace cln