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