6 // Minds our p's and q's. The two last computed convergents.
11 typedef struct pqset_s
*pqset_ptr
;
12 typedef struct pqset_s pqset_t
[1];
14 void pqset_init(pqset_t pq
) {
21 void pqset_clear(pqset_t pq
) {
28 void pqset_neg(pqset_t pq
) {
29 mpz_neg(pq
->p
, pq
->p
);
30 mpz_neg(pq
->q
, pq
->q
);
31 mpz_neg(pq
->pold
, pq
->pold
);
32 mpz_neg(pq
->qold
, pq
->qold
);
35 void pqset_print(pqset_t pq
) {
36 gmp_printf("p's: %Zd %Zd\n", pq
->pold
, pq
->p
);
37 gmp_printf("q's: %Zd %Zd\n", pq
->qold
, pq
->q
);
40 // Compute the next convergent for regular continued fractions.
41 void pqset_regular_recur(pqset_t pq
, mpz_t denom
) {
42 mpz_addmul(pq
->pold
, denom
, pq
->p
);
43 mpz_swap(pq
->pold
, pq
->p
);
44 mpz_addmul(pq
->qold
, denom
, pq
->q
);
45 mpz_swap(pq
->qold
, pq
->q
);
48 // Compute the next convergent for nonregular continued fractions.
49 void pqset_nonregular_recur(pqset_t pq
, mpz_t num
, mpz_t denom
) {
50 mpz_mul(pq
->pold
, pq
->pold
, num
);
51 mpz_addmul(pq
->pold
, pq
->p
, denom
);
52 mpz_swap(pq
->pold
, pq
->p
);
53 mpz_mul(pq
->qold
, pq
->qold
, num
);
54 mpz_addmul(pq
->qold
, pq
->q
, denom
);
55 mpz_swap(pq
->qold
, pq
->q
);
58 // Get rid of nontrivial GCD for {p, q, pold, qold}.
59 // t0 and t1 are temporary variables.
60 void pqset_remove_gcd(pqset_ptr pq
, mpz_t t0
, mpz_t t1
) {
61 mpz_gcd(t0
, pq
->p
, pq
->q
);
62 mpz_gcd(t1
, pq
->pold
, pq
->qold
);
64 if (mpz_cmp_ui(t0
, 1)) {
65 mpz_divexact(pq
->pold
, pq
->pold
, t0
);
66 mpz_divexact(pq
->qold
, pq
->qold
, t0
);
67 mpz_divexact(pq
->p
, pq
->p
, t0
);
68 mpz_divexact(pq
->q
, pq
->q
, t0
);
72 // A Mobius transformation: four coefficients and the input.
73 // TODO: Use an array of size 4.
74 struct mobius_data_s
{
78 typedef struct mobius_data_s
*mobius_data_ptr
;
80 void pqset_set_mobius(pqset_t pq
, mobius_data_ptr md
) {
81 mpz_set(pq
->pold
, md
->b
); mpz_set(pq
->p
, md
->a
);
82 mpz_set(pq
->qold
, md
->d
); mpz_set(pq
->q
, md
->c
);
85 // Compute convergents of Mobius function applied to a regular
86 // continued fraction.
87 static void *mobius_convergent(cf_t cf
) {
88 mobius_data_ptr md
= cf_data(cf
);
89 cf_t input
= md
->input
;
92 pqset_set_mobius(pq
, md
);
98 pqset_regular_recur(pq
, denom
);
113 // Start a thread that, when signalled, computes the convergents of a Mobius
114 // transformation of a continued fraction.
115 cf_t
cf_new_mobius_convergent(cf_t x
, mpz_t a
, mpz_t b
, mpz_t c
, mpz_t d
) {
116 mobius_data_ptr md
= malloc(sizeof(*md
));
117 mpz_init(md
->a
); mpz_init(md
->b
); mpz_init(md
->c
); mpz_init(md
->d
);
118 mpz_set(md
->a
, a
); mpz_set(md
->b
, b
); mpz_set(md
->c
, c
); mpz_set(md
->d
, d
);
120 return cf_new(mobius_convergent
, md
);
123 // Start a thread that, when signalled, computes the convergents of a continued
125 cf_t
cf_new_cf_convergent(cf_t x
) {
127 mpz_init(one
); mpz_init(zero
);
128 mpz_set_ui(one
, 1); mpz_set_ui(zero
, 0);
129 cf_t res
= cf_new_mobius_convergent(x
, one
, zero
, zero
, one
);
130 mpz_clear(one
); mpz_clear(zero
);
134 // Compute nonregular convergents of a Mobius function applied
135 // to a nonregular continued fraction.
136 static void *nonregular_mobius_convergent(cf_t cf
) {
137 mobius_data_ptr md
= cf_data(cf
);
138 cf_t input
= md
->input
;
139 pqset_t pq
; pqset_init(pq
); pqset_set_mobius(pq
, md
);
140 mpz_t num
; mpz_init(num
);
141 mpz_t denom
; mpz_init(denom
);
142 mpz_t t0
, t1
; mpz_init(t0
); mpz_init(t1
);
144 pqset_nonregular_recur(pq
, num
, denom
);
145 pqset_remove_gcd(pq
, t0
, t1
);
151 cf_get(denom
, input
);
155 cf_get(denom
, input
);
162 mpz_clear(md
->a
); mpz_clear(md
->b
); mpz_clear(md
->c
); mpz_clear(md
->d
);
163 mpz_clear(t0
); mpz_clear(t1
);
168 cf_t
cf_new_nonregular_mobius_convergent(cf_t x
, mpz_t a
, mpz_t b
, mpz_t c
, mpz_t d
) {
169 mobius_data_ptr md
= malloc(sizeof(*md
));
170 mpz_init(md
->a
); mpz_init(md
->b
); mpz_init(md
->c
); mpz_init(md
->d
);
171 mpz_set(md
->a
, a
); mpz_set(md
->b
, b
); mpz_set(md
->c
, c
); mpz_set(md
->d
, d
);
173 return cf_new(nonregular_mobius_convergent
, md
);
176 // Input: Mobius transformation and nonregular continued fraction.
177 // Output: Regular continued fraction. Assumes input fraction is well-behaved.
178 static void *mobius_nonregular_throughput(cf_t cf
) {
179 mobius_data_ptr md
= cf_data(cf
);
180 cf_t input
= md
->input
;
181 pqset_t pq
; pqset_init(pq
); pqset_set_mobius(pq
, md
);
182 mpz_t num
; mpz_init(num
);
183 mpz_t denom
; mpz_init(denom
);
184 mpz_t t0
, t1
, t2
; mpz_init(t2
); mpz_init(t1
); mpz_init(t0
);
186 pqset_nonregular_recur(pq
, num
, denom
);
187 pqset_remove_gcd(pq
, t0
, t1
);
189 if (mpz_sgn(pq
->qold
)) {
190 mpz_fdiv_qr(t1
, t0
, pq
->pold
, pq
->qold
);
191 mpz_mul(t2
, t1
, pq
->q
);
193 if (mpz_cmp(t2
, pq
->p
) <= 0) {
194 mpz_add(t2
, t2
, pq
->q
);
195 if (mpz_cmp(t2
, pq
->p
) > 0) {
196 // Output continued fraction term.
198 // Subtract: remainder of p/q.
199 mpz_sub(t2
, t2
, pq
->p
);
200 mpz_sub(t2
, pq
->q
, t2
);
202 mpz_set(pq
->pold
, pq
->qold
);
203 mpz_set(pq
->qold
, t0
);
204 mpz_set(pq
->p
, pq
->q
);
213 cf_get(denom
, input
);
218 cf_get(denom
, input
);
225 mpz_clear(md
->a
); mpz_clear(md
->b
); mpz_clear(md
->c
); mpz_clear(md
->d
);
227 mpz_clear(t2
); mpz_clear(t1
); mpz_clear(t0
);
231 cf_t
cf_new_nonregular_to_cf(cf_t x
, mpz_t a
, mpz_t b
, mpz_t c
, mpz_t d
) {
232 mobius_data_ptr md
= malloc(sizeof(*md
));
233 mpz_init(md
->a
); mpz_init(md
->b
); mpz_init(md
->c
); mpz_init(md
->d
);
234 mpz_set(md
->a
, a
); mpz_set(md
->b
, b
); mpz_set(md
->c
, c
); mpz_set(md
->d
, d
);
236 return cf_new(mobius_nonregular_throughput
, md
);
239 // This seems to be slower than regularizing the continued fraction
240 // and then converting to decimal.
241 static void *nonregular_mobius_decimal(cf_t cf
) {
242 mobius_data_ptr md
= cf_data(cf
);
243 cf_t input
= md
->input
;
244 pqset_t pq
; pqset_init(pq
); pqset_set_mobius(pq
, md
);
245 mpz_t num
; mpz_init(num
);
246 mpz_t denom
; mpz_init(denom
);
247 mpz_t t0
, t1
, t2
; mpz_init(t2
); mpz_init(t1
); mpz_init(t0
);
249 pqset_nonregular_recur(pq
, num
, denom
);
250 pqset_remove_gcd(pq
, t0
, t1
);
252 // If the denominator is zero, we can't do anything yet.
253 if (mpz_sgn(pq
->qold
)) {
254 mpz_fdiv_qr(t1
, t0
, pq
->pold
, pq
->qold
);
255 mpz_mul(t2
, t1
, pq
->q
);
256 if (mpz_cmp(t2
, pq
->p
) <= 0) {
257 mpz_add(t2
, t2
, pq
->q
);
258 if (mpz_cmp(t2
, pq
->p
) > 0) {
259 // Output a decimal digit.
261 // Subtract: remainder of p/q.
262 mpz_sub(t2
, t2
, pq
->p
);
263 mpz_sub(t2
, pq
->q
, t2
);
264 // Multiply numerator by 10.
265 mpz_mul_ui(pq
->pold
, t0
, 10);
266 mpz_mul_ui(pq
->p
, t2
, 10);
274 cf_get(denom
, input
);
279 cf_get(denom
, input
);
285 mpz_clear(t0
); mpz_clear(t1
); mpz_clear(t2
);
286 mpz_clear(md
->a
); mpz_clear(md
->b
); mpz_clear(md
->c
); mpz_clear(md
->d
);
291 cf_t
cf_new_nonregular_mobius_to_decimal(cf_t x
, mpz_t a
, mpz_t b
, mpz_t c
, mpz_t d
) {
292 mobius_data_ptr md
= malloc(sizeof(*md
));
293 mpz_init(md
->a
); mpz_init(md
->b
); mpz_init(md
->c
); mpz_init(md
->d
);
294 mpz_set(md
->a
, a
); mpz_set(md
->b
, b
); mpz_set(md
->c
, c
); mpz_set(md
->d
, d
);
296 return cf_new(nonregular_mobius_decimal
, md
);
299 static void determine_sign(cf_t cf
, pqset_t pq
, mpz_t denom
, cf_t input
) {
300 cf_set_sign(cf
, cf_sign(input
));
301 while (mpz_sgn(pq
->pold
) != mpz_sgn(pq
->p
)
302 || mpz_sgn(pq
->qold
) != mpz_sgn(pq
->q
)) {
303 cf_get(denom
, input
);
304 pqset_regular_recur(pq
, denom
);
306 if (mpz_sgn(pq
->qold
) < 0) {
307 mpz_neg(pq
->qold
, pq
->qold
);
308 mpz_neg(pq
->q
, pq
->q
);
311 if (mpz_sgn(pq
->pold
) < 0) {
312 mpz_neg(pq
->pold
, pq
->pold
);
313 mpz_neg(pq
->p
, pq
->p
);
318 // Input: Mobius transformation and regular continued fraction.
319 // Output: Regular continued fraction.
320 static void *mobius_throughput(cf_t cf
) {
321 mobius_data_ptr md
= cf_data(cf
);
322 cf_t input
= md
->input
;
323 pqset_t pq
; pqset_init(pq
); pqset_set_mobius(pq
, md
);
324 mpz_t denom
; mpz_init(denom
);
325 mpz_t t0
, t1
, t2
; mpz_init(t2
); mpz_init(t1
); mpz_init(t0
);
327 determine_sign(cf
, pq
, denom
, input
);
330 pqset_regular_recur(pq
, denom
);
332 if (mpz_sgn(pq
->qold
)) {
333 mpz_fdiv_qr(t1
, t0
, pq
->pold
, pq
->qold
);
334 mpz_mul(t2
, t1
, pq
->q
);
336 if (mpz_cmp(t2
, pq
->p
) <= 0) {
337 mpz_add(t2
, t2
, pq
->q
);
338 if (mpz_cmp(t2
, pq
->p
) > 0) {
339 // Output continued fraction term.
341 // Subtract: remainder of p/q.
342 mpz_sub(t2
, t2
, pq
->p
);
343 mpz_sub(t2
, pq
->q
, t2
);
345 mpz_set(pq
->pold
, pq
->qold
);
346 mpz_set(pq
->qold
, t0
);
347 mpz_set(pq
->p
, pq
->q
);
357 cf_get(denom
, input
);
362 mpz_clear(t0
); mpz_clear(t1
); mpz_clear(t2
);
363 mpz_clear(md
->a
); mpz_clear(md
->b
); mpz_clear(md
->c
); mpz_clear(md
->d
);
368 cf_t
cf_new_mobius_to_cf(cf_t x
, mpz_t z
[4]) {
369 mobius_data_ptr md
= malloc(sizeof(*md
));
370 mpz_init(md
->a
); mpz_init(md
->b
); mpz_init(md
->c
); mpz_init(md
->d
);
371 mpz_set(md
->a
, z
[0]); mpz_set(md
->b
, z
[1]);
372 mpz_set(md
->c
, z
[2]); mpz_set(md
->d
, z
[3]);
374 return cf_new(mobius_throughput
, md
);
377 // Input: Mobius transformation and regular continued fraction.
378 // Output: Decimal representation. The integer part is given first,
379 // followed by one digit at a time.
380 static void *mobius_decimal(cf_t cf
) {
381 mobius_data_ptr md
= cf_data(cf
);
382 cf_t input
= md
->input
;
383 pqset_t pq
; pqset_init(pq
); pqset_set_mobius(pq
, md
);
384 mpz_t denom
; mpz_init(denom
);
385 mpz_t t0
, t1
, t2
; mpz_init(t2
); mpz_init(t1
); mpz_init(t0
);
387 determine_sign(cf
, pq
, denom
, input
);
389 pqset_regular_recur(pq
, denom
);
391 // If the denominator is zero, we can't do anything yet.
392 if (mpz_sgn(pq
->qold
)) {
393 // Each term except possibly the first is one of {0, ..., 9}.
394 /* Naive attempt to expoit this didn't seem faster:
395 * (and I'd have to handle the first term properly)
398 for (i = 0; i <= 9; i++) {
399 if (mpz_cmp(t0, pq->p) > 0) break;
400 mpz_add(t0, t0, pq->q);
402 mpz_set_ui(pq->pold, i);
403 mpz_sub(t0, t0, pq->p);
404 mpz_sub(t0, pq->q, t0);
405 mpz_mul(pq->qold, pq->pold, pq->qnew);
408 mpz_fdiv_qr(t1
, t0
, pq
->pold
, pq
->qold
);
409 mpz_mul(t2
, t1
, pq
->q
);
410 if (mpz_cmp(t2
, pq
->p
) <= 0) {
411 mpz_add(t2
, t2
, pq
->q
);
412 if (mpz_cmp(t2
, pq
->p
) > 0) {
413 // Output a decimal digit.
415 // Compute t2 = remainder of p/q.
416 mpz_sub(t2
, t2
, pq
->p
);
417 mpz_sub(t2
, pq
->q
, t2
);
418 // Multiply numerator by 10.
419 mpz_mul_ui(pq
->pold
, t0
, 10);
420 mpz_mul_ui(pq
->p
, t2
, 10);
429 cf_get(denom
, input
);
434 mpz_clear(t0
); mpz_clear(t1
); mpz_clear(t2
);
435 mpz_clear(md
->a
); mpz_clear(md
->b
); mpz_clear(md
->c
); mpz_clear(md
->d
);
440 cf_t
cf_new_mobius_to_decimal(cf_t x
, mpz_t a
, mpz_t b
, mpz_t c
, mpz_t d
) {
441 mobius_data_ptr md
= malloc(sizeof(*md
));
442 mpz_init(md
->a
); mpz_init(md
->b
); mpz_init(md
->c
); mpz_init(md
->d
);
443 mpz_set(md
->a
, a
); mpz_set(md
->b
, b
); mpz_set(md
->c
, c
); mpz_set(md
->d
, d
);
445 return cf_new(mobius_decimal
, md
);
448 cf_t
cf_new_cf_to_decimal(cf_t x
) {
450 mpz_init(one
); mpz_init(zero
);
451 mpz_set_ui(one
, 1); mpz_set_ui(zero
, 0);
452 cf_t res
= cf_new_mobius_to_decimal(x
, one
, zero
, zero
, one
);
453 mpz_clear(one
); mpz_clear(zero
);