src/float/transcendental/cl_F_lnx.cc

   1 // lnx().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_F_tran.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/float.h"
  13 #include "cl_low.h"
  14 #include "cl_F.h"
  15 #include "cln/lfloat.h"
  16 #include "cl_LF.h"
  17 #include "cln/integer.h"
  18
  19 #undef MAYBE_INLINE
  20 #define MAYBE_INLINE inline
  21 #include "cl_LF_zerop.cc"
  22 #include "cl_LF_minusp.cc"
  23 #include "cl_LF_exponent.cc"
  24
  25 namespace cln {
  26
  27 // cl_F lnx_naive (const cl_F& x)
  28 // cl_LF lnx_naive (const cl_LF& x)
  29 //
  30 // Methode:
  31 // y:=x-1, e := Exponent aus (decode-float y), d := (float-digits y)
  32 // Bei y=0.0 oder e<=-d liefere y
  33 //   (denn bei e<=-d ist y/2 < 2^(-d)/2 = 2^(-d-1), also
  34 //   0 <= y - ln(x) < y^2/2 < 2^(-d-1)*y
  35 //   also ist ln(x)/y, auf d Bits gerundet, gleich y).
  36 // Bei e<=-sqrt(d) verwende die Potenzreihe
  37 //   ln(x) = sum(j=0..inf,(-1)^j*y^(j+1)/(j+1)):
  38 //   a:=-y, b:=y, i:=1, sum:=0,
  39 //   while (/= sum (setq sum (+ sum (/ b i)))) do i:=i+1, b:=b*a.
  40 //   Ergebnis sum.
  41 // Sonst setze y := sqrt(x), berechne rekursiv z:=ln(y)
  42 //   und liefere 2*z = (scale-float z 1).
  43 // Aufwand: asymptotisch d^0.5*M(d) = d^2.5 .
  44
  45 const cl_LF lnx_naive (const cl_LF& x)
  46 {
  47         var cl_LF y = x-cl_float(1,x);
  48         if (zerop(y)) // y=0.0 -> y als Ergebnis
  49                 return y;
  50         var uintL actuallen = TheLfloat(x)->len;
  51         var uintL d = float_digits(x);
  52         var sintL e = float_exponent(y);
  53         if (e <= -(sintL)d) // e <= -d ?
  54                 return y; // ja -> y als Ergebnis
  55  {      Mutable(cl_LF,x);
  56         var uintL k = 0; // Rekursionszähler k:=0
  57         // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
  58         // angewandt werden.
  59         // Wähle für ln(1+y), naive1: limit_slope = 1.0,
  60         //       für ln(1+y), naive2: limit_slope = 11/16 = 0.7,
  61         //       für atanh(z), naive1: limit_slope = 0.6,
  62         //       für atanh(z), naive1: limit_slope = 0.5.
  63         var sintL e_limit = -1-floor(isqrt(d),2); // -1-floor(sqrt(d))
  64         while (e > e_limit) {
  65                 // e > -1-floor(sqrt(d)) -> muß |y| verkleinern.
  66                 x = sqrt(x); // x := (sqrt x)
  67                 y = x-cl_float(1,x); // y := (- x 1) und
  68                 e = float_exponent(y); // e neu berechnen
  69                 k = k+1; // k:=k+1
  70         }
  71         if (0) {
  72                 // Potenzreihe ln(1+y) anwenden:
  73                 var int i = 1;
  74                 var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
  75                 var cl_LF a = -y;
  76                 var cl_LF b = y;
  77                 if (0) {
  78                         // naive1:
  79                         // floating-point representation
  80                         loop {
  81                                 var cl_LF new_sum = sum + b/(cl_I)i; // (+ sum (/ b i))
  82                                 if (new_sum == sum) // = sum ?
  83                                         break; // ja -> Potenzreihe abbrechen
  84                                 sum = new_sum;
  85                                 b = b*a;
  86                                 i = i+1;
  87                         }
  88                 } else {
  89                         // naive2:
  90                         // floating-point representation with smooth precision reduction
  91                         var cl_LF eps = scale_float(b,-(sintL)d-10);
  92                         loop {
  93                                 var cl_LF new_sum = sum + LF_to_LF(b/(cl_I)i,actuallen); // (+ sum (/ b i))
  94                                 if (new_sum == sum) // = sum ?
  95                                         break; // ja -> Potenzreihe abbrechen
  96                                 sum = new_sum;
  97                                 b = cl_LF_shortenwith(b,eps);
  98                                 b = b*a;
  99                                 i = i+1;
 100                         }
 101                 }
 102                 return scale_float(sum,k); // sum als Ergebnis, wegen Rekursion noch mal 2^k
 103         } else {
 104                 var cl_LF z = y / (x+cl_float(1,x));
 105                 // Potenzreihe atanh(z) anwenden:
 106                 var int i = 1;
 107                 var cl_LF a = square(z); // a = x^2
 108                 var cl_LF b = cl_float(1,x); // b := (float 1 x)
 109                 var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
 110                 if (0) {
 111                         // naive1:
 112                         // floating-point representation
 113                         loop {
 114                                 var cl_LF new_sum = sum + b / (cl_I)i; // (+ sum (/ b i))
 115                                 if (new_sum == sum) // = sum ?
 116                                         break; // ja -> Potenzreihe abbrechen
 117                                 sum = new_sum;
 118                                 b = b*a;
 119                                 i = i+2;
 120                         }
 121                 } else {
 122                         // naive2:
 123                         // floating-point representation with smooth precision reduction
 124                         var cl_LF eps = scale_float(b,-(sintL)d-10);
 125                         loop {
 126                                 var cl_LF new_sum = sum + LF_to_LF(b/(cl_I)i,actuallen); // (+ sum (/ b i))
 127                                 if (new_sum == sum) // = sum ?
 128                                         break; // ja -> Potenzreihe abbrechen
 129                                 sum = new_sum;
 130                                 b = cl_LF_shortenwith(b,eps);
 131                                 b = b*a;
 132                                 i = i+2;
 133                         }
 134                 }
 135                 return scale_float(sum*z,k+1); // 2*sum*z als Ergebnis, wegen Rekursion noch mal 2^k
 136         }
 137 }}
 138 // Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
 139
 140 const cl_F lnx_naive (const cl_F& x)
 141 {
 142         if (longfloatp(x)) {
 143                 DeclareType(cl_LF,x);
 144                 return lnx_naive(x);
 145         }
 146         var cl_F y = x-cl_float(1,x);
 147         if (zerop(y)) // y=0.0 -> y als Ergebnis
 148                 return y;
 149         var uintL d = float_digits(x);
 150         var sintL e = float_exponent(y);
 151         if (e <= -(sintL)d) // e <= -d ?
 152                 return y; // ja -> y als Ergebnis
 153  {      Mutable(cl_F,x);
 154         var uintL k = 0; // Rekursionszähler k:=0
 155         // Bei e <= -1-floor(sqrt(d)) kann die Potenzreihe angewandt werden.
 156         var sintL e_limit = -1-isqrt(d); // -1-floor(sqrt(d))
 157         while (e > e_limit) {
 158                 // e > -1-floor(sqrt(d)) -> muß |y| verkleinern.
 159                 x = sqrt(x); // x := (sqrt x)
 160                 y = x-cl_float(1,x); // y := (- x 1) und
 161                 e = float_exponent(y); // e neu berechnen
 162                 k = k+1; // k:=k+1
 163         }
 164         // Potenzreihe anwenden:
 165         var int i = 1;
 166         var cl_F sum = cl_float(0,x); // sum := (float 0 x)
 167         var cl_F a = -y;
 168         var cl_F b = y;
 169         loop {
 170                 var cl_F new_sum = sum + b/(cl_I)i; // (+ sum (/ b i))
 171                 if (new_sum == sum) // = sum ?
 172                         break; // ja -> Potenzreihe abbrechen
 173                 sum = new_sum;
 174                 b = b*a;
 175                 i = i+1;
 176         }
 177         return scale_float(sum,k); // sum als Ergebnis, wegen Rekursion noch mal 2^k
 178 }}
 179 // Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
 180
 181 const cl_LF lnx_ratseries (const cl_LF& x)
 182 {
 183         // Method:
 184         // Based on the same ideas as expx_ratseries.
 185         // y := 0.
 186         // Loop
 187         //   [x*exp(y) is invariant]
 188         //   x' := x-1. If x' = 0, terminate the loop.
 189         //   Choose approximation y' of log(x) = log(1+x'):
 190         //     If |x'| >= 1/2, set y' = 1/2 * sign(x').
 191         //     If |x'| < 2^-n with n maximal, set
 192         //       y' = truncate(x'*2^(2n))/2^(2n).
 193         //   Set y := y + y' and x := x*exp(-y').
 194         var uintC len = TheLfloat(x)->len;
 195  {      Mutable(cl_LF,x);
 196         var cl_LF y = cl_I_to_LF(0,len);
 197         loop {
 198                 var cl_LF x1 = x + cl_I_to_LF(-1,len);
 199                 var cl_idecoded_float x1_ = integer_decode_float(x1);
 200                 // x1 = (-1)^sign * 2^exponent * mantissa
 201                 if (zerop(x1_.mantissa))
 202                         break;
 203                 var uintL lm = integer_length(x1_.mantissa);
 204                 var uintL me = cl_I_to_UL(- x1_.exponent);
 205                 var cl_I p;
 206                 var uintL lq;
 207                 var cl_boolean last_step = cl_false;
 208                 if (lm >= me) { // |x'| >= 1/2 ?
 209                         p = x1_.sign; // 1 or -1
 210                         lq = 1;
 211                 } else {
 212                         var uintL n = me - lm; // |x'| < 2^-n with n maximal
 213                         // Set p to the first n bits of |x'|:
 214                         if (lm > n) {
 215                                 p = x1_.mantissa >> (lm - n);
 216                                 lq = 2*n;
 217                         } else {
 218                                 p = x1_.mantissa;
 219                                 lq = lm + n;
 220                         }
 221                         if (minusp(x1_.sign)) { p = -p; }
 222                         // If 2*n >= lm = intDsize*len, then within our
 223                         // precision exp(-y') = 1-y', (because |y'^2| < 2^-lm),
 224                         // and we know a priori that the iteration will stop
 225                         // after the next big multiplication. This saves one
 226                         // big multiplication at the end.
 227                         if (2*n >= lm)
 228                                 last_step = cl_true;
 229                 }
 230                 y = y + scale_float(cl_I_to_LF(p,len),-(sintL)lq);
 231                 if (last_step)
 232                         break;
 233                 x = x * cl_exp_aux(-p,lq,len);
 234         }
 235         return y;
 236 }}
 237 // Bit complexity (N = length(x)): O(log(N)^2*M(N)).
 238
 239 // Timings of the above algorithms, on an i486 33 MHz, running Linux,
 240 // applied to x = sqrt(sqrt(2)) = 1.189...
 241 //   N      ln(1+y) ln(1+y) atanh z atanh z   exp
 242 //          naive1  naive2  naive1  naive2  ratseries
 243 //   10     0.019   0.016   0.013   0.012   0.036
 244 //   25     0.077   0.056   0.057   0.040   0.087
 245 //   50     0.30    0.21    0.23    0.15    0.21
 246 //  100     1.24    0.81    0.92    0.59    0.61
 247 //  250     8.8     5.8     6.3     4.3     2.77
 248 //  500    43.9    28.8    29.7    21.0     9.8
 249 // 1000   223     149     144     107      30
 250 // ==> ratseries faster for N >= 110. (N = length before extended by the caller.)
 251
 252 }  // namespace cln