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