12 #include "cl_F_tran.h"
14 #include "cln/integer.h"
15 #include "cln/lfloat.h"
20 const cl_F cos (const cl_F& x)
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.
34 // Rechengenauigkeit erhöhen und durch pi dividieren:
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);
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
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
61 return -cos_r; // q ungerade -> mal -1
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
72 var cl_F s = scale_float(r,-1); // s := r/2
73 cos_r = cl_float(1 - r * s * sinxbyx_naive(s),x);
76 return -cos_r; // q ungerade -> mal -1
82 // Timings of the two algorithms, on an i486 33 MHz, running Linux,
83 // applied to x = sqrt(2)-1 = 0.414...
93 // ==> ratseries faster for N >= 2850.