1 // Compute bihomographic functions of two continued fractions.
11 typedef struct bihom_data_s bihom_data_t
[1];
12 typedef struct bihom_data_s
*bihom_data_ptr
;
14 // In the 3D table, the four convergents are associated with the letter:
17 // e.g. the s convergent is s0/s1.
18 // I work left to right, and top to bottom, unlike Gosper.
20 // This notation is different to that in cf_mobius.c. Here, the same fraction
21 // has the same letter. There, the numerators were p's, and the denominators
29 typedef struct pqrs_s pqrs_t
[1];
31 void pqrs_init(pqrs_t p
) {
32 mpz_init(p
->p0
); mpz_init(p
->p1
);
33 mpz_init(p
->q0
); mpz_init(p
->q1
);
34 mpz_init(p
->r0
); mpz_init(p
->r1
);
35 mpz_init(p
->s0
); mpz_init(p
->s1
);
38 void pqrs_clear(pqrs_t p
) {
39 mpz_clear(p
->p0
); mpz_clear(p
->p1
);
40 mpz_clear(p
->q0
); mpz_clear(p
->q1
);
41 mpz_clear(p
->r0
); mpz_clear(p
->r1
);
42 mpz_clear(p
->s0
); mpz_clear(p
->s1
);
45 void pqrs_set_coeff(pqrs_t p
, mpz_t a
[8]) {
46 mpz_set(p
->p0
, a
[0]); mpz_set(p
->q0
, a
[1]);
47 mpz_set(p
->r0
, a
[2]); mpz_set(p
->s0
, a
[3]);
48 mpz_set(p
->p1
, a
[4]); mpz_set(p
->q1
, a
[5]);
49 mpz_set(p
->r1
, a
[6]); mpz_set(p
->s1
, a
[7]);
52 void pqrs_print(pqrs_t p
) {
53 gmp_printf("%Zd/%Zd %Zd/%Zd\n", p
->s0
, p
->s1
, p
->q0
, p
->q1
);
54 gmp_printf("%Zd/%Zd %Zd/%Zd\n", p
->r0
, p
->r1
, p
->p0
, p
->p1
);
57 static void *bihom(cf_t cf
) {
58 bihom_data_ptr bd
= cf_data(cf
);
61 pqrs_t qr
; // For quotient and remainders.
63 pqrs_set_coeff(p
, bd
->a
);
68 mpz_init(t0
); mpz_init(t1
);
71 mpz_mul(t0
, z
, p
->r0
); mpz_mul(t1
, z
, p
->r1
);
72 mpz_add(t0
, t0
, p
->s0
); mpz_add(t1
, t1
, p
->s1
);
73 mpz_set(p
->s0
, p
->r0
); mpz_set(p
->s1
, p
->r1
);
74 mpz_set(p
->r0
, t0
); mpz_set(p
->r1
, t1
);
76 mpz_mul(t0
, z
, p
->p0
); mpz_mul(t1
, z
, p
->p1
);
77 mpz_add(t0
, t0
, p
->q0
); mpz_add(t1
, t1
, p
->q1
);
78 mpz_set(p
->q0
, p
->p0
); mpz_set(p
->q1
, p
->p1
);
79 mpz_set(p
->p0
, t0
); mpz_set(p
->p1
, t1
);
83 mpz_mul(t0
, z
, p
->q0
); mpz_mul(t1
, z
, p
->q1
);
84 mpz_add(t0
, t0
, p
->s0
); mpz_add(t1
, t1
, p
->s1
);
85 mpz_set(p
->s0
, p
->q0
); mpz_set(p
->s1
, p
->q1
);
86 mpz_set(p
->q0
, t0
); mpz_set(p
->q1
, t1
);
88 mpz_mul(t0
, z
, p
->p0
); mpz_mul(t1
, z
, p
->p1
);
89 mpz_add(t0
, t0
, p
->r0
); mpz_add(t1
, t1
, p
->r1
);
90 mpz_set(p
->r0
, p
->p0
); mpz_set(p
->r1
, p
->p1
);
91 mpz_set(p
->p0
, t0
); mpz_set(p
->p1
, t1
);
94 if (!mpz_sgn(p
->s1
)) {
99 if (!mpz_sgn(p
->p1
)) {
100 if (!mpz_sgn(p
->q1
)) {
107 if (!mpz_sgn(p
->q1
)) {
111 if (!mpz_sgn(p
->r1
)) {
115 mpz_fdiv_qr(qr
->p0
, qr
->p1
, p
->p0
, p
->p1
);
116 // TODO: Optimize by removing other divisions.
117 mpz_fdiv_qr(qr
->q0
, qr
->q1
, p
->q0
, p
->q1
);
118 if (mpz_cmp(qr
->p0
, qr
->q0
)) {
122 mpz_fdiv_qr(qr
->r0
, qr
->r1
, p
->r0
, p
->r1
);
123 if (mpz_cmp(qr
->p0
, qr
->r0
)) {
127 mpz_fdiv_qr(qr
->s0
, qr
->s1
, p
->s0
, p
->s1
);
128 if (mpz_cmp(qr
->s0
, qr
->r0
)) {
129 move_down(); // Either way should work.
133 mpz_set(p
->p0
, p
->p1
); mpz_set(p
->p1
, qr
->p1
);
134 mpz_set(p
->q0
, p
->q1
); mpz_set(p
->q1
, qr
->q1
);
135 mpz_set(p
->r0
, p
->r1
); mpz_set(p
->r1
, qr
->r1
);
136 mpz_set(p
->s0
, p
->s1
); mpz_set(p
->s1
, qr
->s1
);
145 mpz_clear(t0
); mpz_clear(t1
);
149 cf_t
cf_new_bihom(cf_t x
, cf_t y
, mpz_t a
[8]) {
150 bihom_data_ptr p
= malloc(sizeof(*p
));
153 for (int i
= 0; i
< 8; i
++) {
155 mpz_set(p
->a
[i
], a
[i
]);
157 return cf_new(bihom
, p
);