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, pq->p);
 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, mpz_t b, mpz_t c, mpz_t d) {
 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); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 295   md->input = x;
 296   return cf_new(nonregular_mobius_decimal, md);
 297 }
 298
 299 static void determine_sign(cf_t cf, pqset_t pq, mpz_t denom, cf_t input) {
 300   cf_set_sign(cf, cf_sign(input));
 301   while (mpz_sgn(pq->pold) != mpz_sgn(pq->p)
 302       || mpz_sgn(pq->qold) != mpz_sgn(pq->q)) {
 303     cf_get(denom, input);
 304     pqset_regular_recur(pq, denom);
 305   }
 306   if (mpz_sgn(pq->qold) < 0) {
 307     mpz_neg(pq->qold, pq->qold);
 308     mpz_neg(pq->q, pq->q);
 309     cf_flip_sign(cf);
 310   }
 311   if (mpz_sgn(pq->pold) < 0) {
 312     mpz_neg(pq->pold, pq->pold);
 313     mpz_neg(pq->p, pq->p);
 314     cf_flip_sign(cf);
 315   }
 316 }
 317
 318 // Input: Mobius transformation and regular continued fraction.
 319 // Output: Regular continued fraction.
 320 static void *mobius_throughput(cf_t cf) {
 321   mobius_data_ptr md = cf_data(cf);
 322   cf_t input = md->input;
 323   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 324   mpz_t denom; mpz_init(denom);
 325   mpz_t t0, t1, t2; mpz_init(t2); mpz_init(t1); mpz_init(t0);
 326
 327   determine_sign(cf, pq, denom, input);
 328
 329   int recur() {
 330     pqset_regular_recur(pq, denom);
 331
 332     if (mpz_sgn(pq->qold)) {
 333       mpz_fdiv_qr(t1, t0, pq->pold, pq->qold);
 334       mpz_mul(t2, t1, pq->q);
 335
 336       if (mpz_cmp(t2, pq->p) <= 0) {
 337         mpz_add(t2, t2, pq->q);
 338         if (mpz_cmp(t2, pq->p) > 0) {
 339           // Output continued fraction term.
 340           cf_put(cf, t1);
 341           // Subtract: remainder of p/q.
 342           mpz_sub(t2, t2, pq->p);
 343           mpz_sub(t2, pq->q, t2);
 344           // Invert
 345           mpz_set(pq->pold, pq->qold);
 346           mpz_set(pq->qold, t0);
 347           mpz_set(pq->p, pq->q);
 348           mpz_set(pq->q, t2);
 349           return 1;
 350         }
 351       }
 352     }
 353     return 0;
 354   }
 355   while(cf_wait(cf)) {
 356     do {
 357       cf_get(denom, input);
 358     } while(!recur());
 359   }
 360   mpz_clear(denom);
 361   pqset_clear(pq);
 362   mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
 363   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 364   free(md);
 365   return NULL;
 366 }
 367
 368 cf_t cf_new_mobius_to_cf(cf_t x, mpz_t z[4]) {
 369   mobius_data_ptr md = malloc(sizeof(*md));
 370   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 371   mpz_set(md->a, z[0]); mpz_set(md->b, z[1]);
 372   mpz_set(md->c, z[2]); mpz_set(md->d, z[3]);
 373   md->input = x;
 374   return cf_new(mobius_throughput, md);
 375 }
 376
 377 // Input: Mobius transformation and regular continued fraction.
 378 // Output: Decimal representation. The integer part is given first,
 379 // followed by one digit at a time.
 380 static void *mobius_decimal(cf_t cf) {
 381   mobius_data_ptr md = cf_data(cf);
 382   cf_t input = md->input;
 383   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
 384   mpz_t denom; mpz_init(denom);
 385   mpz_t t0, t1, t2; mpz_init(t2); mpz_init(t1); mpz_init(t0);
 386
 387   determine_sign(cf, pq, denom, input);
 388   int recur() {
 389     pqset_regular_recur(pq, denom);
 390
 391     // If the denominator is zero, we can't do anything yet.
 392     if (mpz_sgn(pq->qold)) {
 393       // Each term except possibly the first is one of {0, ..., 9}.
 394       /* Naive attempt to expoit this didn't seem faster:
 395        * (and I'd have to handle the first term properly)
 396       int i;
 397       mpz_set(t0, pq->q);
 398       for (i = 0; i <= 9; i++) {
 399         if (mpz_cmp(t0, pq->p) > 0) break;
 400         mpz_add(t0, t0, pq->q);
 401       }
 402       mpz_set_ui(pq->pold, i);
 403       mpz_sub(t0, t0, pq->p);
 404       mpz_sub(t0, pq->q, t0);
 405       mpz_mul(pq->qold, pq->pold, pq->qnew);
 406       */
 407
 408       mpz_fdiv_qr(t1, t0, pq->pold, pq->qold);
 409       mpz_mul(t2, t1, pq->q);
 410       if (mpz_cmp(t2, pq->p) <= 0) {
 411         mpz_add(t2, t2, pq->q);
 412         if (mpz_cmp(t2, pq->p) > 0) {
 413           // Output a decimal digit.
 414           cf_put(cf, t1);
 415           // Compute t2 = remainder of p/q.
 416           mpz_sub(t2, t2, pq->p);
 417           mpz_sub(t2, pq->q, t2);
 418           // Multiply numerator by 10.
 419           mpz_mul_ui(pq->pold, t0, 10);
 420           mpz_mul_ui(pq->p, t2, 10);
 421           return 1;
 422         }
 423       }
 424     }
 425     return 0;
 426   }
 427   while(cf_wait(cf)) {
 428     do {
 429       cf_get(denom, input);
 430     } while(!recur());
 431   }
 432   mpz_clear(denom);
 433   pqset_clear(pq);
 434   mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
 435   mpz_clear(md->a); mpz_clear(md->b); mpz_clear(md->c); mpz_clear(md->d);
 436   free(md);
 437   return NULL;
 438 }
 439
 440 cf_t cf_new_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) {
 441   mobius_data_ptr md = malloc(sizeof(*md));
 442   mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d);
 443   mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d);
 444   md->input = x;
 445   return cf_new(mobius_decimal, md);
 446 }
 447
 448 cf_t cf_new_cf_to_decimal(cf_t x) {
 449   mpz_t one, zero;
 450   mpz_init(one); mpz_init(zero);
 451   mpz_set_ui(one, 1); mpz_set_ui(zero, 0);
 452   cf_t res = cf_new_mobius_to_decimal(x, one, zero, zero, one);
 453   mpz_clear(one); mpz_clear(zero);
 454   return res;
 455 }