bihom.c

   1 // Compute bihomographic functions of two continued fractions.
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 #include <gmp.h>
   5 #include "cf.h"
   6
   7 struct bihom_data_s {
   8   cf_t x, y;
   9   mpz_t a[8];
  10 };
  11 typedef struct bihom_data_s bihom_data_t[1];
  12 typedef struct bihom_data_s *bihom_data_ptr;
  13
  14 // In the 3D table, the four convergents are associated with the letter:
  15 //   s  q
  16 //   r  p
  17 // e.g. the s convergent is s0/s1.
  18 // I work left to right, and top to bottom, unlike Gosper.
  19 //
  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
  22 // were q's.
  23 struct pqrs_s {
  24   mpz_t p0, p1;
  25   mpz_t q0, q1;
  26   mpz_t r0, r1;
  27   mpz_t s0, s1;
  28 };
  29 typedef struct pqrs_s pqrs_t[1];
  30
  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);
  36 }
  37
  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);
  43 }
  44
  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]);
  50 }
  51
  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);
  55 }
  56
  57 static void *bihom(cf_t cf) {
  58   bihom_data_ptr bd = cf_data(cf);
  59   pqrs_t p;
  60   pqrs_init(p);
  61   pqrs_t qr;  // For quotient and remainders.
  62   pqrs_init(qr);
  63   pqrs_set_coeff(p, bd->a);
  64   cf_t x = bd->x;
  65   cf_t y = bd->y;
  66   mpz_t z, t0, t1;
  67   mpz_init(z);
  68   mpz_init(t0); mpz_init(t1);
  69   void move_down() {
  70     cf_get(z, y);
  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);
  75
  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);
  80   }
  81   void move_right() {
  82     cf_get(z, x);
  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);
  87
  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);
  92   }
  93   // Determine sign.
  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)) {
 100     move_right();
 101     move_down();
 102   }
 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);
 108     cf_flip_sign(cf);
 109   }
 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);
 115     cf_flip_sign(cf);
 116   }
 117
 118   int recur() {
 119     if (!mpz_sgn(p->s1)) {
 120       move_right();
 121       move_down();
 122       return 0;
 123     }
 124     if (!mpz_sgn(p->p1)) {
 125       if (!mpz_sgn(p->q1)) {
 126         move_down();
 127       } else {
 128         move_right();
 129       }
 130       return 0;
 131     }
 132     if (!mpz_sgn(p->q1)) {
 133       move_down();
 134       return 0;
 135     }
 136     if (!mpz_sgn(p->r1)) {
 137       move_right();
 138       return 0;
 139     }
 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)) {
 144       move_down();
 145       return 0;
 146     }
 147     mpz_fdiv_qr(qr->r0, qr->r1, p->r0, p->r1);
 148     if (mpz_cmp(qr->p0, qr->r0)) {
 149       move_right();
 150       return 0;
 151     }
 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.
 155       return 0;
 156     }
 157     cf_put(cf, qr->p0);
 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);
 162     return 1;
 163   }
 164   while(cf_wait(cf)) {
 165     while(!recur());
 166   }
 167   pqrs_clear(p);
 168   pqrs_clear(qr);
 169   mpz_clear(z);
 170   mpz_clear(t0); mpz_clear(t1);
 171   for (int i = 0; i < 8; i++) {
 172     mpz_clear(bd->a[i]);
 173   }
 174   free(bd);
 175   return NULL;
 176 }
 177
 178 cf_t cf_new_bihom(cf_t x, cf_t y, mpz_t a[8]) {
 179   bihom_data_ptr p = malloc(sizeof(*p));
 180   p->x = x;
 181   p->y = y;
 182   for (int i = 0; i < 8; i++) {
 183     mpz_init(p->a[i]);
 184     mpz_set(p->a[i], a[i]);
 185   }
 186   return cf_new(bihom, p);
 187 }
 188
 189 void mpz8_init(mpz_t z[8]) {
 190   int i;
 191   for(i = 0; i < 8; i++) {
 192     mpz_init(z[i]);
 193   }
 194 }
 195
 196 void mpz8_clear(mpz_t z[8]) {
 197   int i;
 198   for(i = 0; i < 8; i++) {
 199     mpz_clear(z[i]);
 200   }
 201 }
 202
 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) {
 206   mpz_set_si(z[0], a);
 207   mpz_set_si(z[1], b);
 208   mpz_set_si(z[2], c);
 209   mpz_set_si(z[3], d);
 210   mpz_set_si(z[4], e);
 211   mpz_set_si(z[5], f);
 212   mpz_set_si(z[6], g);
 213   mpz_set_si(z[7], h);
 214 }
 215
 216 void mpz8_set_add(mpz_t z[8]) {
 217   mpz8_set_int(z,
 218     0, 1, 1, 0,
 219     0, 0, 0, 1);
 220 }
 221
 222 void mpz8_set_sub(mpz_t z[8]) {
 223   mpz8_set_int(z,
 224     0, 1, -1, 0,
 225     0, 0, 0, 1);
 226 }
 227
 228 void mpz8_set_mul(mpz_t z[8]) {
 229   mpz8_set_int(z,
 230     1, 0, 0, 0,
 231     0, 0, 0, 1);
 232 }
 233
 234 void mpz8_set_div(mpz_t z[8]) {
 235   mpz8_set_int(z,
 236     0, 1, 0, 0,
 237     0, 0, 1, 0);
 238 }
 239
 240 cf_t cf_new_add(cf_t x, cf_t y) {
 241   bihom_data_ptr p = malloc(sizeof(*p));
 242   p->x = x;
 243   p->y = y;
 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);
 248
 249   return cf_new(bihom, p);
 250 }
 251
 252 cf_t cf_new_sub(cf_t x, cf_t y) {
 253   bihom_data_ptr p = malloc(sizeof(*p));
 254   p->x = x;
 255   p->y = y;
 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);
 260
 261   return cf_new(bihom, p);
 262 }
 263
 264 cf_t cf_new_mul(cf_t x, cf_t y) {
 265   bihom_data_ptr p = malloc(sizeof(*p));
 266   p->x = x;
 267   p->y = y;
 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);
 271
 272   return cf_new(bihom, p);
 273 }
 274
 275 cf_t cf_new_div(cf_t x, cf_t y) {
 276   bihom_data_ptr p = malloc(sizeof(*p));
 277   p->x = x;
 278   p->y = y;
 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);
 282
 283   return cf_new(bihom, p);
 284 }