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 struct pqrs_s {
  20   mpz_t p0, p1;
  21   mpz_t q0, q1;
  22   mpz_t r0, r1;
  23   mpz_t s0, s1;
  24 };
  25 typedef struct pqrs_s pqrs_t[1];
  26
  27 void pqrs_init(pqrs_t p) {
  28   mpz_init(p->p0); mpz_init(p->p1);
  29   mpz_init(p->q0); mpz_init(p->q1);
  30   mpz_init(p->r0); mpz_init(p->r1);
  31   mpz_init(p->s0); mpz_init(p->s1);
  32 }
  33
  34 void pqrs_clear(pqrs_t p) {
  35   mpz_clear(p->p0); mpz_clear(p->p1);
  36   mpz_clear(p->q0); mpz_clear(p->q1);
  37   mpz_clear(p->r0); mpz_clear(p->r1);
  38   mpz_clear(p->s0); mpz_clear(p->s1);
  39 }
  40
  41 void pqrs_set_coeff(pqrs_t p, mpz_t a[8]) {
  42   mpz_set(p->p0, a[0]); mpz_set(p->q0, a[1]);
  43   mpz_set(p->r0, a[2]); mpz_set(p->s0, a[3]);
  44   mpz_set(p->p1, a[4]); mpz_set(p->q1, a[5]);
  45   mpz_set(p->r1, a[6]); mpz_set(p->s1, a[7]);
  46 }
  47
  48 void pqrs_print(pqrs_t p) {
  49   gmp_printf("%Zd/%Zd %Zd/%Zd\n", p->s0, p->s1, p->q0, p->q1);
  50   gmp_printf("%Zd/%Zd %Zd/%Zd\n", p->r0, p->r1, p->p0, p->p1);
  51 }
  52
  53 static void *bihom(cf_t cf) {
  54   bihom_data_ptr bd = cf_data(cf);
  55   pqrs_t p;
  56   pqrs_init(p);
  57   pqrs_t qr;  // For quotient and remainders.
  58   pqrs_init(qr);
  59   pqrs_set_coeff(p, bd->a);
  60   cf_t x = bd->x;
  61   cf_t y = bd->y;
  62   mpz_t z, t0, t1;
  63   mpz_init(z);
  64   mpz_init(t0); mpz_init(t1);
  65   void move_down() {
  66     cf_get(z, y);
  67     mpz_mul(t0, z, p->r0);  mpz_mul(t1, z, p->r1);
  68     mpz_add(t0, t0, p->s0); mpz_add(t1, t1, p->s1);
  69     mpz_set(p->s0, p->r0);  mpz_set(p->s1, p->r1);
  70     mpz_set(p->r0, t0);     mpz_set(p->r1, t1);
  71
  72     mpz_mul(t0, z, p->p0);  mpz_mul(t1, z, p->p1);
  73     mpz_add(t0, t0, p->q0); mpz_add(t1, t1, p->q1);
  74     mpz_set(p->q0, p->p0);  mpz_set(p->q1, p->p1);
  75     mpz_set(p->p0, t0);     mpz_set(p->p1, t1);
  76   }
  77   void move_right() {
  78     cf_get(z, x);
  79     mpz_mul(t0, z, p->q0);  mpz_mul(t1, z, p->q1);
  80     mpz_add(t0, t0, p->s0); mpz_add(t1, t1, p->s1);
  81     mpz_set(p->s0, p->q0);  mpz_set(p->s1, p->q1);
  82     mpz_set(p->q0, t0);     mpz_set(p->q1, t1);
  83
  84     mpz_mul(t0, z, p->p0);  mpz_mul(t1, z, p->p1);
  85     mpz_add(t0, t0, p->r0); mpz_add(t1, t1, p->r1);
  86     mpz_set(p->r0, p->p0);  mpz_set(p->r1, p->p1);
  87     mpz_set(p->p0, t0);     mpz_set(p->p1, t1);
  88   }
  89   int recur() {
  90     if (!mpz_sgn(p->s1)) {
  91       move_right();
  92       move_down();
  93       return 0;
  94     }
  95     if (!mpz_sgn(p->p1)) {
  96       if (!mpz_sgn(p->q1)) {
  97         move_down();
  98       } else {
  99         move_right();
 100       }
 101       return 0;
 102     }
 103     if (!mpz_sgn(p->q1)) {
 104       move_down();
 105       return 0;
 106     }
 107     if (!mpz_sgn(p->r1)) {
 108       move_right();
 109       return 0;
 110     }
 111     mpz_fdiv_qr(qr->p0, qr->p1, p->p0, p->p1);
 112     // TODO: Optimize by removing other divisions.
 113     mpz_fdiv_qr(qr->q0, qr->q1, p->q0, p->q1);
 114     if (mpz_cmp(qr->p0, qr->q0)) {
 115       move_down();
 116       return 0;
 117     }
 118     mpz_fdiv_qr(qr->r0, qr->r1, p->r0, p->r1);
 119     if (mpz_cmp(qr->p0, qr->r0)) {
 120       move_right();
 121       return 0;
 122     }
 123     mpz_fdiv_qr(qr->s0, qr->s1, p->s0, p->s1);
 124     if (mpz_cmp(qr->q0, qr->q0)) {
 125       move_down();
 126       return 0;
 127     }
 128     cf_put(cf, qr->p0);
 129     mpz_set(p->p0, p->p1); mpz_set(p->p1, qr->p1);
 130     mpz_set(p->q0, p->q1); mpz_set(p->q1, qr->q1);
 131     mpz_set(p->r0, p->r1); mpz_set(p->r1, qr->r1);
 132     mpz_set(p->s0, p->s1); mpz_set(p->s1, qr->s1);
 133     return 1;
 134   }
 135   while(cf_wait(cf)) {
 136     while(!recur());
 137   }
 138   pqrs_clear(p);
 139   pqrs_clear(qr);
 140   mpz_clear(z);
 141   mpz_clear(t0); mpz_clear(t1);
 142   return NULL;
 143 }
 144
 145 cf_t cf_new_bihom(cf_t x, cf_t y, mpz_t a[8]) {
 146   bihom_data_ptr p = malloc(sizeof(*p));
 147   p->x = x;
 148   p->y = y;
 149   for (int i = 0; i < 8; i++) {
 150     mpz_init(p->a[i]);
 151     mpz_set(p->a[i], a[i]);
 152   }
 153   return cf_new(bihom, p);
 154 }