From: Ben Lynn Date: Tue, 16 Sep 2008 09:31:59 +0000 (-0700) Subject: Added sign of a continued fraction. X-Git-Url: https://repo.or.cz/w/frac.git/commitdiff_plain/1ebbdb7d8be4322d8cbec566eb2f3565130c4d37 Added sign of a continued fraction. Determine sign in mobius_decimal(). Added exercise.c: computes 37 digits of sqrt(2) without multiplication in plain C. --- diff --git a/README b/README index e0db81d..99afa0d 100644 --- 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 --- 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 index 0000000..874beee --- /dev/null +++ b/bihom_test.c @@ -0,0 +1,50 @@ +// Test quadratic algorithm. + +#include +#include +#include +#include +#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 --- 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 --- 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 index 0000000..0cc2542 --- /dev/null +++ b/exercise.c @@ -0,0 +1,60 @@ +// Compute 37 digits of sqrt(2). + +#include +#include + +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; +} diff --git a/mobius.c b/mobius.c index 68a2525..23b8373 100644 --- 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 --- 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 --- 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;