12 #include "cl_F_tran.h"
14 #include "cln/integer.h"
15 #include "cln/lfloat.h"
20 const cl_F sin (const cl_F& x)
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.
27 // sin(r) berechnen: r*sqrt(y).
30 // e := Exponent aus (decode-float r), d := (float-digits r)
31 // Bei r=0.0 oder e<=-d/2 liefere 1.0
32 // (denn bei e<=-d/2 ist r^2/2 < 2^(-d)/2 = 2^(-d-1), also
33 // 1 >= cos(r) > 1-r^2/2 > 1-2^(-d-1),
34 // also ist cos(r), auf d Bits gerundet, gleich 1.0).
35 // Sonst sqrt(1-r^2*y).
36 // Falls q == 2,3 mod 4, Vorzeichenwechsel.
38 // Rechengenauigkeit erhöhen und durch pi/2 dividieren:
43 if (TheLfloat(x)->len >= 2750) {
44 var cl_F_div_t q_r = cl_round_pi2(extend(x,TheLfloat(x)->len+1));
46 var cl_LF r = The(cl_LF)(q_r.remainder);
47 var cl_LF_cos_sin_t trig = cl_cossin_ratseries(r);
49 z = cl_float(trig.sin,x);
51 z = cl_float(trig.cos,x);
53 var cl_F_div_t q_r = cl_round_pi2(cl_F_extendsqrt(x));
55 var cl_LF r = The(cl_LF)(q_r.remainder);
56 var cl_LF y = sinx_naive(r); // y := sin(r)^2
59 z = cl_float(sqrt(y),x);
64 if (zerop(r) || (float_exponent(r) <= (-(sintC)float_digits(r))>>1))
65 z = cl_float(1,x); // cos(r) = 1.0
67 z = cl_float(sqrt(1 - y),x); // sqrt(1-y)
71 var cl_F_div_t q_r = cl_round_pi2(cl_F_extendsqrt(x));
73 var cl_F& r = q_r.remainder;
74 var cl_F y = sinxbyx_naive(r); // y := (sin(r)/r)^2
77 z = cl_float(r*sqrt(y),x);
80 if (zerop(r) || (float_exponent(r) <= (-(sintC)float_digits(r))>>1))
81 z = cl_float(1,x); // cos(r) = 1.0
83 z = cl_float(sqrt(1 - square(r)*y),x); // sqrt(1-r^2*y)
86 // evtl. Vorzeichenwechsel:
87 if (cl_I_to_UL(logand(q,2))==0)
93 // Timings of the two algorithms, on an i486 33 MHz, running Linux,
94 // applied to x = sqrt(2)-1 = 0.414...
104 // ==> ratseries faster for N >= 2750.