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