newton.c

   1 // Solve quadratic equations involving continued fractions.
   2 // A slight misnomer: Gosper describes how to use Newton's method, but
   3 // I use binary search instead, and assume unique roots in the right ranges.
   4
   5 #include <stdio.h>
   6 #include <stdlib.h>
   7 #include <gmp.h>
   8 #include "cf.h"
   9
  10 // Represents:
  11 //
  12 // a0 xy + a1 x + a2 y + a3
  13 // ------------------------
  14 // a4 xy - a0 x + a5 y - a2
  15
  16 struct newton_data_s {
  17   cf_t x;
  18   mpz_t a[6];
  19   mpz_t lower;
  20 };
  21 typedef struct newton_data_s newton_data_t[1];
  22 typedef struct newton_data_s *newton_data_ptr;
  23
  24 // In the 3D table, we have two sets of equations. Left is a0, b0, c0,
  25 // Right is a1, b1, c1, and these are the terms involving x:
  26 //
  27 //   a0 y + b0     a1 y + b1
  28 //   ---------     ---------
  29 //   c0 y - a0     c1 y - a1
  30 //
  31 // I work left to right, and top to bottom, unlike Gosper.
  32 //
  33 // This notation is different to that in cf_bihom.c, and similar to cf_mobius.
  34
  35 struct abc_s {
  36   mpz_t a0, a1;
  37   mpz_t b0, b1;
  38   mpz_t c0, c1;
  39 };
  40 typedef struct abc_s abc_t[1];
  41
  42 void abc_init(abc_t p) {
  43   mpz_init(p->a0); mpz_init(p->a1);
  44   mpz_init(p->b0); mpz_init(p->b1);
  45   mpz_init(p->c0); mpz_init(p->c1);
  46 }
  47
  48 void abc_clear(abc_t p) {
  49   mpz_clear(p->a0); mpz_clear(p->a1);
  50   mpz_clear(p->b0); mpz_clear(p->b1);
  51   mpz_clear(p->c0); mpz_clear(p->c1);
  52 }
  53
  54 void abc_set_coeff(abc_t p, mpz_t a[6]) {
  55   mpz_set(p->a0, a[2]); mpz_set(p->b0, a[3]); mpz_set(p->c0, a[5]);
  56   mpz_set(p->a1, a[0]); mpz_set(p->b1, a[1]); mpz_set(p->c1, a[4]);
  57 }
  58
  59 void abc_print(abc_t p) {
  60   gmp_printf("%Zd y + %Zd    %Zd y + %Zd\n", p->a0, p->b0, p->a1, p->b1);
  61   gmp_printf("%Zd y - %Zd    %Zd y - %Zd\n", p->c0, p->a0, p->c1, p->a1);
  62   gmp_printf("(%Zd+sqrt(%Zd^2+%Zd*%Zd))/%Zd\n", p->a0, p->a0, p->b0, p->c0, p->c0);
  63   gmp_printf("%Zd*z^2-2*%Zd*z-%Zd\n", p->c0, p->a0, p->b0);
  64 }
  65
  66 // Finds smallest root greater than given lower bound.
  67 // Assumes there exists exactly one root with this property.
  68 // (Impossible to solve equations if roots have same integer part at
  69 // the moment. Don't do that.)
  70 // TODO: Rewrite so you choose if you want smaller or greater root.
  71 // or so it returns both solutions in an array.
  72 static void *newton(cf_t cf) {
  73   newton_data_ptr nd = cf_data(cf);
  74   abc_t p;
  75   abc_init(p);
  76   abc_set_coeff(p, nd->a);
  77   cf_t x = nd->x;
  78   mpz_t z, z0, z1, one, pow2, t0, t1, t2;
  79   mpz_init(z); mpz_init(z0); mpz_init(z1); mpz_init(pow2);
  80   mpz_init(t0); mpz_init(t1); mpz_init(t2);
  81   mpz_init(one);
  82   mpz_set_ui(one, 1);
  83
  84   void move_right() {
  85     cf_get(z, x);
  86     mpz_mul(t0, z, p->b1);
  87     mpz_add(t0, t0, p->b0);
  88     mpz_set(p->b0, p->b1);
  89     mpz_set(p->b1, t0);
  90
  91     mpz_mul(t0, z, p->a1);  mpz_mul(t1, z, p->c1);
  92     mpz_add(t0, t0, p->a0); mpz_add(t1, t1, p->c0);
  93     mpz_set(p->a0, p->a1);  mpz_set(p->c0, p->c1);
  94     mpz_set(p->a1, t0);     mpz_set(p->c1, t1);
  95   }
  96
  97   void move_down() {
  98     // Recurrence relation with z, then
  99     // Subtract and invert z.
 100     mpz_mul(t0, z, p->a0);
 101     mpz_add(t0, t0, p->b0);
 102     mpz_mul(t1, z, p->c0);
 103     mpz_sub(t1, t1, p->a0);
 104     mpz_mul(p->a0, z, t1);
 105     mpz_set(p->b0, p->c0);
 106     mpz_sub(p->c0, t0, p->a0);
 107     mpz_set(p->a0, t1);
 108
 109     mpz_mul(t0, z, p->a1);
 110     mpz_add(t0, t0, p->b1);
 111     mpz_mul(t1, z, p->c1);
 112     mpz_sub(t1, t1, p->a1);
 113     mpz_mul(p->a1, z, t1);
 114     mpz_set(p->b1, p->c1);
 115     mpz_sub(p->c1, t0, p->a1);
 116     mpz_set(p->a1, t1);
 117   }
 118
 119   int sign_quad() {
 120     // Returns sign of c0 z^2 - 2 a0 z - b0
 121     mpz_mul_si(t0, p->a0, -2);
 122     mpz_addmul(t0, p->c0, z);
 123     mpz_mul(t0, t0, z);
 124     mpz_sub(t0, t0, p->b0);
 125     return mpz_sgn(t0);
 126   }
 127
 128   int sign_quad1() {
 129     // Returns sign of c1 z^2 - 2 a1 z - b1
 130     mpz_mul_si(t0, p->a1, -2);
 131     mpz_addmul(t0, p->c1, z);
 132     mpz_mul(t0, t0, z);
 133     mpz_sub(t0, t0, p->b1);
 134     return mpz_sgn(t0);
 135   }
 136
 137   move_right();  // Get rid of pathological cases.
 138   move_right();
 139   move_right();
 140
 141   // Get integer part, starting search from given lower bound.
 142   void binary_search(mpz_ptr lower) {
 143     while (!mpz_sgn(p->c0)) move_right();
 144     for (;;) {
 145       mpz_set(z0, lower);
 146       mpz_set(z, lower);
 147       int sign = sign_quad();
 148       mpz_set_ui(pow2, 1);
 149       for (;;) {
 150         mpz_add(z, z0, pow2);
 151         if (sign_quad() != sign) break;
 152         mpz_mul_2exp(pow2, pow2, 1);
 153       }
 154       mpz_set(z1, z);
 155
 156       for (;;) {
 157         mpz_add(z, z0, z1);
 158         mpz_div_2exp(z, z, 1);
 159         if (!mpz_cmp(z, z0)) break;
 160         if (sign_quad() == sign) {
 161           mpz_set(z0, z);
 162         } else {
 163           mpz_set(z1, z);
 164         }
 165       }
 166       sign = sign_quad1();
 167       mpz_set(z, z1);
 168       if (sign_quad1() != sign) break;
 169       move_right();
 170     }
 171   }
 172
 173   binary_search(nd->lower);
 174   cf_put(cf, z0);
 175   mpz_set(z, z0);
 176   move_down();
 177
 178   while(cf_wait(cf)) {
 179     binary_search(one);
 180     cf_put(cf, z0);
 181     mpz_set(z, z0);
 182     move_down();
 183   }
 184   abc_clear(p);
 185   mpz_clear(z); mpz_clear(z0); mpz_clear(z1); mpz_clear(pow2);
 186   mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
 187   for (int i = 0; i < 6; i++) mpz_clear(nd->a[i]);
 188   free(nd);
 189   return NULL;
 190 }
 191
 192 cf_t cf_new_newton(cf_t x, mpz_t a[6], mpz_t lower) {
 193   newton_data_ptr p = malloc(sizeof(*p));
 194   p->x = x;
 195   for (int i = 0; i < 6; i++) {
 196     mpz_init(p->a[i]);
 197     mpz_set(p->a[i], a[i]);
 198   }
 199   mpz_init(p->lower);
 200   mpz_set(p->lower, lower);
 201   return cf_new(newton, p);
 202 }
 203
 204 cf_t cf_new_sqrt(cf_t x) {
 205   newton_data_ptr p = malloc(sizeof(*p));
 206   p->x = x;
 207   mpz_init(p->lower);
 208   for (int i = 0; i < 6; i++) mpz_init(p->a[i]);
 209   // Solve y = x/y.
 210   mpz_set_si(p->a[1], 1);
 211   mpz_set_si(p->a[5], 1);
 212   return cf_new(newton, p);
 213 }