src/float/transcendental/cl_F_sin.cc

   1 // sin().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/float.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_F_tran.h"
  13 #include "cl_F.h"
  14 #include "cln/integer.h"
  15 #include "cln/lfloat.h"
  16 #include "cl_LF.h"
  17
  18 namespace cln {
  19
  20 const cl_F sin (const cl_F& x)
  21 {
  22 // Methode:
  23 // Genauigkeit erhöhen,
  24 // (q,r) := (round x (float pi/2 x)), so daß |r|<=pi/4.
  25 // y:=(sin(r)/r)^2 errechnen.
  26 // Falls q gerade:
  27 //   sin(r) berechnen: r*sqrt(y).
  28 // Falls q ungerade:
  29 //   cos(r) berechnen:
  30 //     e := Exponent aus (decode-float r), d := (float-digits r)
  31 //     Bei r=0.0 oder e<=-d/2 liefere 1.0
  32 //       (denn bei e<=-d/2 ist r^2/2 < 2^(-d)/2 = 2^(-d-1), also
  33 //       1 >= cos(r) > 1-r^2/2 > 1-2^(-d-1),
  34 //       also ist cos(r), auf d Bits gerundet, gleich 1.0).
  35 //     Sonst sqrt(1-r^2*y).
  36 // Falls q == 2,3 mod 4, Vorzeichenwechsel.
  37
  38         // Rechengenauigkeit erhöhen und durch pi/2 dividieren:
  39         var cl_F z;
  40         var cl_I q;
  41         if (longfloatp(x)) {
  42                 DeclareType(cl_LF,x);
  43                 if (TheLfloat(x)->len >= 2750) {
  44                         var cl_F_div_t q_r = cl_round_pi2(extend(x,TheLfloat(x)->len+1));
  45                         q = q_r.quotient;
  46                         var cl_LF r = The(cl_LF)(q_r.remainder);
  47                         var cl_LF_cos_sin_t trig = cl_cossin_ratseries(r);
  48                         if (evenp(q))
  49                                 z = cl_float(trig.sin,x);
  50                         else
  51                                 z = cl_float(trig.cos,x);
  52                 } else {
  53                         var cl_F_div_t q_r = cl_round_pi2(cl_F_extendsqrt(x));
  54                         q = q_r.quotient;
  55                         var cl_LF r = The(cl_LF)(q_r.remainder);
  56                         var cl_LF y = sinx_naive(r); // y := sin(r)^2
  57                         if (evenp(q)) {
  58                                 // sin(r) berechnen:
  59                                 z = cl_float(sqrt(y),x);
  60                                 if (minusp(r))
  61                                         z = -z;
  62                         } else {
  63                                 // cos(r) berechnen:
  64                                 if (zerop(r) || (float_exponent(r) <= (-(sintC)float_digits(r))>>1))
  65                                         z = cl_float(1,x); // cos(r) = 1.0
  66                                 else
  67                                         z = cl_float(sqrt(1 - y),x); // sqrt(1-y)
  68                         }
  69                 }
  70         } else {
  71                 var cl_F_div_t q_r = cl_round_pi2(cl_F_extendsqrt(x));
  72                 q = q_r.quotient;
  73                 var cl_F& r = q_r.remainder;
  74                 var cl_F y = sinxbyx_naive(r); // y := (sin(r)/r)^2
  75                 if (evenp(q)) {
  76                         // sin(r) berechnen:
  77                         z = cl_float(r*sqrt(y),x);
  78                 } else {
  79                         // cos(r) berechnen:
  80                         if (zerop(r) || (float_exponent(r) <= (-(sintC)float_digits(r))>>1))
  81                                 z = cl_float(1,x); // cos(r) = 1.0
  82                         else
  83                                 z = cl_float(sqrt(1 - square(r)*y),x); // sqrt(1-r^2*y)
  84                 }
  85         }
  86         // evtl. Vorzeichenwechsel:
  87         if (cl_I_to_UL(logand(q,2))==0)
  88                 return z;
  89         else
  90                 return -z;
  91 }
  92
  93 // Timings of the two algorithms, on an i486 33 MHz, running Linux,
  94 // applied to x = sqrt(2)-1 = 0.414...
  95 //   N      naive  ratseries
  96 //   10     0.010   0.048
  97 //   25     0.035   0.119
  98 //   50     0.12    0.37
  99 //  100     0.44    1.09
 100 //  250     2.8     5.5
 101 //  500    11.6    19.4
 102 // 1000    48      64
 103 // 2500   243     261
 104 // ==> ratseries faster for N >= 2750.
 105
 106 }  // namespace cln