src/float/sfloat/elem/cl_SF_from_RA.cc

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