From: Ben Lynn Date: Wed, 17 Sep 2008 09:04:23 +0000 (-0700) Subject: Replaced Newton's method with binary search. X-Git-Url: https://repo.or.cz/w/frac.git/commitdiff_plain/89898d81e965f3ee74a2b3b5f900b51d87c7cd1e Replaced Newton's method with binary search. --- diff --git a/Makefile b/Makefile index c61c560..3420387 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ CF_OBJS:=cf.o mobius.o famous.o bihom.o taylor.o test.o newton.o TESTS:=bihom_test cf_test famous_test mobius_test newton_test -BINS:=pi pi2 sqrt2 epow +BINS:=pi epow target : $(BINS) diff --git a/README b/README index 99afa0d..9711093 100644 --- a/README +++ b/README @@ -126,15 +126,19 @@ in this code. 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. +=== Disclaimer === + +Let me put my salesmanship on hold, and temper 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. + +Also, Newton's method, and other floating-point algorithms converge +quadratically or faster, thus a good implementation will asymptotically +outperform continued fractions as these converge only linearly. 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: diff --git a/cf.h b/cf.h index f250a9c..aab54a4 100644 --- a/cf.h +++ b/cf.h @@ -47,7 +47,7 @@ cf_t cf_new_nonregular_mobius_convergent(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_ // Output: Regular continued fraction. Assumes input fraction is well-behaved. 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, mpz_t b, mpz_t c, mpz_t d); +cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a[4]); // Well-known continued fraction expansions. // cf_famous.c: @@ -78,6 +78,7 @@ cf_t cf_new_cos1(); // a0 xy + a1 x + a2 y + a3 // y = ------------------------ // a4 xy - a0 x + a5 y - a2 -cf_t cf_new_newton(cf_t x, mpz_t a[6]); +cf_t cf_new_newton(cf_t x, mpz_t a[6], mpz_t lower); +cf_t cf_new_sqrt(cf_t x); #endif // __CF_H__ diff --git a/famous_test.c b/famous_test.c index 7347660..6f9b951 100644 --- a/famous_test.c +++ b/famous_test.c @@ -15,5 +15,21 @@ int main() { CF_NEW_EXPECT_DEC(cf_new_tan1, "1.5574077246549022305"); CF_NEW_EXPECT_DEC(cf_new_sin1, "0.8414709848078965066"); CF_NEW_EXPECT_DEC(cf_new_cos1, "0.5403023058681397174"); + + mpz_t z; + mpz_init(z); + cf_t x; + + mpz_set_si(z, 3); + x = cf_new_tanh(z); + CF_EXPECT_DEC(x, "0.99505475368673045133188018525548847509781385470028"); + cf_free(x); + + mpz_set_ui(z, 2); + x = cf_new_epow(z); + CF_EXPECT_DEC(x, "7.38905609893065022723042746057500781318031557055184"); + cf_free(x); + + mpz_clear(z); return 0; } diff --git a/mobius.c b/mobius.c index 0d00f12..0382587 100644 --- a/mobius.c +++ b/mobius.c @@ -257,7 +257,7 @@ static void *nonregular_mobius_decimal(cf_t cf) { mpz_add(t2, t2, pq->q); if (mpz_cmp(t2, pq->p) > 0) { // Output a decimal digit. - cf_put(cf, pq->p); + cf_put(cf, t1); // Subtract: remainder of p/q. mpz_sub(t2, t2, pq->p); mpz_sub(t2, pq->q, t2); @@ -288,10 +288,11 @@ static void *nonregular_mobius_decimal(cf_t cf) { return NULL; } -cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a, mpz_t b, mpz_t c, mpz_t d) { +cf_t cf_new_nonregular_mobius_to_decimal(cf_t x, mpz_t a[4]) { mobius_data_ptr md = malloc(sizeof(*md)); mpz_init(md->a); mpz_init(md->b); mpz_init(md->c); mpz_init(md->d); - mpz_set(md->a, a); mpz_set(md->b, b); mpz_set(md->c, c); mpz_set(md->d, d); + mpz_set(md->a, a[0]); mpz_set(md->b, a[1]); + mpz_set(md->c, a[2]); mpz_set(md->d, a[3]); md->input = x; return cf_new(nonregular_mobius_decimal, md); } diff --git a/mobius_test.c b/mobius_test.c index 150d840..41f46bc 100644 --- a/mobius_test.c +++ b/mobius_test.c @@ -12,6 +12,34 @@ static void *sqrt2(cf_t cf) { return NULL; } +// Converges extremely slowly. +static void *slow_pi(cf_t cf) { + mpz_t num, denom, t; + mpz_init(num); + mpz_init(denom); + mpz_init(t); + + mpz_set_ui(denom, 3); + cf_put(cf, denom); + + mpz_set_ui(num, 1); + cf_put(cf, num); + mpz_set_ui(t, 8); + + mpz_set_ui(denom, 6); + while(cf_wait(cf)) { + cf_put(cf, denom); + mpz_add(num, num, t); + cf_put(cf, num); + mpz_add_ui(t, t, 8); + } + + mpz_clear(num); + mpz_clear(denom); + mpz_clear(t); + return NULL; +} + int main() { cf_t x, conv; x = cf_new_const(sqrt2); @@ -46,13 +74,37 @@ int main() { cf_free(conv); cf_free(x); - x = cf_new_const(sqrt2); - cf_t mob; mpz_t z[4]; for (int i = 0; i < 4; i++) mpz_init(z[i]); // Identity Mobius transformation. mpz_set_si(z[0], 1); mpz_set_si(z[3], 1); + mpz_t digit; + mpz_init(digit); + + x = cf_new_const(slow_pi); + conv = cf_new_nonregular_mobius_to_decimal(x, z); + cf_get(digit, conv); + EXPECT(!mpz_cmp_ui(digit, 3)); + cf_get(digit, conv); + EXPECT(!mpz_cmp_ui(digit, 1)); + cf_get(digit, conv); + EXPECT(!mpz_cmp_ui(digit, 4)); + cf_get(digit, conv); + EXPECT(!mpz_cmp_ui(digit, 1)); + cf_get(digit, conv); + EXPECT(!mpz_cmp_ui(digit, 5)); + cf_get(digit, conv); + EXPECT(!mpz_cmp_ui(digit, 9)); + cf_get(digit, conv); + EXPECT(!mpz_cmp_ui(digit, 2)); + + mpz_clear(digit); + cf_free(x); + cf_free(conv); + + x = cf_new_const(sqrt2); + cf_t mob; mob = cf_new_mobius_to_cf(x, z); CF_EXPECT_DEC(mob, "1.4142135623730"); cf_free(mob); diff --git a/pi2.c b/pi2.c deleted file mode 100644 index 501f0eb..0000000 --- a/pi2.c +++ /dev/null @@ -1,56 +0,0 @@ -#include -#include -#include -#include "cf.h" - -// Converges extremely slowly. Not worth converting to decimal or -// continued fraction form. -static void *pi_six_sequence(cf_t cf) { - mpz_t num, denom, t; - mpz_init(num); - mpz_init(denom); - mpz_init(t); - - mpz_set_ui(denom, 3); - cf_put(cf, denom); - - mpz_set_ui(num, 1); - cf_put(cf, num); - mpz_set_ui(t, 8); - - mpz_set_ui(denom, 6); - while(cf_wait(cf)) { - cf_put(cf, denom); - mpz_add(num, num, t); - cf_put(cf, num); - mpz_add_ui(t, t, 8); - } - - mpz_clear(num); - mpz_clear(denom); - mpz_clear(t); - return NULL; -} - -int main() { - mpz_t z; - mpz_t a, b, c, d; - mpz_init(z); - cf_t pi, conv; - pi = cf_new(pi_six_sequence, NULL); - 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); - - conv = cf_new_nonregular_mobius_to_decimal(pi, a, b, c, d); - for (int i = 1; i <= 5000; i++) { - cf_get(z, conv); - gmp_printf("%Zd\n",z); - } - printf("\n"); - cf_free(conv); - cf_free(pi); - mpz_clear(z); - mpz_clear(d); mpz_clear(c); mpz_clear(b); mpz_clear(a); - return 0; -}