1 // cos().
2
3 // General includes.
4 #include "base/cl_sysdep.h"
5
6 // Specification.
7 #include "cln/float.h"
8
9
10 // Implementation.
11
12 #include "float/transcendental/cl_F_tran.h"
13 #include "float/cl_F.h"
14 #include "cln/integer.h"
15 #include "cln/lfloat.h"
16 #include "float/lfloat/cl_LF.h"
17
18 namespace cln {
19
cos(const cl_F & x)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) <= (-(sintC)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) <= (-(sintC)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
96