/cln-1.3.2/src/float/lfloat/elem/cl_LF_from_RA.cc

  1// cl_RA_to_LF().
  2
  3// General includes.
  4#include "base/cl_sysdep.h"
  5
  6// Specification.
  7#include "float/lfloat/cl_LF.h"
  8
  9
 10// Implementation.
 11
 12#include "float/lfloat/cl_LF_impl.h"
 13#include "rational/cl_RA.h"
 14#include "cln/integer.h"
 15#include "integer/cl_I.h"
 16#include "float/cl_F.h"
 17
 18namespace cln {
 19
 20const 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