src/float/ffloat/elem/cl_FF_from_RA.cc

   1 // cl_RA_to_FF().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_FF.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_RA.h"
  13 #include "cln/integer.h"
  14 #include "cl_I.h"
  15 #include "cl_F.h"
  16
  17 namespace cln {
  18
  19 const cl_FF cl_RA_to_FF (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+25)*a/b) :
  29 //   Bei n-m>=25 dividiere a durch (ash b (n-m-25)),
  30 //   bei n-m<25 dividiere (ash a (-n+m+25)) durch b.
  31 //   Der erste Wert ist >=2^24, <2^26.
  32 //   Falls er >=2^25 ist, runde 2 Bits weg,
  33 //   falls er <2^25 ist, runde 1 Bit weg.
  34       if (integerp(x)) {
  35         DeclareType(cl_I,x);
  36         return cl_I_to_FF(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 sintL lendiff = (sintL)integer_length(a) // (integer-length a)
  45                           - (sintL)integer_length(b); // (integer-length b)
  46       if (lendiff > FF_exp_high-FF_exp_mid) // Exponent >= n-m > Obergrenze ?
  47         { cl_error_floating_point_overflow(); } // -> Overflow
  48       if (lendiff < FF_exp_low-FF_exp_mid-2) // Exponent <= n-m+2 < Untergrenze ?
  49         { if (underflow_allowed())
  50             { cl_error_floating_point_underflow(); } // -> Underflow
  51             else
  52             { return cl_FF_0; }
  53         }
  54       var cl_I zaehler;
  55       var cl_I nenner;
  56       if (lendiff >= FF_mant_len+2)
  57         // n-m-25>=0
  58         { nenner = ash(b,lendiff - (FF_mant_len+2)); // (ash b n-m-25)
  59           zaehler = a; // a
  60         }
  61         else
  62         { zaehler = ash(a,(FF_mant_len+2) - lendiff); // (ash a -n+m+25)
  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       // 2^24 <= q < 2^26, also ist q Fixnum oder Bignum mit bn_minlength Digits.
  70       var uint32 mant = ((FF_mant_len+3 < cl_value_len)
  71                           ? FN_to_UL(q)
  72                           : cl_I_to_UL(q)
  73                         );
  74       if (mant >= bit(FF_mant_len+2))
  75         // 2^25 <= q < 2^26, schiebe um 2 Bits nach rechts
  76         { var uintL rounding_bits = mant & (bit(2)-1);
  77           lendiff = lendiff+1; // Exponent := n-m+1
  78           mant = mant >> 2;
  79           if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
  80                || ( (rounding_bits == bit(1)) // 10
  81                     && (eq(r,0)) // und genau halbzahlig (r=0)
  82                     && ((mant & bit(0)) ==0) // -> round-to-even
  83              )    )
  84             // abrunden
  85             goto ab;
  86             else
  87             // aufrunden
  88             goto auf;
  89         }
  90         else
  91         { var uintL rounding_bit = mant & bit(0);
  92           mant = mant >> 1;
  93           if ( (rounding_bit == 0) // 0 wird abgerundet
  94                || ( (eq(r,0)) // genau halbzahlig (r=0)
  95                     && ((mant & bit(0)) ==0) // -> round-to-even
  96              )    )
  97             // abrunden
  98             goto ab;
  99             else
 100             // aufrunden
 101             goto auf;
 102         }
 103       auf:
 104       mant += 1;
 105       if (mant >= bit(FF_mant_len+1)) // rounding overflow?
 106         { mant = mant>>1; lendiff = lendiff+1; }
 107       ab:
 108       // Fertig.
 109       return encode_FF(sign,lendiff,mant);
 110 }}
 111
 112 }  // namespace cln