src/float/transcendental/cl_F_cossin.cc

   1 // cl_cos_sin().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_float.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_F_tran.h"
  13 #include "cl_F.h"
  14 #include "cl_integer.h"
  15 #include "cl_lfloat.h"
  16 #include "cl_LF.h"
  17
  18 const cl_cos_sin_t cl_cos_sin (const cl_F& x)
  19 {
  20 // Methode:
  21 // Genauigkeit erhöhen,
  22 // (q,r) := (round x (float pi/2 x)), so daß |r|<=pi/4.
  23 // y:=(sin(r)/r)^2 errechnen.
  24 // cos(r) berechnen:
  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 sqrt(1-r^2*y).
  31 // sin(r) berechnen: r*sqrt(y).
  32 // Genauigkeit wieder verringern.
  33 // Falls q = 0 mod 4: (cos(r), sin(r))
  34 // Falls q = 1 mod 4: (-sin(r), cos(r))
  35 // Falls q = 2 mod 4: (-cos(r), -sin(r))
  36 // Falls q = 3 mod 4: (sin(r), -cos(r))
  37
  38         // Rechengenauigkeit erhöhen und durch pi/2 dividieren:
  39         var cl_F cos_r;
  40         var cl_F sin_r;
  41         var cl_I q;
  42         if (longfloatp(x)) {
  43                 DeclareType(cl_LF,x);
  44                 if (TheLfloat(x)->len >= 2710) {
  45                         var cl_F_div_t q_r = cl_round_pi2(extend(x,TheLfloat(x)->len+1));
  46                         q = q_r.quotient;
  47                         var cl_LF r = The(cl_LF)(q_r.remainder);
  48                         var cl_LF_cos_sin_t trig = cl_cossin_ratseries(r);
  49                         cos_r = cl_float(trig.cos,x);
  50                         sin_r = cl_float(trig.sin,x);
  51                 } else {
  52                         var cl_F_div_t q_r = cl_round_pi2(cl_F_extendsqrt(x));
  53                         q = q_r.quotient;
  54                         var cl_LF r = The(cl_LF)(q_r.remainder);
  55                         var cl_LF y = sinx_naive(r); // y := sin(r)^2
  56                         // erste Komponente cos(r) berechnen:
  57                         if (zerop(r) || (float_exponent(r) <= (-(sintL)float_digits(r))>>1))
  58                                 cos_r = cl_float(1,x); // cos(r) = 1.0
  59                         else
  60                                 cos_r = cl_float(sqrt(1-y),x); // cos(r) = sqrt(1-y)
  61                         // zweite Komponente sin(r) berechnen:
  62                         sin_r = cl_float(sqrt(y),x);
  63                         if (minusp(r))
  64                                 sin_r = - sin_r;
  65                 }
  66         } else {
  67                 var cl_F_div_t q_r = cl_round_pi2(cl_F_extendsqrt(x));
  68                 q = q_r.quotient;
  69                 var cl_F& r = q_r.remainder;
  70                 var cl_F y = sinxbyx_naive(r); // y := (sin(r)/r)^2
  71                 // erste Komponente cos(r) berechnen:
  72                 if (zerop(r) || (float_exponent(r) <= (-(sintL)float_digits(r))>>1))
  73                         cos_r = cl_float(1,x); // cos(r) = 1.0
  74                 else
  75                         cos_r = cl_float(sqrt(1 - square(r)*y),x); // sqrt(1-r^2*y)
  76                 // zweite Komponente sin(r) berechnen:
  77                 sin_r = cl_float(r*sqrt(y),x);
  78         }
  79         // evtl. Vorzeichenwechsel oder Vertauschen:
  80         switch (cl_I_to_UL(logand(q,3))) { // q mod 4
  81                 case 0: return cl_cos_sin_t(cos_r,sin_r);
  82                 case 1: return cl_cos_sin_t(-sin_r,cos_r);
  83                 case 2: return cl_cos_sin_t(-cos_r,-sin_r);
  84                 case 3: return cl_cos_sin_t(sin_r,-cos_r);
  85                 default: NOTREACHED
  86         }
  87 }