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
  94   int recur() {
  95     if (!mpz_sgn(p->s1)) {
  96       move_right();
  97       move_down();
  98       return 0;
  99     }
 100     if (!mpz_sgn(p->p1)) {
 101       if (!mpz_sgn(p->q1)) {
 102         move_down();
 103       } else {
 104         move_right();
 105       }
 106       return 0;
 107     }
 108     if (!mpz_sgn(p->q1)) {
 109       move_down();
 110       return 0;
 111     }
 112     if (!mpz_sgn(p->r1)) {
 113       move_right();
 114       return 0;
 115     }
 116     mpz_fdiv_qr(qr->p0, qr->p1, p->p0, p->p1);
 117     // TODO: Optimize by removing other divisions.
 118     mpz_fdiv_qr(qr->q0, qr->q1, p->q0, p->q1);
 119     if (mpz_cmp(qr->p0, qr->q0)) {
 120       move_down();
 121       return 0;
 122     }
 123     mpz_fdiv_qr(qr->r0, qr->r1, p->r0, p->r1);
 124     if (mpz_cmp(qr->p0, qr->r0)) {
 125       move_right();
 126       return 0;
 127     }
 128     mpz_fdiv_qr(qr->s0, qr->s1, p->s0, p->s1);
 129     if (mpz_cmp(qr->s0, qr->r0)) {
 130       move_down();  // Either way should work.
 131       return 0;
 132     }
 133     cf_put(cf, qr->p0);
 134     mpz_set(p->p0, p->p1); mpz_set(p->p1, qr->p1);
 135     mpz_set(p->q0, p->q1); mpz_set(p->q1, qr->q1);
 136     mpz_set(p->r0, p->r1); mpz_set(p->r1, qr->r1);
 137     mpz_set(p->s0, p->s1); mpz_set(p->s1, qr->s1);
 138     return 1;
 139   }
 140   while(cf_wait(cf)) {
 141     while(!recur());
 142   }
 143   pqrs_clear(p);
 144   pqrs_clear(qr);
 145   mpz_clear(z);
 146   mpz_clear(t0); mpz_clear(t1);
 147   for (int i = 0; i < 8; i++) {
 148     mpz_clear(bd->a[i]);
 149   }
 150   free(bd);
 151   return NULL;
 152 }
 153
 154 cf_t cf_new_bihom(cf_t x, cf_t y, mpz_t a[8]) {
 155   bihom_data_ptr p = malloc(sizeof(*p));
 156   p->x = x;
 157   p->y = y;
 158   for (int i = 0; i < 8; i++) {
 159     mpz_init(p->a[i]);
 160     mpz_set(p->a[i], a[i]);
 161   }
 162   return cf_new(bihom, p);
 163 }