src/float/lfloat/elem/cl_LF_from_RA.cc

   1 // cl_RA_to_LF().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_LF.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_LF_impl.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 const cl_LF cl_RA_to_LF (const cl_RA& x, uintC len)
  21 {
  22 // Methode:
  23 // x ganz -> klar.
  24 // x = +/- a/b mit Integers a,b>0:
  25 //   Sei k,m so gewählt, daß
  26 //     2^(k-1) <= a < 2^k, 2^(m-1) <= b < 2^m.
  27 //   Dann ist 2^(k-m-1) < a/b < 2^(k-m+1).
  28 //   Ergebnis-Vorzeichen := Vorzeichen von x.
  29 //   Berechne k=(integer-length a) und m=(integer-length b).
  30 //   Ergebnis-Exponent := k-m.
  31 //   Ergebnis-Mantisse:
  32 //     Berechne floor(2^(-k+m+16n+1)*a/b) :
  33 //       Bei k-m>=16n+1 dividiere a durch (ash b (k-m-16n-1)),
  34 //       bei k-m<16n+1 dividiere (ash a (-k+m+16n+1)) durch b.
  35 //     Der erste Wert ist >=2^16n, <2^(16n+2).
  36 //     Falls er >=2^(16n+1) ist, erhöhe Exponent um 1,
  37 //       runde 2 Bits weg und schiebe dabei um 2 Bits nach rechts;
  38 //     falls er <2^(16n+1) ist,
  39 //       runde 1 Bit weg und schiebe dabei um 1 Bit nach rechts.
  40 // NB: Wenn a und b länger sind als len, ist dieser Algorithmus weniger
  41 //     effizient, als cl_float(a,len)/cl_float(b,len) zu berechnen. Aber
  42 //     es ist wichtig, dass cl_RA_to_LF nicht mehr als 0.5 ulp Fehler hat,
  43 //     deswegen belassen wir es beim ineffizienten aber exakten Algorithmus.
  44 //     Wenn es auf Rundungsfehler nicht ankommt, muss der Aufrufer im Fall
  45 //           ceiling(integer_length(a),intDsize) >= len
  46 //        && ceiling(integer_length(b),intDsize) >= len
  47 //     einen anderen Algorithmus wählen.
  48       if (integerp(x)) {
  49         DeclareType(cl_I,x);
  50         return cl_I_to_LF(x,len);
  51       }
  52  {    // x Ratio
  53       DeclareType(cl_RT,x);
  54       var cl_I a = numerator(x); // +/- a
  55       var const cl_I& b = denominator(x); // b
  56       var cl_signean sign = -(cl_signean)minusp(a); // Vorzeichen
  57       if (!(sign==0)) { a = -a; } // Betrag nehmen, liefert a
  58       var sintC lendiff = (sintC)integer_length(a) // (integer-length a)
  59                           - (sintC)integer_length(b); // (integer-length b)
  60       // |lendiff| < intDsize*2^intCsize. Da für LF-Exponenten ein sintL zur
  61       // Verfügung steht, braucht man keinen Test auf Overflow oder Underflow.
  62       var uintC difflimit = intDsize*len + 1; // 16n+1
  63       var cl_I zaehler;
  64       var cl_I nenner;
  65       if (lendiff > (sintC)difflimit)
  66         // 0 <= k-m-16n-1 < k < intDsize*2^intCsize
  67         { nenner = ash(b,(uintC)(lendiff - difflimit));
  68           zaehler = a;
  69         }
  70         else
  71         // 0 < -k+m+16n+1 <= m+1 + 16n < intDsize*2^intCsize + intDsize*2^intCsize
  72         { zaehler = ash(a,(uintC)(difflimit - lendiff)); // (ash a -k+m+16n+1)
  73           nenner = b; // b
  74         }
  75       // Division zaehler/nenner durchführen:
  76       var cl_I_div_t q_r = cl_divide(zaehler,nenner);
  77       var cl_I& q = q_r.quotient;
  78       var cl_I& r = q_r.remainder;
  79       // 2^16n <= q < 2^(16n+2), also ist q Bignum mit n+1 Digits.
  80       var Lfloat y = allocate_lfloat(len,lendiff+LF_exp_mid,sign); // neues Long-Float
  81       var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
  82       {var uintD* q_MSDptr = arrayMSDptr(TheBignum(q)->data,len+1);
  83        if (mspref(q_MSDptr,0) == 1) // erstes Digit =1 oder =2,3 ?
  84          // 2^16n <= q < 2^(16n+1), also 2^(k-m-1) < a/b < 2^(k-m).
  85          { // Mantisse mit einer Schiebeschleife um 1 Bit nach rechts füllen:
  86            var uintD rounding_bit =
  87              shiftrightcopy_loop_msp(q_MSDptr mspop 1,y_mantMSDptr,len,1,1);
  88            if ( (rounding_bit == 0) // herausgeschobenes Bit =0 -> abrunden
  89                 || ( eq(r,0) // =1 und Rest r > 0 -> aufrunden
  90                      // round-to-even
  91                      && ((mspref(y_mantMSDptr,len-1) & bit(0)) ==0)
  92               )    )
  93              goto ab; // abrunden
  94              else
  95              goto auf; // aufrunden
  96          }
  97          else
  98          // 2^(16n+1) <= q < 2^(16n+2), also 2^(k-m) < a/b < 2^(k-m+1).
  99          { // Mantisse mit einer Schiebeschleife um 2 Bit nach rechts füllen:
 100            var uintD rounding_bits =
 101              shiftrightcopy_loop_msp(q_MSDptr mspop 1,y_mantMSDptr,len,2,mspref(q_MSDptr,0));
 102            (TheLfloat(y)->expo)++; // Exponenten incrementieren auf k-m+1
 103            if ( ((sintD)rounding_bits >= 0) // herausgeschobenes Bit =0 -> abrunden
 104                 || ( ((rounding_bits & bit(intDsize-2)) ==0) // =1 und nächstes Bit =1 oder Rest r > 0 -> aufrunden
 105                      && eq(r,0)
 106                      // round-to-even
 107                      && ((mspref(y_mantMSDptr,len-1) & bit(0)) ==0)
 108               )    )
 109              goto ab; // abrunden
 110              else
 111              goto auf; // aufrunden
 112          }
 113       }
 114       auf: // aufrunden
 115         { if ( inc_loop_lsp(y_mantMSDptr mspop len,len) )
 116             // Übertrag durchs Aufrunden
 117             { mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
 118               (TheLfloat(y)->expo)++; // Exponenten incrementieren
 119         }   }
 120       ab: // abrunden
 121       return y;
 122 }}
 123
 124 // Timings on an i486 33 MHz, running Linux, in 0.01 sec.
 125 // First timing:  cl_I_to_LF(numerator,len)/cl_I_to_LF(denominator,len)
 126 // Second timing: cl_RA_to_LF(x,len)
 127 // with len = 100.
 128 //     num_length    50          70          100         200         500
 129 // den_length
 130 //
 131 //         50     1.86 0.97   1.84 0.97   1.85 0.96   1.86 1.86   1.85 7.14
 132 //
 133 //         70     1.86 1.33   1.85 1.31   1.85 1.32   1.84 1.84   1.85 7.13
 134 //
 135 //        100     1.85 1.85   1.86 1.85   1.85 1.84   1.84 1.84   1.86 7.13
 136 //
 137 //        200     1.85 3.61   1.84 3.61   1.85 3.59   1.85 3.59   1.87 7.12
 138 //
 139 //        500     1.84 7.44   1.84 7.55   1.85 7.56   1.84 7.66   1.86 7.63
 140 //
 141 // We see that cl_RA_to_LF is faster only if
 142 //            num_length < 2*len && den_length < len
 143 // whereas cl_I_to_LF(numerator,len)/cl_I_to_LF(denominator,len) is faster if
 144 //            num_length > 2*len || den_length > len
 145
 146 }  // namespace cln