cf_mobius.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <gmp.h>
   4 #include "cf.h"
   5
   6 // Minds our p's and q's. The two last computed convergents,
   7 // and room to compute the next one.
   8 struct pqset_s {
   9   mpz_t pold, p, pnew;
  10   mpz_t qold, q, qnew;
  11 };
  12 typedef struct pqset_s *pqset_ptr;
  13 typedef struct pqset_s pqset_t[1];
  14
  15 void pqset_init(pqset_t pq) {
  16   mpz_init(pq->pold);
  17   mpz_init(pq->qold);
  18   mpz_init(pq->p);
  19   mpz_init(pq->q);
  20   mpz_init(pq->pnew);
  21   mpz_init(pq->qnew);
  22 }
  23
  24 void pqset_clear(pqset_t pq) {
  25   mpz_clear(pq->pold);
  26   mpz_clear(pq->qold);
  27   mpz_clear(pq->p);
  28   mpz_clear(pq->q);
  29   mpz_clear(pq->pnew);
  30   mpz_clear(pq->qnew);
  31 }
  32
  33 void pqset_print(pqset_t pq) {
  34   gmp_printf("p's: %Zd %Zd %Zd\n", pq->pold, pq->p, pq->pnew);
  35   gmp_printf("q's: %Zd %Zd %Zd\n", pq->qold, pq->q, pq->qnew);
  36 }
  37
  38 // Compute the next convergent for regular continued fractions.
  39 void pqset_regular_recur(pqset_t pq, mpz_t denom) {
  40   mpz_mul(pq->pnew, pq->p, denom);
  41   mpz_add(pq->pnew, pq->pnew, pq->pold);
  42   mpz_mul(pq->qnew, pq->q, denom);
  43   mpz_add(pq->qnew, pq->qnew, pq->qold);
  44 }
  45
  46 // Compute the next convergent for nonregular continued fractions.
  47 void pqset_nonregular_recur(pqset_t pq, mpz_t num, mpz_t denom) {
  48   mpz_mul(pq->pnew, pq->p, denom);
  49   mpz_addmul(pq->pnew, pq->pold, num);
  50   mpz_mul(pq->qnew, pq->q, denom);
  51   mpz_addmul(pq->qnew, pq->qold, num);
  52 }
  53
  54 // Use pold, qold as temporary variables.
  55 // Get rid of nontrivial GCD for {p, q, pnew, qnew}.
  56 void pqset_remove_gcd(pqset_ptr pq) {
  57   mpz_gcd(pq->pold, pq->pnew, pq->qnew);
  58   mpz_gcd(pq->qold, pq->p, pq->q);
  59   mpz_gcd(pq->pold, pq->pold, pq->qold);
  60   if (mpz_cmp_ui(pq->pold, 1)) {
  61     mpz_divexact(pq->pnew, pq->pnew, pq->pold);
  62     mpz_divexact(pq->qnew, pq->qnew, pq->pold);
  63     mpz_divexact(pq->p, pq->p, pq->pold);
  64     mpz_divexact(pq->q, pq->q, pq->pold);
  65   }
  66 }
  67
  68 void pqset_next(pqset_ptr pq) {
  69   mpz_set(pq->pold, pq->p); mpz_set(pq->p, pq->pnew);
  70   mpz_set(pq->qold, pq->q); mpz_set(pq->q, pq->qnew);
  71 }
  72
  73 struct mobius_data_s {
  74   cf_t input;
  75   mpz_t a, b, c, d;
  76 };
  77 typedef struct mobius_data_s *mobius_data_ptr;
  78
  79 void pqset_set_mobius(pqset_t pq, mobius_data_ptr md) {
  80   mpz_set(pq->pold, md->b); mpz_set(pq->p, md->a);
  81   mpz_set(pq->qold, md->d); mpz_set(pq->q, md->c);
  82 }
  83
  84 static void *mobius(cf_t cf) {
  85   mobius_data_ptr md = cf_data(cf);
  86   cf_t input = md->input;
  87   pqset_t pq;
  88   pqset_init(pq);
  89   pqset_set_mobius(pq, md);
  90
  91   mpz_t denom;
  92   mpz_init(denom);
  93   while(cf_wait(cf)) {
  94     cf_get(denom, input);
  95     pqset_regular_recur(pq, denom);
  96
  97     cf_put(cf, pq->pnew);
  98     cf_put(cf, pq->qnew);
  99     pqset_next(pq);
 100   }
 101   mpz_clear(denom);
 102   pqset_clear(pq);
 103
 104   mpz_clear(md->a);
 105   mpz_clear(md->b);
 106   mpz_clear(md->c);
 107   mpz_clear(md->d);
 108   free(md);
 109   return NULL;
 110 }
 111 cf_t cf_new_mobius_convergent(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 112   mobius_data_ptr md = malloc(sizeof(*md));
 113   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 114   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 115   md->input = x;
 116   return cf_new(mobius, md);
 117 }
 118
 119 cf_t cf_new_convergent(cf_t x) {
 120   mpz_t one, zero;
 121   mpz_init(one); mpz_init(zero);
 122   mpz_set_ui(one, 1); mpz_set_ui(zero, 0);
 123   cf_t res = cf_new_mobius_convergent(x, one, zero, zero, one);
 124   mpz_clear(one); mpz_clear(zero);
 125   return res;
 126 }
 127
 128 static void *nonregular_mobius_convergent(cf_t cf) {
 129   mobius_data_ptr md = cf_data(cf);
 130   cf_t input = md->input;
 131   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 132   mpz_t num; mpz_init(num);
 133   mpz_t denom; mpz_init(denom);
 134   void recur() {
 135     pqset_nonregular_recur(pq, num, denom);
 136     pqset_remove_gcd(pq);
 137
 138     cf_put(cf, pq->pnew);
 139     cf_put(cf, pq->qnew);
 140     pqset_next(pq);
 141   }
 142   mpz_set_ui(num, 1);
 143   cf_get(denom, input);
 144   recur();
 145   while(cf_wait(cf)) {
 146     cf_get(num, input);
 147     cf_get(denom, input);
 148     recur();
 149   }
 150   mpz_clear(num);
 151   mpz_clear(denom);
 152   pqset_clear(pq);
 153
 154   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 155   free(md);
 156   return NULL;
 157 }
 158 cf_t cf_new_nonregular_mobius_convergent(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 159   mobius_data_ptr md = malloc(sizeof(*md));
 160   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 161   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 162   md->input = x;
 163   return cf_new(nonregular_mobius_convergent, md);
 164 }
 165
 166 static void *mobius_nonregular_throughput(cf_t cf) {
 167   mobius_data_ptr md = cf_data(cf);
 168   cf_t input = md->input;
 169   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 170   mpz_t num; mpz_init(num);
 171   mpz_t denom; mpz_init(denom);
 172   mpz_t t0, t1; mpz_init(t1); mpz_init(t0);
 173   int recur() {
 174     pqset_nonregular_recur(pq, num, denom);
 175     pqset_remove_gcd(pq);
 176     if (mpz_sgn(pq->q)) {
 177       mpz_fdiv_qr(pq->pold, t0, pq->p, pq->q);
 178       // mpz_fdiv_qr(qold, t1, pnew, qnew);
 179       mpz_mul(pq->qold, pq->pold, pq->qnew);
 180
 181       // TODO: What if we're comparing equals?
 182       if (mpz_cmp(pq->qold, pq->pnew) < 0) {
 183         mpz_add(pq->qold, pq->qold, pq->qnew);
 184         if (mpz_cmp(pq->qold, pq->pnew) > 0) {
 185           // Output continued fraction term.
 186           cf_put(cf, pq->pold);
 187           // Compute remainder of pnew/qnew.
 188           mpz_sub(t1, pq->qold, pq->pnew);
 189           mpz_sub(t1, pq->qnew, t1);
 190
 191           mpz_set(pq->pold, pq->q);
 192           mpz_set(pq->qold, t0);
 193           mpz_set(pq->p, pq->qnew);
 194           mpz_set(pq->q, t1);
 195           return 1;
 196         }
 197       }
 198     }
 199     pqset_next(pq);
 200     return 0;
 201   }
 202   mpz_set_ui(num, 1);
 203   cf_get(denom, input);
 204   recur();
 205   while(cf_wait(cf)) {
 206     do {
 207       cf_get(num, input);
 208       cf_get(denom, input);
 209     } while(!recur());
 210   }
 211   mpz_clear(num);
 212   mpz_clear(denom);
 213   pqset_clear(pq);
 214
 215   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 216   free(md);
 217   mpz_clear(t1); mpz_clear(t0);
 218   return NULL;
 219 }
 220
 221 cf_t cf_new_nonregular_to_cf(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 222   mobius_data_ptr md = malloc(sizeof(*md));
 223   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 224   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 225   md->input = x;
 226   return cf_new(mobius_nonregular_throughput, md);
 227 }
 228
 229 static void *mobius_decimal(cf_t cf) {
 230   mobius_data_ptr md = cf_data(cf);
 231   cf_t input = md->input;
 232   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 233   mpz_t denom; mpz_init(denom);
 234   mpz_t t0, t1; mpz_init(t1); mpz_init(t0);
 235   int recur() {
 236     pqset_regular_recur(pq, denom);
 237
 238     // If the denominator is zero, we can't do anything yet.
 239     if (mpz_sgn(pq->q)) {
 240       // The answer is one of {0, ..., 9}.
 241       /* Naive attempt to expoit this didn't work well:
 242       int i;
 243       mpz_set(t0, pq->q);
 244       for (i = 0; i <= 9; i++) {
 245         if (mpz_cmp(t0, pq->p) > 0) break;
 246         mpz_add(t0, t0, pq->q);
 247       }
 248       mpz_set_ui(pq->pold, i);
 249       mpz_sub(t0, t0, pq->p);
 250       mpz_sub(t0, pq->q, t0);
 251       mpz_mul(pq->qold, pq->pold, pq->qnew);
 252       */
 253
 254       mpz_fdiv_qr(pq->pold, t0, pq->p, pq->q);
 255       mpz_mul(pq->qold, pq->pold, pq->qnew);
 256       // TODO: What if we're comparing equals?
 257       if (mpz_cmp(pq->qold, pq->pnew) < 0) {
 258         mpz_add(pq->qold, pq->qold, pq->qnew);
 259         if (mpz_cmp(pq->qold, pq->pnew) > 0) {
 260           // Output a decimal digit.
 261           cf_put(cf, pq->pold);
 262           // Compute t1 = remainder of pnew/qnew.
 263           mpz_sub(t1, pq->qold, pq->pnew);
 264           mpz_sub(t1, pq->qnew, t1);
 265           // Multiply numerator by 10.
 266           mpz_mul_ui(pq->pold, t0, 10);
 267           mpz_set(pq->qold, pq->q);
 268           mpz_mul_ui(pq->p, t1, 10);
 269           mpz_set(pq->q, pq->qnew);
 270           return 1;
 271         }
 272       }
 273     }
 274     pqset_next(pq);
 275     return 0;
 276   }
 277
 278   while(cf_wait(cf)) {
 279     do {
 280       cf_get(denom, input);
 281     } while(!recur());
 282   }
 283   mpz_clear(denom);
 284   pqset_clear(pq);
 285   mpz_clear(t0); mpz_clear(t1);
 286   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 287   free(md);
 288   return NULL;
 289 }
 290 cf_t cf_new_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 291   mobius_data_ptr md = malloc(sizeof(*md));
 292   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 293   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 294   md->input = x;
 295   return cf_new(mobius_decimal, md);
 296 }
 297
 298 cf_t cf_new_cf_to_decimal(cf_t x) {
 299   mpz_t one, zero;
 300   mpz_init(one); mpz_init(zero);
 301   mpz_set_ui(one, 1); mpz_set_ui(zero, 0);
 302   cf_t res = cf_new_mobius_to_decimal(x, one, zero, zero, one);
 303   mpz_clear(one); mpz_clear(zero);
 304   return res;
 305 }
 306
 307 static void *nonregular_mobius_decimal(cf_t cf) {
 308   mobius_data_ptr md = cf_data(cf);
 309   cf_t input = md->input;
 310   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 311   mpz_t num; mpz_init(num);
 312   mpz_t denom; mpz_init(denom);
 313   mpz_t t0, t1; mpz_init(t1); mpz_init(t0);
 314   int recur() {
 315     pqset_nonregular_recur(pq, num, denom);
 316     pqset_remove_gcd(pq);
 317
 318     // If the denominator is zero, we can't do anything yet.
 319     if (mpz_sgn(pq->q)) {
 320       mpz_fdiv_qr(pq->pold, t0, pq->p, pq->q);
 321       mpz_mul(pq->qold, pq->pold, pq->qnew);
 322       // TODO: What if we're comparing equals?
 323       if (mpz_cmp(pq->qold, pq->pnew) < 0) {
 324         mpz_add(pq->qold, pq->qold, pq->qnew);
 325         if (mpz_cmp(pq->qold, pq->pnew) > 0) {
 326           // Output a decimal digit.
 327           cf_put(cf, pq->pold);
 328           // Compute t1 = remainder of pnew/qnew.
 329           mpz_sub(t1, pq->qold, pq->pnew);
 330           mpz_sub(t1, pq->qnew, t1);
 331           // Multiply numerator by 10.
 332           mpz_mul_ui(pq->pold, t0, 10);
 333           mpz_set(pq->qold, pq->q);
 334           mpz_mul_ui(pq->p, t1, 10);
 335           mpz_set(pq->q, pq->qnew);
 336           return 1;
 337         }
 338       }
 339     }
 340     pqset_next(pq);
 341     return 0;
 342   }
 343   mpz_set_ui(num, 1);
 344   cf_get(denom, input);
 345   recur();
 346   while(cf_wait(cf)) {
 347     do {
 348       cf_get(num, input);
 349       cf_get(denom, input);
 350     } while(!recur());
 351   }
 352   mpz_clear(num);
 353   mpz_clear(denom);
 354   pqset_clear(pq);
 355   mpz_clear(t0); mpz_clear(t1);
 356   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 357   free(md);
 358   return NULL;
 359 }
 360
 361 cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 362   mobius_data_ptr md = malloc(sizeof(*md));
 363   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 364   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 365   md->input = x;
 366   return cf_new(nonregular_mobius_decimal, md);
 367 }