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_cf_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 // Input: Mobius transformation and nonregular continued fraction.
 177 // Output: Regular continued fraction. Assumes input fraction is well-behaved.
 178 static void *mobius_nonregular_throughput(cf_t cf) {
 179   mobius_data_ptr md = cf_data(cf);
 180   cf_t input = md->input;
 181   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 182   mpz_t num; mpz_init(num);
 183   mpz_t denom; mpz_init(denom);
 184   mpz_t t0, t1, t2; mpz_init(t2); mpz_init(t1); mpz_init(t0);
 185   int recur() {
 186     pqset_nonregular_recur(pq, num, denom);
 187     pqset_remove_gcd(pq, t0, t1);
 188
 189     if (mpz_sgn(pq->qold)) {
 190       mpz_fdiv_qr(t1, t0, pq->pold, pq->qold);
 191       mpz_mul(t2, t1, pq->q);
 192
 193       if (mpz_cmp(t2, pq->p) <= 0) {
 194         mpz_add(t2, t2, pq->q);
 195         if (mpz_cmp(t2, pq->p) > 0) {
 196           // Output continued fraction term.
 197           cf_put(cf, t1);
 198           // Subtract: remainder of p/q.
 199           mpz_sub(t2, t2, pq->p);
 200           mpz_sub(t2, pq->q, t2);
 201           // Invert
 202           mpz_set(pq->pold, pq->qold);
 203           mpz_set(pq->qold, t0);
 204           mpz_set(pq->p, pq->q);
 205           mpz_set(pq->q, t2);
 206           return 1;
 207         }
 208       }
 209     }
 210     return 0;
 211   }
 212   mpz_set_ui(num, 1);
 213   cf_get(denom, input);
 214   recur();
 215   while(cf_wait(cf)) {
 216     do {
 217       cf_get(num, input);
 218       cf_get(denom, input);
 219     } while(!recur());
 220   }
 221   mpz_clear(num);
 222   mpz_clear(denom);
 223   pqset_clear(pq);
 224
 225   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 226   free(md);
 227   mpz_clear(t2); mpz_clear(t1); mpz_clear(t0);
 228   return NULL;
 229 }
 230
 231 cf_t cf_new_nonregular_to_cf(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 232   mobius_data_ptr md = malloc(sizeof(*md));
 233   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 234   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 235   md->input = x;
 236   return cf_new(mobius_nonregular_throughput, md);
 237 }
 238
 239 // This seems to be slower than regularizing the continued fraction
 240 // and then converting to decimal.
 241 static void *nonregular_mobius_decimal(cf_t cf) {
 242   mobius_data_ptr md = cf_data(cf);
 243   cf_t input = md->input;
 244   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 245   mpz_t num; mpz_init(num);
 246   mpz_t denom; mpz_init(denom);
 247   mpz_t t0, t1, t2; mpz_init(t2); mpz_init(t1); mpz_init(t0);
 248   int recur() {
 249     pqset_nonregular_recur(pq, num, denom);
 250     pqset_remove_gcd(pq, t0, t1);
 251
 252     // If the denominator is zero, we can't do anything yet.
 253     if (mpz_sgn(pq->qold)) {
 254       mpz_fdiv_qr(t1, t0, pq->pold, pq->qold);
 255       mpz_mul(t2, t1, pq->q);
 256       if (mpz_cmp(t2, pq->p) <= 0) {
 257         mpz_add(t2, t2, pq->q);
 258         if (mpz_cmp(t2, pq->p) > 0) {
 259           // Output a decimal digit.
 260           cf_put(cf, t1);
 261           // Subtract: remainder of p/q.
 262           mpz_sub(t2, t2, pq->p);
 263           mpz_sub(t2, pq->q, t2);
 264           // Multiply numerator by 10.
 265           mpz_mul_ui(pq->pold, t0, 10);
 266           mpz_mul_ui(pq->p, t2, 10);
 267           return 1;
 268         }
 269       }
 270     }
 271     return 0;
 272   }
 273   mpz_set_ui(num, 1);
 274   cf_get(denom, input);
 275   recur();
 276   while(cf_wait(cf)) {
 277     do {
 278       cf_get(num, input);
 279       cf_get(denom, input);
 280     } while(!recur());
 281   }
 282   mpz_clear(num);
 283   mpz_clear(denom);
 284   pqset_clear(pq);
 285   mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
 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
 291 cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a[4]) {
 292   mobius_data_ptr md = malloc(sizeof(*md));
 293   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 294   mpz_set(md->a, a[0]); mpz_set(md->b, a[1]);
 295   mpz_set(md->c, a[2]); mpz_set(md->d, a[3]);
 296   md->input = x;
 297   return cf_new(nonregular_mobius_decimal, md);
 298 }
 299
 300 static void determine_sign(cf_t cf, pqset_t pq, mpz_t denom, cf_t input) {
 301   cf_set_sign(cf, cf_sign(input));
 302   while (mpz_sgn(pq->pold) != mpz_sgn(pq->p)
 303       || mpz_sgn(pq->qold) != mpz_sgn(pq->q)) {
 304     cf_get(denom, input);
 305     pqset_regular_recur(pq, denom);
 306   }
 307   if (mpz_sgn(pq->qold) < 0) {
 308     mpz_neg(pq->qold, pq->qold);
 309     mpz_neg(pq->q, pq->q);
 310     cf_flip_sign(cf);
 311   }
 312   if (mpz_sgn(pq->pold) < 0) {
 313     mpz_neg(pq->pold, pq->pold);
 314     mpz_neg(pq->p, pq->p);
 315     cf_flip_sign(cf);
 316   }
 317 }
 318
 319 // Input: Mobius transformation and regular continued fraction.
 320 // Output: Regular continued fraction.
 321 static void *mobius_throughput(cf_t cf) {
 322   mobius_data_ptr md = cf_data(cf);
 323   cf_t input = md->input;
 324   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 325   mpz_t denom; mpz_init(denom);
 326   mpz_t t0, t1, t2; mpz_init(t2); mpz_init(t1); mpz_init(t0);
 327
 328   determine_sign(cf, pq, denom, input);
 329
 330   int recur() {
 331     pqset_regular_recur(pq, denom);
 332
 333     if (mpz_sgn(pq->qold)) {
 334       mpz_fdiv_qr(t1, t0, pq->pold, pq->qold);
 335       mpz_mul(t2, t1, pq->q);
 336
 337       if (mpz_cmp(t2, pq->p) <= 0) {
 338         mpz_add(t2, t2, pq->q);
 339         if (mpz_cmp(t2, pq->p) > 0) {
 340           // Output continued fraction term.
 341           cf_put(cf, t1);
 342           // Subtract: remainder of p/q.
 343           mpz_sub(t2, t2, pq->p);
 344           mpz_sub(t2, pq->q, t2);
 345           // Invert
 346           mpz_set(pq->pold, pq->qold);
 347           mpz_set(pq->qold, t0);
 348           mpz_set(pq->p, pq->q);
 349           mpz_set(pq->q, t2);
 350           return 1;
 351         }
 352       }
 353     }
 354     return 0;
 355   }
 356   while(cf_wait(cf)) {
 357     do {
 358       cf_get(denom, input);
 359     } while(!recur());
 360   }
 361   mpz_clear(denom);
 362   pqset_clear(pq);
 363   mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
 364   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 365   free(md);
 366   return NULL;
 367 }
 368
 369 cf_t cf_new_mobius_to_cf(cf_t x, mpz_t z[4]) {
 370   mobius_data_ptr md = malloc(sizeof(*md));
 371   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 372   mpz_set(md->a, z[0]); mpz_set(md->b, z[1]);
 373   mpz_set(md->c, z[2]); mpz_set(md->d, z[3]);
 374   md->input = x;
 375   return cf_new(mobius_throughput, md);
 376 }
 377
 378 // Input: Mobius transformation and regular continued fraction.
 379 // Output: Decimal representation. The integer part is given first,
 380 // followed by one digit at a time.
 381 static void *mobius_decimal(cf_t cf) {
 382   mobius_data_ptr md = cf_data(cf);
 383   cf_t input = md->input;
 384   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 385   mpz_t denom; mpz_init(denom);
 386   mpz_t t0, t1, t2; mpz_init(t2); mpz_init(t1); mpz_init(t0);
 387
 388   determine_sign(cf, pq, denom, input);
 389   int recur() {
 390     pqset_regular_recur(pq, denom);
 391
 392     // If the denominator is zero, we can't do anything yet.
 393     if (mpz_sgn(pq->qold)) {
 394       // Each term except possibly the first is one of {0, ..., 9}.
 395       /* Naive attempt to expoit this didn't seem faster:
 396        * (and I'd have to handle the first term properly)
 397       int i;
 398       mpz_set(t0, pq->q);
 399       for (i = 0; i <= 9; i++) {
 400         if (mpz_cmp(t0, pq->p) > 0) break;
 401         mpz_add(t0, t0, pq->q);
 402       }
 403       mpz_set_ui(pq->pold, i);
 404       mpz_sub(t0, t0, pq->p);
 405       mpz_sub(t0, pq->q, t0);
 406       mpz_mul(pq->qold, pq->pold, pq->qnew);
 407       */
 408
 409       mpz_fdiv_qr(t1, t0, pq->pold, pq->qold);
 410       mpz_mul(t2, t1, pq->q);
 411       if (mpz_cmp(t2, pq->p) <= 0) {
 412         mpz_add(t2, t2, pq->q);
 413         if (mpz_cmp(t2, pq->p) > 0) {
 414           // Output a decimal digit.
 415           cf_put(cf, t1);
 416           // Compute t2 = remainder of p/q.
 417           mpz_sub(t2, t2, pq->p);
 418           mpz_sub(t2, pq->q, t2);
 419           // Multiply numerator by 10.
 420           mpz_mul_ui(pq->pold, t0, 10);
 421           mpz_mul_ui(pq->p, t2, 10);
 422           return 1;
 423         }
 424       }
 425     }
 426     return 0;
 427   }
 428   while(cf_wait(cf)) {
 429     do {
 430       cf_get(denom, input);
 431     } while(!recur());
 432   }
 433   mpz_clear(denom);
 434   pqset_clear(pq);
 435   mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
 436   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 437   free(md);
 438   return NULL;
 439 }
 440
 441 cf_t cf_new_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 442   mobius_data_ptr md = malloc(sizeof(*md));
 443   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 444   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 445   md->input = x;
 446   return cf_new(mobius_decimal, md);
 447 }
 448
 449 cf_t cf_new_cf_to_decimal(cf_t x) {
 450   mpz_t one, zero;
 451   mpz_init(one); mpz_init(zero);
 452   mpz_set_ui(one, 1); mpz_set_ui(zero, 0);
 453   cf_t res = cf_new_mobius_to_decimal(x, one, zero, zero, one);
 454   mpz_clear(one); mpz_clear(zero);
 455   return res;
 456 }
 457
 458 struct funarg_s {
 459   void *(*fun)(cf_t);
 460   mpz_t arg;
 461 };
 462 typedef struct funarg_s *funarg_ptr;
 463
 464 static void *one_arg_nonregular(cf_t cf) {
 465   funarg_ptr p = cf_data(cf);
 466   mpz_t a, b, c, d;
 467   mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
 468   mpz_set_ui(a, 1); mpz_set_ui(b, 0);
 469   mpz_set_ui(c, 0); mpz_set_ui(d, 1);
 470   mpz_ptr copy = malloc(sizeof(*copy));
 471   mpz_init(copy);
 472   mpz_set(copy, p->arg);
 473   cf_t nonregular = cf_new(p->fun, copy);
 474   cf_t conv = cf_new_nonregular_to_cf(nonregular, a, b, c, d);
 475   mpz_t z;
 476   mpz_init(z);
 477   while(cf_wait(cf)) {
 478     cf_get(z, conv);
 479     cf_put(cf, z);
 480   }
 481   mpz_clear(z);
 482   cf_free(conv);
 483   cf_free(nonregular);
 484   mpz_clear(a); mpz_clear(b); mpz_clear(c); mpz_clear(d);
 485
 486   mpz_clear(p->arg);
 487   free(p);
 488   return NULL;
 489 }
 490
 491 cf_t cf_new_one_arg(void *(*fun)(cf_t), mpz_t z) {
 492   mpz_ptr p = malloc(sizeof(*p));
 493   mpz_init(p);
 494   mpz_set(p, z);
 495   return cf_new(fun, p);
 496 }
 497
 498 cf_t cf_new_one_arg_nonregular(void *(*fun)(cf_t), mpz_t z) {
 499   funarg_ptr p = malloc(sizeof(*p));
 500   p->fun = fun;
 501   mpz_init(p->arg);
 502   mpz_set(p->arg, z);
 503   return cf_new(one_arg_nonregular, p);
 504 }
 505
 506 static void *nonreg_const(cf_t cf) {
 507   funarg_ptr p = cf_data(cf);
 508   mpz_t a, b, c, d;
 509   mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
 510   mpz_set_ui(a, 1); mpz_set_ui(b, 0);
 511   mpz_set_ui(c, 0); mpz_set_ui(d, 1);
 512   cf_t nonregular = cf_new(p->fun, NULL);
 513   cf_t conv = cf_new_nonregular_to_cf(nonregular, a, b, c, d);
 514   mpz_clear(a); mpz_clear(b); mpz_clear(c); mpz_clear(d);
 515   mpz_t z;
 516   mpz_init(z);
 517   while(cf_wait(cf)) {
 518     cf_get(z, conv);
 519     cf_put(cf, z);
 520   }
 521   mpz_clear(z);
 522   cf_free(conv);
 523   cf_free(nonregular);
 524   free(p);
 525   return NULL;
 526 }
 527
 528 cf_t cf_new_const_nonregular(void *(*fun)(cf_t)) {
 529   funarg_ptr p = malloc(sizeof(*p));
 530   p->fun = fun;
 531   return cf_new(nonreg_const, p);
 532 }