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 while (mpz_sgn(p
->p1
) != mpz_sgn(p
->q1
)
95 || mpz_sgn(p
->q1
) != mpz_sgn(p
->r1
)
96 || mpz_sgn(p
->r1
) != mpz_sgn(p
->s1
)
97 || mpz_sgn(p
->p0
) != mpz_sgn(p
->q0
)
98 || mpz_sgn(p
->q0
) != mpz_sgn(p
->r0
)
99 || mpz_sgn(p
->r0
) != mpz_sgn(p
->s0
)) {
103 if (mpz_sgn(p
->p0
) < 0) {
104 mpz_neg(p
->p0
, p
->p0
);
105 mpz_neg(p
->q0
, p
->q0
);
106 mpz_neg(p
->r0
, p
->r0
);
107 mpz_neg(p
->s0
, p
->s0
);
110 if (mpz_sgn(p
->p1
) < 0) {
111 mpz_neg(p
->p1
, p
->p1
);
112 mpz_neg(p
->q1
, p
->q1
);
113 mpz_neg(p
->r1
, p
->r1
);
114 mpz_neg(p
->s1
, p
->s1
);
119 if (!mpz_sgn(p
->s1
)) {
124 if (!mpz_sgn(p
->p1
)) {
125 if (!mpz_sgn(p
->q1
)) {
132 if (!mpz_sgn(p
->q1
)) {
136 if (!mpz_sgn(p
->r1
)) {
140 mpz_fdiv_qr(qr
->p0
, qr
->p1
, p
->p0
, p
->p1
);
141 // TODO: Optimize by removing other divisions.
142 mpz_fdiv_qr(qr
->q0
, qr
->q1
, p
->q0
, p
->q1
);
143 if (mpz_cmp(qr
->p0
, qr
->q0
)) {
147 mpz_fdiv_qr(qr
->r0
, qr
->r1
, p
->r0
, p
->r1
);
148 if (mpz_cmp(qr
->p0
, qr
->r0
)) {
152 mpz_fdiv_qr(qr
->s0
, qr
->s1
, p
->s0
, p
->s1
);
153 if (mpz_cmp(qr
->s0
, qr
->r0
)) {
154 move_down(); // Either way should work.
158 mpz_set(p
->p0
, p
->p1
); mpz_set(p
->p1
, qr
->p1
);
159 mpz_set(p
->q0
, p
->q1
); mpz_set(p
->q1
, qr
->q1
);
160 mpz_set(p
->r0
, p
->r1
); mpz_set(p
->r1
, qr
->r1
);
161 mpz_set(p
->s0
, p
->s1
); mpz_set(p
->s1
, qr
->s1
);
170 mpz_clear(t0
); mpz_clear(t1
);
171 for (int i
= 0; i
< 8; i
++) {
178 cf_t
cf_new_bihom(cf_t x
, cf_t y
, mpz_t a
[8]) {
179 bihom_data_ptr p
= malloc(sizeof(*p
));
182 for (int i
= 0; i
< 8; i
++) {
184 mpz_set(p
->a
[i
], a
[i
]);
186 return cf_new(bihom
, p
);
189 void mpz8_init(mpz_t z
[8]) {
191 for(i
= 0; i
< 8; i
++) {
196 void mpz8_clear(mpz_t z
[8]) {
198 for(i
= 0; i
< 8; i
++) {
203 void mpz8_set_int(mpz_t z
[8],
204 int a
, int b
, int c
, int d
,
205 int e
, int f
, int g
, int h
) {
216 void mpz8_set_add(mpz_t z
[8]) {
222 void mpz8_set_sub(mpz_t z
[8]) {
228 void mpz8_set_mul(mpz_t z
[8]) {
234 void mpz8_set_div(mpz_t z
[8]) {
240 cf_t
cf_new_add(cf_t x
, cf_t y
) {
241 bihom_data_ptr p
= malloc(sizeof(*p
));
244 for (int i
= 0; i
< 8; i
++) mpz_init(p
->a
[i
]);
245 mpz_set_si(p
->a
[1], 1);
246 mpz_set_si(p
->a
[2], 1);
247 mpz_set_si(p
->a
[7], 1);
249 return cf_new(bihom
, p
);
252 cf_t
cf_new_sub(cf_t x
, cf_t y
) {
253 bihom_data_ptr p
= malloc(sizeof(*p
));
256 for (int i
= 0; i
< 8; i
++) mpz_init(p
->a
[i
]);
257 mpz_set_si(p
->a
[1], 1);
258 mpz_set_si(p
->a
[2], -1);
259 mpz_set_si(p
->a
[7], 1);
261 return cf_new(bihom
, p
);
264 cf_t
cf_new_mul(cf_t x
, cf_t y
) {
265 bihom_data_ptr p
= malloc(sizeof(*p
));
268 for (int i
= 0; i
< 8; i
++) mpz_init(p
->a
[i
]);
269 mpz_set_si(p
->a
[0], 1);
270 mpz_set_si(p
->a
[7], 1);
272 return cf_new(bihom
, p
);
275 cf_t
cf_new_div(cf_t x
, cf_t y
) {
276 bihom_data_ptr p
= malloc(sizeof(*p
));
279 for (int i
= 0; i
< 8; i
++) mpz_init(p
->a
[i
]);
280 mpz_set_si(p
->a
[1], 1);
281 mpz_set_si(p
->a
[6], 1);
283 return cf_new(bihom
, p
);