From ee04081c68933664d6b200bf6216fa11373df907 Mon Sep 17 00:00:00 2001 From: Ben Lynn Date: Thu, 18 Sep 2008 03:20:48 -0700 Subject: [PATCH] Compute HAKMEM constant. --- Makefile | 9 ++-- README | 2 +- bihom.c | 97 +++++++++++++++++++++++++++++++++++++++++++ bihom_test.c | 48 ++++++++++++---------- cf.h | 20 +++++++++ famous.c | 80 +++++++++--------------------------- hakmem.c | 125 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ mobius.c | 76 ++++++++++++++++++++++++++++++++++ newton.c | 130 +++++++++++++++++++++++++++++++++++++++++++++++++++++++--- newton_test.c | 14 ++++--- nonregular.c | 4 ++ 11 files changed, 507 insertions(+), 98 deletions(-) create mode 100644 hakmem.c create mode 100644 nonregular.c diff --git a/Makefile b/Makefile index 522fc9a..d8ad84f 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 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 --- 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 --- 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); +} diff --git a/bihom_test.c b/bihom_test.c index 85dc545..36fce20 100644 --- a/bihom_test.c +++ b/bihom_test.c @@ -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 --- 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__ diff --git a/famous.c b/famous.c index 48e279f..df33b5a 100644 --- 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 #include #include #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 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 +#include +#include +#include +#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; +} diff --git a/mobius.c b/mobius.c index 0382587..ee62187 100644 --- 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); +} diff --git a/newton.c b/newton.c index aeb8ce8..6b9775d 100644 --- 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 #include #include #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; +} diff --git a/newton_test.c b/newton_test.c index 52ce9d5..2c2bc7c 100644 --- a/newton_test.c +++ b/newton_test.c @@ -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 index 0000000..badbbf0 --- /dev/null +++ b/nonregular.c @@ -0,0 +1,4 @@ +#include +#include +#include +#include "cf.h" -- 2.11.4.GIT