Added sign of a continued fraction.
authorBen Lynn <benlynn@gmail.com>
Tue, 16 Sep 2008 09:31:59 +0000 (16 02:31 -0700)
committerBen Lynn <benlynn@gmail.com>
Tue, 16 Sep 2008 09:45:18 +0000 (16 02:45 -0700)
Determine sign in mobius_decimal().

Added exercise.c: computes 37 digits of sqrt(2) without multiplication
in plain C.

README
bihom.c
bihom_test.c [new file with mode: 0644]
cf.c
cf.h
exercise.c [new file with mode: 0644]
mobius.c
pi.c
sqrt2.c

diff --git a/README b/README
index e0db81d..99afa0d 100644 (file)
--- a/README
+++ b/README
@@ -1,7 +1,9 @@
 = Continued Fraction Library =
 Ben Lynn, September 2008
 
-== Introduction ==
+Frac is a C library for computing with continued fractions.
+
+== Continued Fractions ==
 
 1/3 has no finite floating point representation.  How can a standard tool for
 working with real numbers be unable to handle a number a young child could
@@ -30,8 +32,6 @@ If only we had used continued fractions, for then we could arbitrarily
 increase the precision of our answer without redoing any work. Furthermore,
 the error analysis becomes trivial.
 
-== Disadvantages ==
-
 We must pay a price for the advantages of continued fractions.  We
 have a representation from which we can instantly infer error bounds. A
 representation for which we may arbitrarily increase precision without redoing
@@ -41,11 +41,13 @@ whose output must conform to this wondrous representation, are clumsy and
 convoluted, and that deriving continued fraction expansions for many
 reals is time-consuming?
 
-These drawbacks are almost certainly why continued fractions remain obscure.
-Yet surely there are some problems that scream for continued fractions?
-Imagine a real-time application where on a timer interrupt, all continued
-fractions simply output the last computed convergent.  Under heavy system load,
-approximations are coarser but the show goes on.
+These drawbacks are almost certainly why these fascinating creatures remain
+obscure.  Yet surely there are some problems that scream for continued
+fractions?  How about a screen saver that computes more and more digits of pi,
+each time picking up where it left off?  Or imagine a real-time application
+where all continued fractions simply output the last computed convergent on a
+timer interrupt.  Under heavy system load, approximations are coarser but the
+show goes on.
 
 == Installation ==
 
@@ -123,3 +125,17 @@ in this code.
   fractions a "tee"-like function is more apt: rather than spawning a thread
   for each term held back, we spawn and watch over two threads on invocation to
   service the two demand channels, storing all backlogged terms in one place.
+
+=== Increasing floating-point precision ===
+
+A disclaimer to some of my anti-floating-point rhetoric: increasing the
+precision of floating-point operations without redoing work is in fact possible
+for many common cases. For example, Newton's method is self-correcting, meaning
+we may use low precision at first, and increase it for future iterations.
+Even pi enjoys such methods: there is a formula computing any hexadecimal
+digit of pi without computing any previous digits, though it requires about the
+same amount of work.
+
+However, I'd prefer continued fractions to Newton's method when writing code
+to compute digits of, say, sqrt(2). Precision decisions are made automatically:
+one simply needs enough bits to hold the last computed convergents.
diff --git a/bihom.c b/bihom.c
index 4180774..652040d 100644 (file)
--- a/bihom.c
+++ b/bihom.c
@@ -90,6 +90,7 @@ static void *bihom(cf_t cf) {
     mpz_set(p->r0, p->p0);  mpz_set(p->r1, p->p1);
     mpz_set(p->p0, t0);     mpz_set(p->p1, t1);
   }
