src/float/transcendental/cl_F_cossin.cc

   1 // cos_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 cos_sin_t cos_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 // cos(r) berechnen:
  27 //   e := Exponent aus (decode-float r), d := (float-digits r)
  28 //   Bei r=0.0 oder e<=-d/2 liefere 1.0
  29 //     (denn bei e<=-d/2 ist r^2/2 < 2^(-d)/2 = 2^(-d-1), also
  30 //     1 >= cos(r) > 1-r^2/2 > 1-2^(-d-1),
  31 //     also ist cos(r), auf d Bits gerundet, gleich 1.0).
  32 //   Sonst sqrt(1-r^2*y).
  33 // sin(r) berechnen: r*sqrt(y).
  34 // Genauigkeit wieder verringern.
  35 // Falls q = 0 mod 4: (cos(r), sin(r))
  36 // Falls q = 1 mod 4: (-sin(r), cos(r))
  37 // Falls q = 2 mod 4: (-cos(r), -sin(r))
  38 // Falls q = 3 mod 4: (sin(r), -cos(r))
  39
  40         // Rechengenauigkeit erhöhen und durch pi/2 dividieren:
  41         var cl_F cos_r;
  42         var cl_F sin_r;
  43         var cl_I q;
  44         if (longfloatp(x)) {
  45                 DeclareType(cl_LF,x);
  46                 if (TheLfloat(x)->len >= 2710) {
  47                         var cl_F_div_t q_r = cl_round_pi2(extend(x,TheLfloat(x)->len+1));
  48                         q = q_r.quotient;
  49                         var cl_LF r = The(cl_LF)(q_r.remainder);
  50                         var cl_LF_cos_sin_t trig = cl_cossin_ratseries(r);
  51                         cos_r = cl_float(trig.cos,x);
  52                         sin_r = cl_float(trig.sin,x);
  53                 } else {
  54                         var cl_F_div_t q_r = cl_round_pi2(cl_F_extendsqrt(x));
  55                         q = q_r.quotient;
  56                         var cl_LF r = The(cl_LF)(q_r.remainder);
  57                         var cl_LF y = sinx_naive(r); // y := sin(r)^2
  58                         // erste Komponente cos(r) berechnen:
  59                         if (zerop(r) || (float_exponent(r) <= (-(sintL)float_digits(r))>>1))
  60                                 cos_r = cl_float(1,x); // cos(r) = 1.0
  61                         else
  62                                 cos_r = cl_float(sqrt(1-y),x); // cos(r) = sqrt(1-y)
  63                         // zweite Komponente sin(r) berechnen:
  64                         sin_r = cl_float(sqrt(y),x);
  65                         if (minusp(r))
  66                                 sin_r = - sin_r;
  67                 }
  68         } else {
  69                 var cl_F_div_t q_r = cl_round_pi2(cl_F_extendsqrt(x));
  70                 q = q_r.quotient;
  71                 var cl_F& r = q_r.remainder;
  72                 var cl_F y = sinxbyx_naive(r); // y := (sin(r)/r)^2
  73                 // erste Komponente cos(r) berechnen:
  74                 if (zerop(r) || (float_exponent(r) <= (-(sintL)float_digits(r))>>1))
  75                         cos_r = cl_float(1,x); // cos(r) = 1.0
  76                 else
  77                         cos_r = cl_float(sqrt(1 - square(r)*y),x); // sqrt(1-r^2*y)
  78                 // zweite Komponente sin(r) berechnen:
  79                 sin_r = cl_float(r*sqrt(y),x);
  80         }
  81         // evtl. Vorzeichenwechsel oder Vertauschen:
  82         switch (cl_I_to_UL(logand(q,3))) { // q mod 4
  83                 case 0: return cos_sin_t(cos_r,sin_r);
  84                 case 1: return cos_sin_t(-sin_r,cos_r);
  85                 case 2: return cos_sin_t(-cos_r,-sin_r);
  86                 case 3: return cos_sin_t(sin_r,-cos_r);
  87                 default: NOTREACHED
  88         }
  89 }
  90
  91 }  // namespace cln