src/float/transcendental/cl_F_sinx.cc

   1 // sinxbyx(), sinx().
   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_exponent.cc"
  23
  24 namespace cln {
  25
  26 // sinxbyx is mainly for cl_SF, cl_FF, cl_DF, where we want to avoid underflow.
  27
  28 const cl_F sinxbyx_naive (const cl_F& x)
  29 {
  30 // Methode:
  31 // e := Exponent aus (decode-float x), d := (float-digits x)
  32 // Bei x=0.0 oder e<=-d/2 liefere 1.0
  33 //   (denn bei e<=-d/2 ist x^2/6 < x^2/4 < 2^(-d)/4 = 2^(-d-2), also
  34 //   1 >= sin(x)/x > 1-x^2/6 > 1-2^(-d-2), also 1 >= (sin(x)/x)^2 > 1-2^(-d-1),
  35 //   also ist (sin(x)/x)^2, auf d Bits gerundet, gleich 1.0).
  36 // Bei e<=-sqrt(d) verwende die Potenzreihe
  37 //   sin(x)/x = sum(j=0..inf,(-x^2)^j/(2j+1)!):
  38 //   a:=-x^2, b:=1, i:=1, sum:=0,
  39 //   while (/= sum (setq sum (+ sum b))) do b:=b*a/((i+1)*(i+2)), i:=i+2.
  40 //   Ergebnis sum^2.
  41 // Sonst setze y := x/2 = (scale-float x -1),
  42 //   berechne rekursiv z:=(sin(y)/y)^2 und liefere z*(1-y^2*z).
  43 // [Die Grenze sqrt(d) ergibt sich so:
  44 //  Man braucht bei der Potenzreihe mit x=2^-k etwa j Glieder, mit
  45 //  k*j*ln 2 + j*(ln j - 1) = d, und der Aufwand beträgt etwa 2.8*(j/2)
  46 //  Multiplikationen von d-Bit-Zahlen. Bei Halbierungen bis x=2^-k ist der
  47 //  Gesamtaufwand etwa 2*(k+e)+1.4*j(k). Dieses minimieren nach k: Soll sein
  48 //  -1.4 = d/dk j(k) = (d/dj k(j))^-1 = - j^2/(d+j)*ln 2, also j^2=2(d+j),
  49 //  grob j=sqrt(2d) und damit k=sqrt(d).]
  50 // Aufwand: asymptotisch d^2.5 .
  51
  52         if (zerop(x))
  53                 return cl_float(1,x);
  54         var uintC d = float_digits(x);
  55         var sintE e = float_exponent(x);
  56         if (e <= (-(sintC)d)>>1) // e <= (-d)/2 <==> e <= -ceiling(d/2) ?
  57                 return cl_float(1,x); // ja -> 1.0 als Ergebnis
  58  {      Mutable(cl_F,x);
  59         // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
  60         // angewandt werden. limit_slope = 1.0 ist sehr schlecht (ca. 30%
  61         // zu schlecht). Gute Werte bei N limbs:
  62         //    N     limit_slope
  63         //    5       0.15-0.30
  64         //   10          0.25
  65         //   25       0.25-0.35
  66         //   50       0.35-0.40
  67         //  100       0.40-0.45
  68         //  200       0.40-0.45
  69         // Wähle limit_slope = 13/32 = 0.4.
  70         var sintL e_limit = -1-floor(isqrtC(d)*13,32); // -1-floor(sqrt(d))
  71         if (e > e_limit) {
  72                 // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
  73                 x = scale_float(x,e_limit-e);
  74                 // Neuer Exponent = e_limit.
  75         }
  76         var cl_F x2 = square(x);        // x^2
  77         // Potenzreihe anwenden:
  78         var cl_F a = - x2; // a := -x^2
  79         var int i = 1;
  80         var cl_F b = cl_float(1,x); // b := (float 1 x)
  81         var cl_F sum = cl_float(0,x); // sum := (float 0 x)
  82         loop {
  83                 var cl_F new_sum = sum + b;
  84                 if (new_sum == sum) // = sum ?
  85                         break; // ja -> Potenzreihe abbrechen
  86                 sum = new_sum;
  87                 b = (b*a)/(cl_I)((i+1)*(i+2));
  88                 i = i+2;
  89         }
  90         var cl_F z = square(sum); // sum^2 als Ergebnis
  91         while (e > e_limit) {
  92                 z = z - x2 * square(z);
  93                 x2 = scale_float(x2,2); // x^2 := x^2*4
  94                 e--;
  95         }
  96         return z;
  97 }}
  98 // Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
  99
 100 const cl_LF sinx_naive (const cl_LF& x)
 101 {
 102 // Methode:
 103 // e := Exponent aus (decode-float x), d := (float-digits x)
 104 // Bei x=0.0 oder e<=-d/2 liefere x
 105 //   (denn bei e<=-d/2 ist x^2/6 < x^2/4 < 2^(-d)/4 = 2^(-d-2), also
 106 //   1 >= sin(x)/x > 1-x^2/6 > 1-2^(-d-2), also ist sin(x)^2, auf d Bits
 107 //   gerundet, gleich x).
 108 // Bei e<=-sqrt(d) verwende die Potenzreihe
 109 //   sin(x) = sum(j=0..inf,x*(-x^2)^j/(2j+1)!):
 110 //   a:=-x^2, b:=x, i:=1, sum:=0,
 111 //   while (/= sum (setq sum (+ sum b))) do b:=b*a/((i+1)*(i+2)), i:=i+2.
 112 //   Ergebnis sum^2.
 113 // Sonst setze y := x/2 = (scale-float x -1),
 114 //   berechne rekursiv z:=sin(y)^2 und liefere 4*z*(1-z) = 1-(1-2*z)^2.
 115 // [Die Grenze sqrt(d) ergibt sich so:
 116 //  Man braucht bei der Potenzreihe mit x=2^-k etwa j Glieder, mit
 117 //  k*j*ln 2 + j*(ln j - 1) = d, und der Aufwand beträgt etwa 2.8*(j/2)
 118 //  Multiplikationen von d-Bit-Zahlen. Bei Halbierungen bis x=2^-k ist der
 119 //  Gesamtaufwand etwa 2*(k+e)+1.4*j(k). Dieses minimieren nach k: Soll sein
 120 //  -1.4 = d/dk j(k) = (d/dj k(j))^-1 = - j^2/(d+j)*ln 2, also j^2=2(d+j),
 121 //  grob j=sqrt(2d) und damit k=sqrt(d).]
 122 // Aufwand: asymptotisch d^2.5 .
 123
 124         if (zerop(x))
 125                 return x;
 126         var uintC actuallen = TheLfloat(x)->len;
 127         var uintC d = float_digits(x);
 128         var sintE e = float_exponent(x);
 129         if (e <= (-(sintC)d)>>1) // e <= (-d)/2 <==> e <= -ceiling(d/2) ?
 130                 return square(x); // ja -> x^2 als Ergebnis
 131  {      Mutable(cl_LF,x);
 132         var sintE ee = e;
 133         // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
 134         // angewandt werden. limit_slope = 1.0 ist schlecht (ca. 10% zu
 135         // schlecht). Ein guter Wert für naive1 ist limit_slope = 0.6,
 136         // für naive3 aber limit_slope = 0.5.
 137         var sintL e_limit = -1-floor(isqrtC(d),2); // -1-floor(sqrt(d))
 138         if (e > e_limit) {
 139                 // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
 140                 x = scale_float(x,e_limit-e);
 141                 ee = e_limit; // Neuer Exponent = e_limit.
 142         }
 143         var cl_LF x2 = square(x); // x^2
 144         // Potenzreihe anwenden:
 145         var cl_LF powser_value;
 146         var cl_LF a = - x2; // a := -x^2
 147         var int i = 1;
 148         if (0) {
 149                 // naive1:
 150                 // fixed-point representation
 151                 d = d-ee; // fixed-point representation with d mantissa bits
 152                 var cl_I b = round1(scale_float(x,d)); // b := x
 153                 var cl_I sum = 0; // sum := (float 0 x)
 154                 loop {
 155                         if (b == 0) break;
 156                         sum = sum + b;
 157                         b = round1(round1(The(cl_LF)(b*a)),(cl_I)((i+1)*(i+2)));
 158                         i = i+2;
 159                 }
 160                 powser_value = scale_float(cl_float(sum,x),-(sintC)d);
 161         } else if (actuallen <= 7) { // Break-even-Point before extendsqrt: N<=6
 162                 // naive2:
 163                 // floating-point representation
 164                 var cl_LF b = x; // b := x
 165                 var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
 166                 loop {
 167                         var cl_LF new_sum = sum + b;
 168                         if (new_sum == sum) // = sum ?
 169                                 break; // ja -> Potenzreihe abbrechen
 170                         sum = new_sum;
 171                         b = (b*a)/(cl_I)((i+1)*(i+2));
 172                         i = i+2;
 173                 }
 174                 powser_value = sum;
 175         } else {
 176                 // naive3:
 177                 // floating-point representation with smooth precision reduction
 178                 var cl_LF b = x; // b := x
 179                 var cl_LF eps = scale_float(b,-(sintC)d-10);
 180                 var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
 181                 loop {
 182                         var cl_LF new_sum = sum + LF_to_LF(b,actuallen);
 183                         if (new_sum == sum) // = sum ?
 184                                 break; // ja -> Potenzreihe abbrechen
 185                         sum = new_sum;
 186                         b = cl_LF_shortenwith(b,eps);
 187                         b = (b*a)/(cl_I)((i+1)*(i+2));
 188                         i = i+2;
 189                 }
 190                 powser_value = sum;
 191         }
 192         var cl_LF z = square(powser_value); // sin^2 als Ergebnis
 193         while (e > e_limit) {
 194                 z = cl_float(1,x) - square(cl_float(1,x) - scale_float(z,1)); // z := 1-(1-2*z)^2
 195                 e--;
 196         }
 197         return z;
 198 }}
 199 // Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
 200
 201 // Timings of the three variants, on an i486 33 MHz, running Linux,
 202 // applied to x = sqrt(2)-1 = 0.414...
 203 //   N     naive1  naive2  naive3  ratseries
 204 //    4     0.0064  0.0048  0.0049  0.023
 205 //    6     0.0081  0.0064  0.0065  0.031
 206 //    8     0.0103  0.0085  0.0083  0.038
 207 //   10     0.012   0.011   0.010   0.048
 208 //   25     0.043   0.047   0.035   0.119
 209 //   50     0.15    0.17    0.12    0.37
 210 //  100     0.54    0.67    0.44    1.09
 211 //  250     3.5     4.4     2.8     5.5
 212 //  500    14.7    18.5    11.6    19.4
 213 // 1000    61      78      48      64
 214 // 2500   315     361     243     261
 215 // 2700                   265     270
 216 // 3000                   294     282
 217 // 3500                   339     303
 218 // ==> ratseries faster for N >= 2750.
 219
 220 }  // namespace cln