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