6 // sin 1 from Taylor series
7 static void *sin1_expansion(cf_t cf
) {
8 mpz_t p0
, q0
; // p0/q0 = last convergent
9 mpz_t p1
, q1
; // p1/q1 = convergent
10 mpq_t r0
; // lower bound for sin 1
11 mpq_t r1
; // upper bound for sin 1
12 mpz_t t0
, t1
, t2
; // Temporary variables.
15 mpq_init(r0
); mpq_init(r1
);
16 mpz_init(p0
); mpz_init(p1
); mpz_init(q0
); mpz_init(q1
);
17 mpz_init(t0
); mpz_init(t1
); mpz_init(t2
);
20 // Zero causes all sorts of problems.
21 // Output it and forget about it.
24 // Initialize convergents.
28 // sin 1 = 1 - 1/6 + ...
32 int k
= 1; // Number of Taylor series terms we've computed,
33 // numbered from 0, unlike continued fraction terms.
34 int i
= 2; // Number of convergent we're computing.
38 mpz_mul_ui(mpq_numref(r1
), mpq_numref(r0
), (2 * k
) * (2 * k
+ 1));
39 mpz_mul_ui(mpq_denref(r1
), mpq_denref(r0
), (2 * k
) * (2 * k
+ 1));
40 mpz_add_ui(mpq_numref(r1
), mpq_numref(r1
), 1);
42 mpz_mul_ui(mpq_numref(r0
), mpq_numref(r1
), (2 * k
) * (2 * k
+ 1));
43 mpz_mul_ui(mpq_denref(r0
), mpq_denref(r1
), (2 * k
) * (2 * k
+ 1));
44 mpz_sub_ui(mpq_numref(r0
), mpq_numref(r0
), 1);
48 // Find largest x such that (p1 x + p0)/(q1 x + q0)
49 // <= r0 for odd i, and >= r1 for even i
50 // => x < (q0 num(r) - p0 den(r))/(p1 den(r) - q1 num(r))
51 // for the appropriate r.
52 mpq_ptr r
= i
& 1 ? r0
: r1
;
53 mpz_mul(t0
, q0
, mpq_numref(r
));
54 mpz_mul(t1
, p0
, mpq_denref(r
));
56 mpz_mul(t0
, p1
, mpq_denref(r
));
57 mpz_mul(t1
, q1
, mpq_numref(r
));
59 mpz_fdiv_q(t0
, t2
, t1
);
61 // Then x is the next convergent provided substituting x + 1 overshoots
62 // the other bound. (If not, we need better bounds to decide.)
66 mpz_add_ui(t1
, t0
, 1);
67 mpz_set(mpq_numref(tq0
), p0
);
68 mpz_addmul(mpq_numref(tq0
), p1
, t1
);
69 mpz_set(mpq_denref(tq0
), q0
);
70 mpz_addmul(mpq_denref(tq0
), q1
, t1
);
72 if (mpq_cmp(tq0
, r1
) > 0) break;
73 // Check (p1(x+1) + p0)/(q1(x+1) + q0) > r1
75 if (mpq_cmp(tq0
, r0
) < 0) break;
76 // Check (p1(x+1) + p0)/(q1(x+1) + q0) < r0
83 mpz_addmul(p0
, p1
, t0
);
84 mpz_addmul(q0
, q1
, t0
);
89 mpq_clear(r0
); mpq_clear(r1
);
90 mpz_clear(p0
); mpz_clear(p1
); mpz_clear(q0
); mpz_clear(q1
);
91 mpz_clear(t0
); mpz_clear(t1
); mpz_clear(t2
);
97 return cf_new_const(sin1_expansion
);
100 // cos 1 from Taylor series
101 static void *cos1_expansion(cf_t cf
) {
102 mpz_t p0
, q0
; // p0/q0 = last convergent
103 mpz_t p1
, q1
; // p1/q1 = convergent
104 mpq_t r0
; // lower bound for sin 1
105 mpq_t r1
; // upper bound for sin 1
106 mpz_t t0
, t1
, t2
; // Temporary variables.
109 mpq_init(r0
); mpq_init(r1
);
110 mpz_init(p0
); mpz_init(p1
); mpz_init(q0
); mpz_init(q1
);
111 mpz_init(t0
); mpz_init(t1
); mpz_init(t2
);
114 // Zero causes all sorts of problems.
115 // Output it and forget about it.
118 // Initialize convergents.
122 // cos 1 = 1 - 1/2 + ...
123 mpq_set_ui(r0
, 1, 2);
124 mpq_set_ui(r1
, 1, 1);
126 int k
= 1; // Number of Taylor series terms we've computed,
127 // numbered from 0, unlike continued fraction terms.
128 int i
= 2; // Number of convergent we're computing.
132 mpz_mul_ui(mpq_numref(r1
), mpq_numref(r0
), (2 * k
) * (2 * k
- 1));
133 mpz_mul_ui(mpq_denref(r1
), mpq_denref(r0
), (2 * k
) * (2 * k
- 1));
134 mpz_add_ui(mpq_numref(r1
), mpq_numref(r1
), 1);
136 mpz_mul_ui(mpq_numref(r0
), mpq_numref(r1
), (2 * k
) * (2 * k
- 1));
137 mpz_mul_ui(mpq_denref(r0
), mpq_denref(r1
), (2 * k
) * (2 * k
- 1));
138 mpz_sub_ui(mpq_numref(r0
), mpq_numref(r0
), 1);
140 while (cf_wait(cf
)) {
142 // Find largest x such that (p1 x + p0)/(q1 x + q0)
143 // <= r0 for odd i, and >= r1 for even i
144 // => x < (q0 num(r) - p0 den(r))/(p1 den(r) - q1 num(r))
145 // for the appropriate r.
146 mpq_ptr r
= i
& 1 ? r0
: r1
;
147 mpz_mul(t0
, q0
, mpq_numref(r
));
148 mpz_mul(t1
, p0
, mpq_denref(r
));
150 mpz_mul(t0
, p1
, mpq_denref(r
));
151 mpz_mul(t1
, q1
, mpq_numref(r
));
153 mpz_fdiv_q(t0
, t2
, t1
);
155 // Then x is the next convergent provided substituting x + 1 overshoots
156 // the other bound. (If not, we need better bounds to decide.)
160 mpz_add_ui(t1
, t0
, 1);
161 mpz_set(mpq_numref(tq0
), p0
);
162 mpz_addmul(mpq_numref(tq0
), p1
, t1
);
163 mpz_set(mpq_denref(tq0
), q0
);
164 mpz_addmul(mpq_denref(tq0
), q1
, t1
);
166 if (mpq_cmp(tq0
, r1
) > 0) break;
167 // Check (p1(x+1) + p0)/(q1(x+1) + q0) > r1
169 if (mpq_cmp(tq0
, r0
) < 0) break;
170 // Check (p1(x+1) + p0)/(q1(x+1) + q0) < r0
177 mpz_addmul(p0
, p1
, t0
);
178 mpz_addmul(q0
, q1
, t0
);
183 mpq_clear(r0
); mpq_clear(r1
);
184 mpz_clear(p0
); mpz_clear(p1
); mpz_clear(q0
); mpz_clear(q1
);
185 mpz_clear(t0
); mpz_clear(t1
); mpz_clear(t2
);
191 return cf_new_const(cos1_expansion
);