+
   int recur() {
     if (!mpz_sgn(p->s1)) {
       move_right();
diff --git a/bihom_test.c b/bihom_test.c
new file mode 100644 (file)
index 0000000..874beee
--- /dev/null
@@ -0,0 +1,50 @@
+// Test quadratic algorithm.
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gmp.h>
+#include "cf.h"
+#include "test.h"
+
+int main() {
+  // I used bc to get the reference values.
+  // 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_EXPECT_DEC(b, "58598744820488384738");
+
+  cf_free(b);
+  cf_free(e);
+  cf_free(pi);
+
+  // Check 2 sin 1 cos 2 = sin 2 = 0.909...
+  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);
+
+  CF_EXPECT_DEC(b, "09092974268256816953");
+
+  cf_free(b);
+  cf_free(c1);
+  cf_free(s1);
+  for (int i = 0; i < 8; i++) {
+    mpz_clear(addarray[i]);
+  }
+  return 0;
+}
diff --git a/cf.c b/cf.c
index 51eca4f..2bad1b5 100644 (file)
--- a/cf.c
+++ b/cf.c
@@ -28,6 +28,12 @@ struct cf_s {
   pthread_mutex_t chan_mu;
   channel_ptr chan, next;
 
+  // We break the CSP model slightly here: the sign of the continued
+  // fraction is read directly from a variable, not over channels.
+  // Only 'thread' may write to 'sign', and it should do so before the first
+  // cf_put(). Other threads should only read it after they have
+  // called cf_get() at least once.
+  int sign;
   int quitflag;
   void *data;
 };
@@ -38,6 +44,14 @@ void *cf_data(cf_t cf) {
   return cf->data;
 }
 
+int cf_sign(cf_t cf) {
+  return cf->sign;
+}
+
+int cf_flip_sign(cf_t cf) {
+  return cf->sign = -cf->sign;
+}
+
 // A bit like cooperative multitasking. Continued fractions are expected
 // to call this as often as practical, and on a return value of 0,
 // to drop everything and stop.
@@ -127,6 +141,7 @@ void cf_get(mpz_t z, cf_t cf) {
 
 cf_t cf_new(void *(*func)(cf_t), void *data) {
   cf_t cf = malloc(sizeof(*cf));
+  cf->sign = 1;
   cf->chan = NULL;
   cf->next = NULL;
   cf->quitflag = 0;
diff --git a/cf.h b/cf.h
index 13b379c..1fd64e3 100644 (file)
--- a/cf.h
+++ b/cf.h
@@ -14,6 +14,8 @@ static inline cf_t cf_new_const(void *(*func)(cf_t)) {
 }
 void cf_free(cf_t cf);
 
+int cf_sign(cf_t cf);
+int cf_flip_sign(cf_t cf);
 void cf_get(mpz_t z, cf_t cf);
 void cf_put(cf_t cf, mpz_t z);
 void cf_put_int(cf_t cf, int n);
diff --git a/exercise.c b/exercise.c
new file mode 100644 (file)
index 0000000..0cc2542
--- /dev/null
@@ -0,0 +1,60 @@
+// Compute 37 digits of sqrt(2).
+
+#include <stdio.h>
+#include <stdlib.h>
+
+int main(int argc, char **argv) {
+  // 4 byte precision gives 18 digits,
+  // 8 byte precision gives 37 digits.
+  int n = 37;
+  if (argc > 1) {
+    n = atoi(argv[1]);
+    if (n <= 0) n = 37;
+  }
+  unsigned long long p0, p1, q0, q1, r0, r1;
+
+  // sqrt(2) = [1; 2, 2, ...]
+  // Handle first term here.
+  p0 = p1 = q1 = 1;
+  q0 = 0;
+
+  while (n > 0) {
+    // Continued fraction recurrence relations.
+    r0 = p0 + (p1 << 1);
+    p0 = p1;
+    p1 = r0;
+
+    r0 = q0 + (q1 << 1);
+    q0 = q1;
+    q1 = r0;
+
+    // digit <- floor(p0 / q0)
+    // r0 = q0 * (digit + 1)
+    // r1 = q1 * (digit + 1)
+    int digit = 0;
+    r0 = q0;
+    r1 = q1;
+    while (digit < 9 && r0 <= p0) {
+      r0 += q0;
+      r1 += q1;
+      digit++;
+    }
+
+    if (r1 > p1) {
+      r1 -= q1;
+      if (r1 <= p1) {
+       n--;
+       putchar(digit + '0');
+       p1 -= r1;
+       r0 -= q0;
+       p0 -= r0;
+
+       // p0 <- 10 * p0, p1 <- 10 * p1
+       p0 = ((p0 << 2) + p0) << 1;
+       p1 = ((p1 << 2) + p1) << 1;
+      }
+    }
+  }
+  putchar('\n');
+  return 0;
+}
index 68a2525..23b8373 100644 (file)
--- a/mobius.c
+++ b/mobius.c
@@ -240,13 +240,32 @@ static void *mobius_decimal(cf_t cf) {
   pqset_t pq; pqset_init(pq); pqset_set_mobius(pq, md);
   mpz_t denom; mpz_init(denom);
   mpz_t t0, t1, t2; mpz_init(t2); mpz_init(t1); mpz_init(t0);
+
+  // Determine the sign.
+  while (mpz_sgn(pq->pold) != mpz_sgn(pq->p)
+      || mpz_sgn(pq->qold) != mpz_sgn(pq->q)) {
+    cf_get(denom, input);
+    pqset_regular_recur(pq, denom);
+  }
+  if (mpz_sgn(pq->qold) < 0) {
+    mpz_neg(pq->qold, pq->qold);
+    mpz_neg(pq->q, pq->q);
+    cf_flip_sign(cf);
+  }
+  if (mpz_sgn(pq->pold) < 0) {
+    mpz_neg(pq->pold, pq->pold);
+    mpz_neg(pq->p, pq->p);
+    cf_flip_sign(cf);
+  }
+
   int recur() {
     pqset_regular_recur(pq, denom);
 
     // If the denominator is zero, we can't do anything yet.
     if (mpz_sgn(pq->qold)) {
-      // The answer is one of {0, ..., 9}.
-      /* Naive attempt to expoit this didn't work well:
+      // Each term except possibly the first is one of {0, ..., 9}.
+      /* Naive attempt to expoit this didn't seem faster:
+       * (and I'd have to handle the first term properly)
       int i;
       mpz_set(t0, pq->q);
       for (i = 0; i <= 9; i++) {
@@ -278,7 +297,6 @@ static void *mobius_decimal(cf_t cf) {
     }
     return 0;
   }
-
   while(cf_wait(cf)) {
     do {
       cf_get(denom, input);
diff --git a/pi.c b/pi.c
index bbb47ab..eaf467f 100644 (file)
--- a/pi.c
+++ b/pi.c
@@ -9,7 +9,8 @@ int main(int argc, char **argv) {
   mpz_init(z);
   cf_t pi, conv;
   pi = cf_new_pi();
-  int n = 1000;
+  int n = 2037 + 1;  // To outdo Metropolis, Reitwieser and von Neumann's
+                     // 1949 ENIAC record.
   if (argc > 1) {
     n = atoi(argv[1]);
     if (n <= 0) n = 100;
diff --git a/sqrt2.c b/sqrt2.c
index 6762460..6109ff4 100644 (file)
--- a/sqrt2.c
+++ b/sqrt2.c
@@ -30,6 +30,24 @@ int main() {
   }
   cf_free(conv);
   cf_free(x);
+
+  x = cf_new_sqrt2();
+  mpz_t a, b, c, d;
+  mpz_init(a); mpz_init(b); mpz_init(c); mpz_init(d);
+  mpz_set_si(a, -7);
+  mpz_set_si(b, 5);
+  mpz_set_si(c, -12);
+  mpz_set_si(d, 13);
+  cf_t dec = cf_new_mobius_to_decimal(x, a, b, c, d);
+  cf_get(p, dec);
+  putchar(cf_sign(dec) > 0 ? ' ' : '-');
+  gmp_printf("%Zd.", p);
+  for (int i = 0; i < 10; i++) {
+    cf_get(p, dec);
+    gmp_printf("%Zd", p);
+  }
+  putchar('\n');
+  cf_free(x);
   mpz_clear(p);
   mpz_clear(q);
   return 0;