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 "cl_integer.h"
  14 #include "cl_I.h"
  15 #include "cl_F.h"
  16
  17 const cl_FF cl_RA_to_FF (const cl_RA& x)
  18 {
  19 // Methode:
  20 // x ganz -> klar.
  21 // x = +/- a/b mit Integers a,b>0:
  22 //   Seien n,m so gewählt, daß
  23 //     2^(n-1) <= a < 2^n, 2^(m-1) <= b < 2^m.
  24 //   Dann ist 2^(n-m-1) < a/b < 2^(n-m+1).
  25 //   Berechne n=(integer-length a) und m=(integer-length b) und
  26 //   floor(2^(-n+m+25)*a/b) :
  27 //   Bei n-m>=25 dividiere a durch (ash b (n-m-25)),
  28 //   bei n-m<25 dividiere (ash a (-n+m+25)) durch b.
  29 //   Der erste Wert ist >=2^24, <2^26.
  30 //   Falls er >=2^25 ist, runde 2 Bits weg,
  31 //   falls er <2^25 ist, runde 1 Bit weg.
  32       if (integerp(x)) {
  33         DeclareType(cl_I,x);
  34         return cl_I_to_FF(x);
  35       }
  36  {    // x Ratio
  37       DeclareType(cl_RT,x);
  38       var cl_I a = numerator(x); // +/- a
  39       var const cl_I& b = denominator(x); // b
  40       var cl_signean sign = -(cl_signean)minusp(a); // Vorzeichen
  41       if (!(sign==0)) { a = -a; } // Betrag nehmen, liefert a
  42       var sintL lendiff = (sintL)integer_length(a) // (integer-length a)
  43                           - (sintL)integer_length(b); // (integer-length b)
  44       if (lendiff > FF_exp_high-FF_exp_mid) // Exponent >= n-m > Obergrenze ?
  45         { cl_error_floating_point_overflow(); } // -> Overflow
  46       if (lendiff < FF_exp_low-FF_exp_mid-2) // Exponent <= n-m+2 < Untergrenze ?
  47         { if (underflow_allowed())
  48             { cl_error_floating_point_underflow(); } // -> Underflow
  49             else
  50             { return cl_FF_0; }
  51         }
  52       var cl_I zaehler;
  53       var cl_I nenner;
  54       if (lendiff >= FF_mant_len+2)
  55         // n-m-25>=0
  56         { nenner = ash(b,lendiff - (FF_mant_len+2)); // (ash b n-m-25)
  57           zaehler = a; // a
  58         }
  59         else
  60         { zaehler = ash(a,(FF_mant_len+2) - lendiff); // (ash a -n+m+25)
  61           nenner = b; // b
  62         }
  63       // Division zaehler/nenner durchführen:
  64       var cl_I_div_t q_r = cl_divide(zaehler,nenner);
  65       var cl_I& q = q_r.quotient;
  66       var cl_I& r = q_r.remainder;
  67       // 2^24 <= q < 2^26, also ist q Fixnum oder Bignum mit bn_minlength Digits.
  68       var uint32 mant = ((FF_mant_len+3 < cl_value_len)
  69                           ? FN_to_UL(q)
  70                           : cl_I_to_UL(q)
  71                         );
  72       if (mant >= bit(FF_mant_len+2))
  73         // 2^25 <= q < 2^26, schiebe um 2 Bits nach rechts
  74         { var uintL rounding_bits = mant & (bit(2)-1);
  75           lendiff = lendiff+1; // Exponent := n-m+1
  76           mant = mant >> 2;
  77           if ( (rounding_bits < bit(1)) // 00,01 werden abgerundet
  78                || ( (rounding_bits == bit(1)) // 10
  79                     && (eq(r,0)) // und genau halbzahlig (r=0)
  80                     && ((mant & bit(0)) ==0) // -> round-to-even
  81              )    )
  82             // abrunden
  83             goto ab;
  84             else
  85             // aufrunden
  86             goto auf;
  87         }
  88         else
  89         { var uintL rounding_bit = mant & bit(0);
  90           mant = mant >> 1;
  91           if ( (rounding_bit == 0) // 0 wird abgerundet
  92                || ( (eq(r,0)) // genau halbzahlig (r=0)
  93                     && ((mant & bit(0)) ==0) // -> round-to-even
  94              )    )
  95             // abrunden
  96             goto ab;
  97             else
  98             // aufrunden
  99             goto auf;
 100         }
 101       auf:
 102       mant += 1;
 103       if (mant >= bit(FF_mant_len+1)) // rounding overflow?
 104         { mant = mant>>1; lendiff = lendiff+1; }
 105       ab:
 106       // Fertig.
 107       return encode_FF(sign,lendiff,mant);
 108 }}