Compute HAKMEM constant.
authorBen Lynn <benlynn@gmail.com>
Thu, 18 Sep 2008 10:20:48 +0000 (18 03:20 -0700)
committerBen Lynn <benlynn@gmail.com>
Thu, 18 Sep 2008 10:22:54 +0000 (18 03:22 -0700)
Makefile
README
bihom.c
bihom_test.c
cf.h
famous.c
hakmem.c [new file with mode: 0644]
mobius.c
newton.c
newton_test.c
nonregular.c [new file with mode: 0644]

index 522fc9a..d8ad84f 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,7 @@
 
 CF_OBJS:=cf.o mobius.o famous.o bihom.o taylor.o test.o newton.o tee.o
 TESTS:=bihom_test cf_test famous_test mobius_test newton_test tee_test
-BINS:=pi epow
+BINS:=pi epow hakmem
 
 target : $(BINS)
 
@@ -10,11 +10,12 @@ libfrac.a : $(CF_OBJS)
        ar rvs $@ $^
 
 %.o : %.c
-       gcc -g -std=c99 -Wall -c -o $@ $<
+       gcc -O3 -std=c99 -Wall -c -o $@ $<
+       #gcc -g -std=c99 -Wall -c -o $@ $<
 
 % : %.c libfrac.a
-       #gcc -O3 -std=c99 -Wall -o $@ $< -lgmp -lpthread -lfrac -L .
-       gcc -g -std=c99 -Wall -o $@ $< -lgmp -lpthread -lfrac -L .
+       gcc -O3 -std=c99 -Wall -o $@ $< -lgmp -lpthread -lfrac -L .
+       #gcc -g -std=c99 -Wall -o $@ $< -lgmp -lpthread -lfrac -L .
 
 test: $(TESTS)
 
diff --git a/README b/README
index 9711093..bc0a2d6 100644 (file)
--- a/README
+++ b/README
@@ -62,7 +62,7 @@ fraction. It takes longer than I had hoped, but it's faster than
 
 There's no API documentation; see the source.
 
-== Demand channels ==
+== Design ==
 
 Threads are the future, if not already the present. Once, single core CPUs grew
 faster at an astounding pace. But limits are being approached, and engineers
diff --git a/bihom.c b/bihom.c
index 96c8942..b7007d1 100644 (file)
--- a/bihom.c
+++ b/bihom.c
@@ -185,3 +185,100 @@ cf_t cf_new_bihom(cf_t x, cf_t y, mpz_t a[8]) {
   }
   return cf_new(bihom, p);
 }
