src/float/transcendental/cl_F_cos.cc

   1 // cos().
   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 cos (const cl_F& x)
  21 {
  22 // Methode:
  23 // Genauigkeit erhöhen,
  24 // (q,r) := (round x (float pi x)), so daß |r|<=pi/2.
  25 // e := Exponent aus (decode-float r), d := (float-digits r)
  26 // Bei r=0.0 oder e<=-d/2 liefere 1.0
  27 //   (denn bei e<=-d/2 ist r^2/2 < 2^(-d)/2 = 2^(-d-1), also
  28 //   1 >= cos(r) > 1-r^2/2 > 1-2^(-d-1),
  29 //   also ist cos(r), auf d Bits gerundet, gleich 1.0).
  30 // Sonst s := r/2 = (scale-float r -1),
  31 //   (sin(s)/s)^2 errechnen, cos(r) = 1-r*s*(sin(s)/s)^2 errechnen.
  32 // Falls q ungerade: Vorzeichenwechsel.
  33
  34         // Rechengenauigkeit erhöhen und durch pi dividieren:
  35         var cl_F cos_r;
  36         if (longfloatp(x)) {
  37                 DeclareType(cl_LF,x);
  38                 if (TheLfloat(x)->len >= 2850) {
  39                         var cl_F_div_t q_r = cl_round_pi2(extend(x,TheLfloat(x)->len+1));
  40                         var cl_I& q = q_r.quotient;
  41                         var cl_LF r = The(cl_LF)(q_r.remainder);
  42                         var cl_LF_cos_sin_t trig = cl_cossin_ratseries(r);
  43                         switch (cl_I_to_UL(logand(q,3))) { // q mod 4
  44                                 case 0: return cl_float(trig.cos,x);
  45                                 case 1: return -cl_float(trig.sin,x);
  46                                 case 2: return -cl_float(trig.cos,x);
  47                                 case 3: return cl_float(trig.sin,x);
  48                                 default: NOTREACHED
  49                         }
  50                 } else {
  51                         var cl_F_div_t q_r = cl_round_pi(cl_F_extendsqrt(x));
  52                         var cl_I& q = q_r.quotient;
  53                         var cl_LF r = The(cl_LF)(q_r.remainder);
  54                         if (zerop(r) || (float_exponent(r) <= (-(sintL)float_digits(r))>>1))
  55                                 cos_r = cl_float(1,x); // (cos r) = 1.0
  56                         else {
  57                                 var cl_LF s = scale_float(r,-1); // s := r/2
  58                                 cos_r = cl_float(1-scale_float(sinx_naive(s),1),x); // cos(2s) = 1-2*sin(s)^2
  59                         }
  60                         if (oddp(q))
  61                                 return -cos_r; // q ungerade -> mal -1
  62                         else
  63                                 return cos_r;
  64                 }
  65         } else {
  66                 var cl_F_div_t q_r = cl_round_pi(cl_F_extendsqrt(x));
  67                 var cl_I& q = q_r.quotient;
  68                 var cl_F& r = q_r.remainder;
  69                 if (zerop(r) || (float_exponent(r) <= (-(sintL)float_digits(r))>>1))
  70                         cos_r = cl_float(1,x); // (cos r) = 1.0
  71                 else {
  72                         var cl_F s = scale_float(r,-1); // s := r/2
  73                         cos_r = cl_float(1 - r * s * sinxbyx_naive(s),x);
  74                 }
  75                 if (oddp(q))
  76                         return -cos_r; // q ungerade -> mal -1
  77                 else
  78                         return cos_r;
  79         }
  80 }
  81
  82 // Timings of the two algorithms, on an i486 33 MHz, running Linux,
  83 // applied to x = sqrt(2)-1 = 0.414...
  84 //   N      naive  ratseries
  85 //   10     0.009   0.049
  86 //   25     0.033   0.137
  87 //   50     0.11    0.37
  88 //  100     0.41    1.15
  89 //  250     2.7     5.5
  90 //  500    11.1    19.4
  91 // 1000    46      64
  92 // 2500   239     260
  93 // ==> ratseries faster for N >= 2850.
  94
  95 }  // namespace cln