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   return cf_new_mobius_convergent(x, one, zero, zero, one);
 124   mpz_clear(one); mpz_clear(zero);
 125 }
 126
 127 static void *nonregular_mobius_convergent(cf_t cf) {
 128   mobius_data_ptr md = cf_data(cf);
 129   cf_t input = md->input;
 130   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 131   mpz_t num; mpz_init(num);
 132   mpz_t denom; mpz_init(denom);
 133   void recur() {
 134     pqset_nonregular_recur(pq, num, denom);
 135     pqset_remove_gcd(pq);
 136
 137     cf_put(cf, pq->pnew);
 138     cf_put(cf, pq->qnew);
 139     pqset_next(pq);
 140   }
 141   mpz_set_ui(num, 1);
 142   cf_get(denom, input);
 143   recur();
 144   while(cf_wait(cf)) {
 145     cf_get(num, input);
 146     cf_get(denom, input);
 147     recur();
 148   }
 149   mpz_clear(num);
 150   mpz_clear(denom);
 151   pqset_clear(pq);
 152
 153   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 154   free(md);
 155   return NULL;
 156 }
 157 cf_t cf_new_nonregular_mobius_convergent(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 158   mobius_data_ptr md = malloc(sizeof(*md));
 159   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 160   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 161   md->input = x;
 162   return cf_new(nonregular_mobius_convergent, md);
 163 }
 164
 165 static void *mobius_nonregular_throughput(cf_t cf) {
 166   mobius_data_ptr md = cf_data(cf);
 167   cf_t input = md->input;
 168   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 169   mpz_t num; mpz_init(num);
 170   mpz_t denom; mpz_init(denom);
 171   mpz_t t0, t1; mpz_init(t1); mpz_init(t0);
 172   int recur() {
 173     pqset_nonregular_recur(pq, num, denom);
 174     pqset_remove_gcd(pq);
 175     if (mpz_sgn(pq->q)) {
 176       mpz_fdiv_qr(pq->pold, t0, pq->p, pq->q);
 177       // mpz_fdiv_qr(qold, t1, pnew, qnew);
 178       mpz_mul(pq->qold, pq->pold, pq->qnew);
 179
 180       // TODO: What if we're comparing equals?
 181       if (mpz_cmp(pq->qold, pq->pnew) < 0) {
 182         mpz_add(pq->qold, pq->qold, pq->qnew);
 183         if (mpz_cmp(pq->qold, pq->pnew) > 0) {
 184           // Output continued fraction term.
 185           cf_put(cf, pq->pold);
 186           // Compute remainder of pnew/qnew.
 187           mpz_sub(t1, pq->qold, pq->pnew);
 188           mpz_sub(t1, pq->qnew, t1);
 189
 190           mpz_set(pq->pold, pq->q);
 191           mpz_set(pq->qold, t0);
 192           mpz_set(pq->p, pq->qnew);
 193           mpz_set(pq->q, t1);
 194           return 1;
 195         }
 196       }
 197     }
 198     pqset_next(pq);
 199     return 0;
 200   }
 201   mpz_set_ui(num, 1);
 202   cf_get(denom, input);
 203   recur();
 204   while(cf_wait(cf)) {
 205     do {
 206       cf_get(num, input);
 207       cf_get(denom, input);
 208     } while(!recur());
 209   }
 210   mpz_clear(denom);
 211   pqset_clear(pq);
 212
 213   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 214   free(md);
 215   mpz_clear(t1); mpz_clear(t0);
 216   return NULL;
 217 }
 218 cf_t cf_new_nonregular_to_cf(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 219   mobius_data_ptr md = malloc(sizeof(*md));
 220   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 221   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 222   md->input = x;
 223   return cf_new(mobius_nonregular_throughput, md);
 224 }
 225
 226 static void *mobius_decimal(cf_t cf) {
 227   mobius_data_ptr md = cf_data(cf);
 228   cf_t input = md->input;
 229   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 230   mpz_t denom; mpz_init(denom);
 231   mpz_t t0, t1; mpz_init(t1); mpz_init(t0);
 232   int recur() {
 233     pqset_regular_recur(pq, denom);
 234
 235     // If the denominator is zero, we can't do anything yet.
 236     if (mpz_sgn(pq->q)) {
 237       // The answer is one of {0, ..., 9}.
 238       /* Naive attempt to expoit this didn't work well:
 239       int i;
 240       mpz_set(t0, pq->q);
 241       for (i = 0; i <= 9; i++) {
 242         if (mpz_cmp(t0, pq->p) > 0) break;
 243         mpz_add(t0, t0, pq->q);
 244       }
 245       mpz_set_ui(pq->pold, i);
 246       mpz_sub(t0, t0, pq->p);
 247       mpz_sub(t0, pq->q, t0);
 248       mpz_mul(pq->qold, pq->pold, pq->qnew);
 249       */
 250
 251       mpz_fdiv_qr(pq->pold, t0, pq->p, pq->q);
 252       mpz_mul(pq->qold, pq->pold, pq->qnew);
 253       // TODO: What if we're comparing equals?
 254       if (mpz_cmp(pq->qold, pq->pnew) < 0) {
 255         mpz_add(pq->qold, pq->qold, pq->qnew);
 256         if (mpz_cmp(pq->qold, pq->pnew) > 0) {
 257           // Output a decimal digit.
 258           cf_put(cf, pq->pold);
 259           // Compute t1 = remainder of pnew/qnew.
 260           mpz_sub(t1, pq->qold, pq->pnew);
 261           mpz_sub(t1, pq->qnew, t1);
 262           // Multiply numerator by 10.
 263           mpz_mul_ui(pq->pold, t0, 10);
 264           mpz_set(pq->qold, pq->q);
 265           mpz_mul_ui(pq->p, t1, 10);
 266           mpz_set(pq->q, pq->qnew);
 267           return 1;
 268         }
 269       }
 270     }
 271     pqset_next(pq);
 272     return 0;
 273   }
 274
 275   while(cf_wait(cf)) {
 276     do {
 277       cf_get(denom, input);
 278     } while(!recur());
 279   }
 280   mpz_clear(denom);
 281   pqset_clear(pq);
 282   mpz_clear(t0); mpz_clear(t1);
 283   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 284   free(md);
 285   return NULL;
 286 }
 287 cf_t cf_new_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 288   mobius_data_ptr md = malloc(sizeof(*md));
 289   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 290   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 291   md->input = x;
 292   return cf_new(mobius_decimal, md);
 293 }
 294
 295 cf_t cf_new_cf_to_decimal(cf_t x) {
 296   mpz_t one, zero;
 297   mpz_init(one); mpz_init(zero);
 298   mpz_set_ui(one, 1); mpz_set_ui(zero, 0);
 299   return cf_new_mobius_to_decimal(x, one, zero, zero, one);
 300   mpz_clear(one); mpz_clear(zero);
 301 }
 302
 303 static void *nonregular_mobius_decimal(cf_t cf) {
 304   mobius_data_ptr md = cf_data(cf);
 305   cf_t input = md->input;
 306   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 307   mpz_t num; mpz_init(num);
 308   mpz_t denom; mpz_init(denom);
 309   mpz_t t0, t1; mpz_init(t1); mpz_init(t0);
 310   int recur() {
 311     pqset_nonregular_recur(pq, num, denom);
 312     pqset_remove_gcd(pq);
 313
 314     // If the denominator is zero, we can't do anything yet.
 315     if (mpz_sgn(pq->q)) {
 316       mpz_fdiv_qr(pq->pold, t0, pq->p, pq->q);
 317       mpz_mul(pq->qold, pq->pold, pq->qnew);
 318       // TODO: What if we're comparing equals?
 319       if (mpz_cmp(pq->qold, pq->pnew) < 0) {
 320         mpz_add(pq->qold, pq->qold, pq->qnew);
 321         if (mpz_cmp(pq->qold, pq->pnew) > 0) {
 322           // Output a decimal digit.
 323           cf_put(cf, pq->pold);
 324           // Compute t1 = remainder of pnew/qnew.
 325           mpz_sub(t1, pq->qold, pq->pnew);
 326           mpz_sub(t1, pq->qnew, t1);
 327           // Multiply numerator by 10.
 328           mpz_mul_ui(pq->pold, t0, 10);
 329           mpz_set(pq->qold, pq->q);
 330           mpz_mul_ui(pq->p, t1, 10);
 331           mpz_set(pq->q, pq->qnew);
 332           return 1;
 333         }
 334       }
 335     }
 336     pqset_next(pq);
 337     return 0;
 338   }
 339   mpz_set_ui(num, 1);
 340   cf_get(denom, input);
 341   recur();
 342   while(cf_wait(cf)) {
 343     do {
 344       cf_get(num, input);
 345       cf_get(denom, input);
 346     } while(!recur());
 347   }
 348   mpz_clear(num);
 349   mpz_clear(denom);
 350   pqset_clear(pq);
 351   mpz_clear(t0); mpz_clear(t1);
 352   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 353   free(md);
 354   return NULL;
 355 }
 356
 357 cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 358   mobius_data_ptr md = malloc(sizeof(*md));
 359   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 360   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 361   md->input = x;
 362   return cf_new(nonregular_mobius_decimal, md);
 363 }