+
+void mpz8_init(mpz_t z[8]) {
+  int i;
+  for(i = 0; i < 8; i++) {
+    mpz_init(z[i]);
+  }
+}
+
+void mpz8_clear(mpz_t z[8]) {
+  int i;
+  for(i = 0; i < 8; i++) {
+    mpz_clear(z[i]);
+  }
+}
+
+void mpz8_set_int(mpz_t z[8],
+    int a, int b, int c, int d,
+    int e, int f, int g, int h) {
+  mpz_set_si(z[0], a);
+  mpz_set_si(z[1], b);
+  mpz_set_si(z[2], c);
+  mpz_set_si(z[3], d);
+  mpz_set_si(z[4], e);
+  mpz_set_si(z[5], f);
+  mpz_set_si(z[6], g);
+  mpz_set_si(z[7], h);
+}
+
+void mpz8_set_add(mpz_t z[8]) {
+  mpz8_set_int(z,
+    0, 1, 1, 0,
+    0, 0, 0, 1);
+}
+
+void mpz8_set_sub(mpz_t z[8]) {
+  mpz8_set_int(z,
+    0, 1, -1, 0,
+    0, 0, 0, 1);
+}
+
+void mpz8_set_mul(mpz_t z[8]) {
+  mpz8_set_int(z,
+    1, 0, 0, 0,
+    0, 0, 0, 1);
+}
+
+void mpz8_set_div(mpz_t z[8]) {
+  mpz8_set_int(z,
+    0, 1, 0, 0,
+    0, 0, 1, 0);
+}
+
+cf_t cf_new_add(cf_t x, cf_t y) {
+  bihom_data_ptr p = malloc(sizeof(*p));
+  p->x = x;
+  p->y = y;
+  for (int i = 0; i < 8; i++) mpz_init(p->a[i]);
+  mpz_set_si(p->a[1], 1);
+  mpz_set_si(p->a[2], 1);
+  mpz_set_si(p->a[7], 1);
+
+  return cf_new(bihom, p);
+}
+
+cf_t cf_new_sub(cf_t x, cf_t y) {
+  bihom_data_ptr p = malloc(sizeof(*p));
+  p->x = x;
+  p->y = y;
+  for (int i = 0; i < 8; i++) mpz_init(p->a[i]);
+  mpz_set_si(p->a[1], 1);
+  mpz_set_si(p->a[2], -1);
+  mpz_set_si(p->a[7], 1);
+
+  return cf_new(bihom, p);
+}
+
+cf_t cf_new_mul(cf_t x, cf_t y) {
+  bihom_data_ptr p = malloc(sizeof(*p));
+  p->x = x;
+  p->y = y;
+  for (int i = 0; i < 8; i++) mpz_init(p->a[i]);
+  mpz_set_si(p->a[0], 1);
+  mpz_set_si(p->a[7], 1);
+
+  return cf_new(bihom, p);
+}
+
+cf_t cf_new_div(cf_t x, cf_t y) {
+  bihom_data_ptr p = malloc(sizeof(*p));
+  p->x = x;
+  p->y = y;
+  for (int i = 0; i < 8; i++) mpz_init(p->a[i]);
+  mpz_set_si(p->a[1], 1);
+  mpz_set_si(p->a[6], 1);
+
+  return cf_new(bihom, p);
+}
index 85dc545..36fce20 100644 (file)
@@ -12,15 +12,8 @@ int main() {
   // Check e + pi
   cf_t e = cf_new_e();
   cf_t pi = cf_new_pi();
-  mpz_t addarray[8];
-  for (int i = 0; i < 8; i++) {
-    mpz_init(addarray[i]);
-  }
-  mpz_set_ui(addarray[1], 1);
-  mpz_set_ui(addarray[2], 1);
-  mpz_set_ui(addarray[7], 1);
 
-  cf_t b = cf_new_bihom(e, pi, addarray);
+  cf_t b = cf_new_add(e, pi);
 
   CF_EXPECT_DEC(b, "5.8598744820488384738");
 
@@ -28,15 +21,28 @@ int main() {
   cf_free(e);
   cf_free(pi);
 
-  // Check 2 sin 1 cos 2 = sin 2 = 0.909...
+  // Check e * pi
+  e = cf_new_e();
+  pi = cf_new_pi();
+
+  b = cf_new_mul(e, pi);
+
+  CF_EXPECT_DEC(b, "8.53973422267356706546");
+
+  cf_free(b);
+  cf_free(e);
+  cf_free(pi);
+
+  // Check 2 sin 1 cos 1 = sin 2.
   cf_t s1, c1;
   s1 = cf_new_sin1();
   c1 = cf_new_cos1();
-  mpz_set_ui(addarray[0], 2);
-  mpz_set_ui(addarray[1], 0);
-  mpz_set_ui(addarray[2], 0);
-  mpz_set_si(addarray[7], 1);
-  b = cf_new_bihom(s1, c1, addarray);
+  mpz_t a[8];
+  mpz8_init(a);
+  mpz8_set_int(a,
+      2, 0, 0, 0,
+      0, 0, 0, 1);
+  b = cf_new_bihom(s1, c1, a);
 
   CF_EXPECT_DEC(b, "0.9092974268256816953");
 
@@ -44,14 +50,14 @@ int main() {
   cf_free(c1);
   cf_free(s1);
 
-  // Check 2 (cos 1)^2 - 1 = cos 2 =
+  // Check 2 (cos 1)^2 - 1 = cos 2.
   c1 = cf_new_cos1();
   cf_t t[2];
   cf_tee(t, c1);
-  mpz_set_si(addarray[0], 2);
-  mpz_set_si(addarray[3], -1);
-  mpz_set_si(addarray[7], 1);
-  b = cf_new_bihom(t[0], t[1], addarray);
+  mpz8_set_int(a,
+      2, 0, 0, -1,
+      0, 0, 0, 1);
+  b = cf_new_bihom(t[0], t[1], a);
 
   CF_EXPECT_DEC(b, "-0.41614683654714238699");
 
@@ -60,8 +66,6 @@ int main() {
   cf_free(b);
   cf_free(c1);
 
-  for (int i = 0; i < 8; i++) {
-    mpz_clear(addarray[i]);
-  }
+  mpz8_clear(a);
   return 0;
 }
diff --git a/cf.h b/cf.h
index 4c1474d..7226ff7 100644 (file)
--- a/cf.h
+++ b/cf.h
@@ -56,10 +56,14 @@ cf_t cf_new_nonregular_to_cf(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d);
 // Does both of the above at once. Seems slow.
 cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a[4]);
 
+cf_t cf_new_const_nonregular(void *(*fun)(cf_t));
+cf_t cf_new_one_arg_nonregular(void *(*fun)(cf_t), mpz_t z);
+
 // Well-known continued fraction expansions.
 // cf_famous.c:
 // e:
 cf_t cf_new_sqrt2();
+cf_t cf_new_sqrt5();
 cf_t cf_new_e();
 cf_t cf_new_pi();
 cf_t cf_new_tan1();
@@ -74,6 +78,20 @@ cf_t cf_new_tan(mpz_t z);
 
 // Gosper's method for computing bihomographic functions of continued fractions.
 cf_t cf_new_bihom(cf_t x, cf_t y, mpz_t a[8]);
+cf_t cf_new_add(cf_t x, cf_t y);
+cf_t cf_new_sub(cf_t x, cf_t y);
+cf_t cf_new_mul(cf_t x, cf_t y);
+cf_t cf_new_div(cf_t x, cf_t y);
+
+void mpz8_init(mpz_t z[8]);
+void mpz8_clear(mpz_t z[8]);
+void mpz8_set_int(mpz_t z[8],
+    int a, int b, int c, int d,
+    int e, int f, int g, int h);
+void mpz8_set_add(mpz_t z[8]);
+void mpz8_set_sub(mpz_t z[8]);
+void mpz8_set_mul(mpz_t z[8]);
+void mpz8_set_div(mpz_t z[8]);
 
 // From taylor.c:
 cf_t cf_new_sin1();
@@ -87,5 +105,7 @@ cf_t cf_new_cos1();
 //     a4 xy - a0 x + a5 y - a2
 cf_t cf_new_newton(cf_t x, mpz_t a[6], mpz_t lower);
 cf_t cf_new_sqrt(cf_t x);
+cf_t cf_new_sqrt_int(int a, int b);
+cf_t cf_new_sqrt_pq(mpz_t zp, mpz_t zq);
 
 #endif  // __CF_H__
index 48e279f..df33b5a 100644 (file)
--- a/famous.c
+++ b/famous.c
@@ -1,19 +1,27 @@
+// TODO: Use cf_new_const_nonregular instead of regularized_pi by allowing
+// arbitrary starting Mobius function.
 #include <stdio.h>
 #include <stdlib.h>
 #include <gmp.h>
 #include "cf.h"
 
-// sqrt(2) = [1; 2, 2, ...]
-static void *sqrt2(cf_t cf) {
-  cf_put_int(cf, 1);
+// sqrt(n^2 + 1) = [n; 2n, 2n, ...]
+static void *sqrt_easy(cf_t cf) {
+  unsigned int n = (unsigned int) cf_data(cf);
+  cf_put_int(cf, n);
+  n += n;
   while(cf_wait(cf)) {
-    cf_put_int(cf, 2);
+    cf_put_int(cf, n);
   }
   return NULL;
 }
 
 cf_t cf_new_sqrt2() {
-  return cf_new_const(sqrt2);
+  return cf_new(sqrt_easy, (void *) 1);
+}
+
+cf_t cf_new_sqrt5() {
+  return cf_new(sqrt_easy, (void *) 2);
 }
 
 // e = [2; 1, 2, 1, 1, 4, 1, ...]
@@ -76,10 +84,6 @@ static void *regularized_pi(cf_t cf) {
   return NULL;
 }
 
-cf_t cf_new_pi() {
-  return cf_new_const(regularized_pi);
-}
-
 // tan 1 = [1; 1, 1, 3, 1, 5, ...] 
 static void *tan1_expansion(cf_t cf) {
   int odd = 1;
@@ -163,62 +167,18 @@ static void *gauss_tan_expansion(cf_t cf) {
   return NULL;
 }
 
-cf_t cf_new_one_arg(void *(*fun)(cf_t), mpz_t z) {
-  mpz_ptr p = malloc(sizeof(*p));
-  mpz_init(p);
-  mpz_set(p, z);
-  return cf_new(fun, p);
-}
-
-struct funarg_s {
-  void *(*fun)(cf_t);
-  mpz_t arg;
-};
-typedef struct funarg_s *funarg_ptr;
-
-static void *one_arg_nonreg(cf_t cf) {
-  funarg_ptr p = cf_data(cf);
-  mpz_t a, b, c, d;
-  mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
-  mpz_set_ui(a, 1); mpz_set_ui(b, 0);
-  mpz_set_ui(c, 0); mpz_set_ui(d, 1);
-  mpz_ptr copy = malloc(sizeof(*copy));
-  mpz_init(copy);
-  mpz_set(copy, p->arg);
-  cf_t nonreg = cf_new(p->fun, copy);
-  cf_t conv = cf_new_nonregular_to_cf(nonreg, a, b, c, d);
-  mpz_t z;
-  mpz_init(z);
-  while(cf_wait(cf)) {
-    cf_get(z, conv);
-    cf_put(cf, z);
-  }
-  mpz_clear(z);
-  cf_free(conv);
-  cf_free(nonreg);
-  mpz_clear(a); mpz_clear(b); mpz_clear(c); mpz_clear(d);
-
-  mpz_clear(p->arg);
-  free(p);
-  return NULL;
-}
-
-cf_t cf_new_one_arg_nonreg(void *(*fun)(cf_t), mpz_t z) {
-  funarg_ptr p = malloc(sizeof(*p));
-  p->fun = fun;
-  mpz_init(p->arg);
-  mpz_set(p->arg, z);
-  return cf_new(one_arg_nonreg, p);
-}
-
 cf_t cf_new_epow(mpz_t pow) {
-  return cf_new_one_arg_nonreg(exp_expansion, pow);
+  return cf_new_one_arg_nonregular(exp_expansion, pow);
 }
 
 cf_t cf_new_tanh(mpz_t z) {
-  return cf_new_one_arg_nonreg(gauss_tanh_expansion, z);
+  return cf_new_one_arg_nonregular(gauss_tanh_expansion, z);
 }
 
 cf_t cf_new_tan(mpz_t z) {
-  return cf_new_one_arg_nonreg(gauss_tan_expansion, z);
+  return cf_new_one_arg_nonregular(gauss_tan_expansion, z);
+}
+
+cf_t cf_new_pi() {
+  return cf_new_const(regularized_pi);
 }
diff --git a/hakmem.c b/hakmem.c
new file mode 100644 (file)
index 0000000..3dc6d97
--- /dev/null
+++ b/hakmem.c
@@ -0,0 +1,125 @@
+// Compute the HAKMEM constant.
+//   (sqrt(3/pi^2 + e)/(tanh(sqrt(5))-sin(69))
+// which is
+//  $ echo "pi=4*a(1)
+//          x=2*sqrt(5)
+//          (sqrt(3/pi^2+e(1)))/((e(x)-1)/(e(x)+1)-s(69))" | bc -l
+// Almost exclusively uses techniques described by Gosper.
+// Differences:
+//  - Taylor series for cos(1) to find its continued fraction expansion.
+//  - Binary search instead of Newton's method when finding integer part of
+//  solution of quadratic.
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "cf.h"
+
+void *tanhsqrt5denom(cf_t cf) {
+  int odd = 1;
+  while(cf_wait(cf)) {
+    cf_put_int(cf, odd);
+    odd += 2;
+    cf_put_int(cf, 5);
+  }
+  return NULL;
+}
+
+void cf_dump(cf_t cf, int n) {
+  mpz_t z;
+  mpz_init(z);
+  cf_t dec = cf_new_cf_to_decimal(cf);
+  for (int i = 0; i < n; i++) {
+    cf_get(z, dec);
+    gmp_printf("%Zd", z);
+    if (!(i % 5)) putchar(' ');
+    if (!(i % 50)) putchar('\n');
+  }
+  printf("\n");
+  cf_free(dec);
+  mpz_clear(z);
+}
+
+int main() {
+  int n = 500;
+  cf_t c[7];
+  cf_t t[6][2];
+  mpz_t b[8];
+  int i;
+
+  mpz8_init(b);
+  mpz8_set_int(b,
+      2, 0, 0, -1,
+      0, 0, 0, 1);
+
+  c[0] = cf_new_cos1();
+  for (i = 1; i < 7; i++) {
+    cf_tee(t[i - 1], c[i - 1]);
+    c[i] = cf_new_bihom(t[i - 1][0], t[i - 1][1], b);
+  }
+  // c[6] = cos 64
+
+  cf_t c1 = cf_new_cos1();
+  cf_t ct0[2], ct1[2], ct2[2];
+  cf_tee(ct1, c1);
+  cf_tee(ct2, ct1[1]);
+  cf_t c2 = cf_new_mul(ct2[0], ct2[1]);
+  cf_tee(ct0, c2);
+  mpz8_set_int(b,
+      16, -20, 0, 5,  // cos 5n (5th Chebyshev polynomial)
+      0, 0, 0, 1);
+  cf_t c3 = cf_new_bihom(ct0[0], ct0[1], b);
+  cf_t c5 = cf_new_mul(ct1[0], c3);
+  // c5 = cos 5
+
+  cf_t c64t[2], c64t1[2];
+  cf_tee(c64t, c[6]);
+  cf_tee(c64t1, c64t[1]);
+  mpz8_set_int(b,
+      -1, 0, 0, 1,
+       0, 0, 0, 1);
+  cf_t s1 = cf_new_bihom(c64t1[0], c64t1[1], b);
+  cf_t s64 = cf_new_sqrt(s1);
+  // s64 = sin 64, c64t[0] = untouched cos 64
+
+  cf_t c5t[2], c5t1[2];
+  cf_tee(c5t, c5);
+  cf_tee(c5t1, c5t[1]);
+  cf_t s3 = cf_new_bihom(c5t1[0], c5t1[1], b);
+  cf_t s5 = cf_new_sqrt(s3);
+  // s5 = sin 5, c5t[0] = untouched cos 5
+
+  cf_t cf0 = cf_new_mul(s5, c64t[0]);
+  cf_t cf1 = cf_new_mul(s64, c5t[0]);
+  cf_t s69 = cf_new_sub(cf0, cf1); // TODO: Respect signs.
+                                  // This is supposed to be an addition,
+                                  // and the result should be negative.
+  // s69 = sin 69
+
+  cf_t sqrt5 = cf_new_sqrt5();
+  cf_t ts5d = cf_new_const_nonregular(tanhsqrt5denom);
+  cf_t ts5 = cf_new_div(sqrt5, ts5d);
+  
+  cf_t den = cf_new_add(ts5, s69); // TODO: This should be a subtract, but
+                                   // again, we don't handle signs properly.
+
+  cf_t e = cf_new_e();
+  cf_t pi = cf_new_pi();
+  cf_t pitee[2];
+  cf_tee(pitee, pi);
+  mpz8_set_int(b,
+      0, 0, 0, 3,
+      1, 0, 0, 0);
+  // tops = Three-Over-Pi-Squared.
+  cf_t tops = cf_new_bihom(pitee[0], pitee[1], b);
+
+  mpz8_set_add(b);
+  cf_t sum = cf_new_bihom(tops, e, b);
+  cf_t num = cf_new_sqrt(sum);
+  cf_t hakmem_constant = cf_new_div(num, den);
+  cf_dump(hakmem_constant, n);
+
+  mpz8_clear(b);
+  return 0;
+}
index 0382587..ee62187 100644 (file)
--- a/mobius.c
+++ b/mobius.c
@@ -454,3 +454,79 @@ cf_t cf_new_cf_to_decimal(cf_t x) {
   mpz_clear(one); mpz_clear(zero);
   return res;
 }
+
+struct funarg_s {
+  void *(*fun)(cf_t);
+  mpz_t arg;
+};
+typedef struct funarg_s *funarg_ptr;
+
+static void *one_arg_nonregular(cf_t cf) {
+  funarg_ptr p = cf_data(cf);
+  mpz_t a, b, c, d;
+  mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
+  mpz_set_ui(a, 1); mpz_set_ui(b, 0);
+  mpz_set_ui(c, 0); mpz_set_ui(d, 1);
+  mpz_ptr copy = malloc(sizeof(*copy));
+  mpz_init(copy);
+  mpz_set(copy, p->arg);
+  cf_t nonregular = cf_new(p->fun, copy);
+  cf_t conv = cf_new_nonregular_to_cf(nonregular, a, b, c, d);
+  mpz_t z;
+  mpz_init(z);
+  while(cf_wait(cf)) {
+    cf_get(z, conv);
+    cf_put(cf, z);
+  }
+  mpz_clear(z);
+  cf_free(conv);
+  cf_free(nonregular);
+  mpz_clear(a); mpz_clear(b); mpz_clear(c); mpz_clear(d);
+
+  mpz_clear(p->arg);
+  free(p);
+  return NULL;
+}
+
+cf_t cf_new_one_arg(void *(*fun)(cf_t), mpz_t z) {
+  mpz_ptr p = malloc(sizeof(*p));
+  mpz_init(p);
+  mpz_set(p, z);
+  return cf_new(fun, p);
+}
+
+cf_t cf_new_one_arg_nonregular(void *(*fun)(cf_t), mpz_t z) {
+  funarg_ptr p = malloc(sizeof(*p));
+  p->fun = fun;
+  mpz_init(p->arg);
+  mpz_set(p->arg, z);
+  return cf_new(one_arg_nonregular, p);
+}
+
+static void *nonreg_const(cf_t cf) {
+  funarg_ptr p = cf_data(cf);
+  mpz_t a, b, c, d;
+  mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
+  mpz_set_ui(a, 1); mpz_set_ui(b, 0);
+  mpz_set_ui(c, 0); mpz_set_ui(d, 1);
+  cf_t nonregular = cf_new(p->fun, NULL);
+  cf_t conv = cf_new_nonregular_to_cf(nonregular, a, b, c, d);
+  mpz_clear(a); mpz_clear(b); mpz_clear(c); mpz_clear(d);
+  mpz_t z;
+  mpz_init(z);
+  while(cf_wait(cf)) {
+    cf_get(z, conv);
+    cf_put(cf, z);
+  }
+  mpz_clear(z);
+  cf_free(conv);
+  cf_free(nonregular);
+  free(p);
+  return NULL;
+}
+
+cf_t cf_new_const_nonregular(void *(*fun)(cf_t)) {
+  funarg_ptr p = malloc(sizeof(*p));
+  p->fun = fun;
+  return cf_new(nonreg_const, p);
+}
index aeb8ce8..6b9775d 100644 (file)
--- a/newton.c
+++ b/newton.c
@@ -1,18 +1,20 @@
 // Solve quadratic equations involving continued fractions.
 // A slight misnomer: Gosper describes how to use Newton's method, but
 // I use binary search instead, and assume unique roots in the right ranges.
+//
+// For technical reasons (see Gosper), we write quadratics as
+//
+//     a0 xy + a1 x + a2 y + a3
+// y = ------------------------
+//     a4 xy - a0 x + a5 y - a2
+//
+// where y is the variable and x is some fixed continued fraction constant.
 
 #include <stdio.h>
 #include <stdlib.h>
 #include <gmp.h>
 #include "cf.h"
 
-// Represents:
-//
-// a0 xy + a1 x + a2 y + a3
-// ------------------------
-// a4 xy - a0 x + a5 y - a2
-
 struct newton_data_s {
   cf_t x;
   mpz_t a[6];
@@ -184,7 +186,9 @@ static void *newton(cf_t cf) {
   abc_clear(p);
   mpz_clear(z); mpz_clear(z0); mpz_clear(z1); mpz_clear(pow2);
   mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
+  mpz_clear(one);
   for (int i = 0; i < 6; i++) mpz_clear(nd->a[i]);
+  mpz_clear(nd->lower);
   free(nd);
   return NULL;
 }
@@ -211,3 +215,117 @@ cf_t cf_new_sqrt(cf_t x) {
   mpz_set_si(p->a[5], 1);
   return cf_new(newton, p);
 }
+
+struct newton_int_data_s {
+  mpz_t coeff[3];
+  mpz_t lower;
+};
+typedef struct newton_int_data_s *newton_int_data_ptr;
+typedef struct newton_int_data_s newton_int_data_t[1];
+
+// Finds smallest root greater than given lower bound of quadratic equation
+// with integer coefficients. If the coefficient of x is odd, we double
+// all the coefficients in our representation so we can find integers a, b, c
+// such that the quadratic may be written: c y^2 - 2 a y - b = 0.
+//
+// Assumes there exists exactly one root above given bound.
+// (Impossible to solve equations if roots have same integer part at
+// the moment. Don't do that.)
+// TODO: Rewrite so you choose if you want smaller or greater root.
+// or so it returns both solutions in an array.
+static void *newton_integer(cf_t cf) {
+  newton_int_data_ptr p = cf_data(cf);
+  mpz_t a, b, c;
+  mpz_init(a); mpz_init(b); mpz_init(c);
+  mpz_set(a, p->coeff[0]);
+  mpz_set(b, p->coeff[1]);
+  mpz_set(c, p->coeff[2]);
+
+  mpz_t z, z0, z1, one, pow2, t0, t1, t2;
+  mpz_init(z); mpz_init(z0); mpz_init(z1); mpz_init(pow2);
+  mpz_init(t0); mpz_init(t1); mpz_init(t2);
+  mpz_init(one);
+  mpz_set_ui(one, 1);
+
+  void move_down() {
+    // Recurrence relation with z, then
+    // Subtract and invert z.
+    mpz_mul(t0, z, a);
+    mpz_add(t0, t0, b);
+    mpz_mul(t1, z, c);
+    mpz_sub(t1, t1, a);
+    mpz_mul(a, z, t1);
+    mpz_set(b, c);
+    mpz_sub(c, t0, a);
+    mpz_set(a, t1);
+  }
+
+  int sign_quad() {
+    // Returns sign of c z^2 - 2 a z - b
+    mpz_mul_si(t0, a, -2);
+    mpz_addmul(t0, c, z);
+    mpz_mul(t0, t0, z);
+    mpz_sub(t0, t0, b);
+    return mpz_sgn(t0);
+  }
+
+  // Get integer part, starting search from given lower bound.
+  void binary_search(mpz_ptr lower) {
+    mpz_set(z0, lower);
+    mpz_set(z, lower);
+    int sign = sign_quad();
+    mpz_set_ui(pow2, 1);
+    for (;;) {
+      mpz_add(z, z0, pow2);
+      if (sign_quad() != sign) break;
+      mpz_mul_2exp(pow2, pow2, 1);
+    }
+    mpz_set(z1, z);
+
+    for (;;) {
+      mpz_add(z, z0, z1);
+      mpz_div_2exp(z, z, 1);
+      if (!mpz_cmp(z, z0)) break;
+      if (sign_quad() == sign) {
+       mpz_set(z0, z);
+      } else {
+       mpz_set(z1, z);
+      }
+    }
+    cf_put(cf, z0);
+    mpz_set(z, z0);
+    move_down();
+  }
+
+  binary_search(p->lower);
+
+  while(cf_wait(cf)) {
+    binary_search(one);
+  }
+  mpz_clear(z); mpz_clear(z0); mpz_clear(z1); mpz_clear(pow2);
+  mpz_clear(t0); mpz_clear(t1); mpz_clear(t2);
+  mpz_clear(one);
+  for (int i = 0; i < 3; i++) mpz_clear(p->coeff[i]);
+  mpz_clear(p->lower);
+  free(p);
+  return NULL;
+}
+
+cf_t cf_new_sqrt_pq(mpz_t zp, mpz_t zq) {
+  newton_int_data_ptr p = malloc(sizeof(*p));
+  mpz_init(p->lower);
+  for (int i = 0; i < 3; i++) mpz_init(p->coeff[i]);
+  mpz_set(p->coeff[1], zp);
+  mpz_set(p->coeff[2], zq);
+  return cf_new(newton_integer, p);
+}
+
+cf_t cf_new_sqrt_int(int a, int b) {
+  mpz_t p, q;
+  mpz_init(p); mpz_init(q);
+  mpz_set_si(p, a);
+  mpz_set_si(q, b);
+  cf_t res = cf_new_sqrt_pq(p, q);
+  mpz_clear(p); mpz_clear(q);
+  return res;
+}
index 52ce9d5..2c2bc7c 100644 (file)
@@ -37,10 +37,10 @@ int main() {
   cf_t s1[2];
   cf_tee(s1, x);
   mpz_t b[8];
-  for (i = 0; i < 8; i++) mpz_init(b[i]);
-  mpz_set_si(b[0], -1);
-  mpz_set_si(b[3], 1);
-  mpz_set_si(b[7], 1);
+  mpz8_init(b);
+  mpz8_set_int(b,
+      -1, 0, 0, 1,
+       0, 0, 0, 1);
 
   cf_t bi = cf_new_bihom(s1[0], s1[1], b);
   cf_t n = cf_new_sqrt(bi);
@@ -52,7 +52,11 @@ int main() {
   cf_free(bi);
   cf_free(n);
 
-  for (i = 0; i < 8; i++) mpz_clear(b[i]);
+  x = cf_new_sqrt_int(355, 113);
+  CF_EXPECT_DEC(x, "1.77245392615830279609");
+  cf_free(x);
+
+  mpz8_clear(b);
   for (i = 0; i < 6; i++) mpz_clear(a[i]);
   return 0;
 }
diff --git a/nonregular.c b/nonregular.c
new file mode 100644 (file)
index 0000000..badbbf0
--- /dev/null
@@ -0,0 +1,4 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <gmp.h>
+#include "cf.h"