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 }