From 005aece5d7d02bd861c75f49368169be5d847147 Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 16:44:01 +0800 Subject: [PATCH 01/10] Add quadratic-ring binary splitting benchmark for D2 --- validation/benchmark_d2_bs_vs_classics.c | 320 +++++++++++++++++++++++ 1 file changed, 320 insertions(+) create mode 100644 validation/benchmark_d2_bs_vs_classics.c diff --git a/validation/benchmark_d2_bs_vs_classics.c b/validation/benchmark_d2_bs_vs_classics.c new file mode 100644 index 0000000..0b08718 --- /dev/null +++ b/validation/benchmark_d2_bs_vs_classics.c @@ -0,0 +1,320 @@ +#define main original_d2_benchmark_main +#include "benchmark_original_d2_basic_vs_chud.c" +#undef main + +/* + * Accelerated comparison only. + * + * D2 uses a balanced exact product tree in the quadratic ring generated by + * x = H2/J, where x^2 + H1*x + H2 = 0. This is the direct analogue of + * Chudnovsky binary splitting: the tree topology is unchanged, while scalar + * integer numerators are replaced by degree-2 algebraic numerators. + * + * The tree carries the common denominator Q, the ratio product P, the partial + * hypergeometric sum S, and its theta-weighted sum U = sum n*t_n. No modular + * polynomial, no (z,u,v) transport, and no repeated modular Newton solve is + * used. + */ + +typedef struct { + mpz_t a; + mpz_t b; +} d2_ring_t; + +typedef struct { + mpz_t q; + d2_ring_t p; + d2_ring_t s; + d2_ring_t u; +} d2_seg_t; + +static void d2_ring_init(d2_ring_t *r) { + mpz_init(r->a); + mpz_init(r->b); +} + +static void d2_ring_clear(d2_ring_t *r) { + mpz_clear(r->a); + mpz_clear(r->b); +} + +static void d2_seg_init(d2_seg_t *s) { + mpz_init(s->q); + d2_ring_init(&s->p); + d2_ring_init(&s->s); + d2_ring_init(&s->u); +} + +static void d2_seg_clear(d2_seg_t *s) { + mpz_clear(s->q); + d2_ring_clear(&s->p); + d2_ring_clear(&s->s); + d2_ring_clear(&s->u); +} + +/* + * (a+b*x)(c+d*x), x^2 = -H1*x-H2. + * Three large products are used for the bilinear part; products by H1/H2 + * have one fixed ~100-bit operand. + */ +static void d2_ring_mul(d2_ring_t *out, + const d2_ring_t *x, + const d2_ring_t *y, + const mpz_t h1, + const mpz_t h2, + mpz_t m0, + mpz_t m1, + mpz_t m2, + mpz_t t0, + mpz_t t1) { + mpz_mul(m0, x->a, y->a); + mpz_mul(m1, x->b, y->b); + + mpz_add(t0, x->a, x->b); + mpz_add(t1, y->a, y->b); + mpz_mul(m2, t0, t1); + + mpz_mul(t0, h2, m1); + mpz_sub(out->a, m0, t0); + + mpz_sub(m2, m2, m0); + mpz_sub(m2, m2, m1); + mpz_mul(t0, h1, m1); + mpz_sub(out->b, m2, t0); +} + +static void d2_leaf(d2_seg_t *out, + unsigned long n, + const mpz_t h2) { + mpz_t num, den, g; + mpz_inits(num, den, g, NULL); + + /* r_n = t_{n+1}/t_n = A_n*1728*x / (72*(n+1)^3*H2). */ + mpz_set_ui(num, 6UL * n + 1UL); + mpz_mul_ui(num, num, 2UL * n + 1UL); + mpz_mul_ui(num, num, 6UL * n + 5UL); + mpz_mul_ui(num, num, 1728UL); + + mpz_set_ui(den, n + 1UL); + mpz_mul_ui(den, den, n + 1UL); + mpz_mul_ui(den, den, n + 1UL); + mpz_mul_ui(den, den, 72UL); + mpz_mul(den, den, h2); + + /* Cancel the scalar part before it enters the tree. */ + mpz_gcd(g, num, den); + mpz_divexact(num, num, g); + mpz_divexact(out->q, den, g); + + mpz_set_ui(out->p.a, 0UL); + mpz_set(out->p.b, num); + + mpz_set(out->s.a, out->q); + mpz_set_ui(out->s.b, 0UL); + + mpz_mul_ui(out->u.a, out->q, n); + mpz_set_ui(out->u.b, 0UL); + + mpz_clears(num, den, g, NULL); +} + +static void d2_bs_build(d2_seg_t *out, + unsigned long a, + unsigned long b, + const mpz_t h1, + const mpz_t h2) { + if (b - a == 1UL) { + d2_leaf(out, a, h2); + return; + } + + unsigned long m = a + (b - a) / 2UL; + d2_seg_t left, right; + d2_seg_init(&left); + d2_seg_init(&right); + d2_bs_build(&left, a, m, h1, h2); + d2_bs_build(&right, m, b, h1, h2); + + mpz_t m0, m1, m2, t0, t1; + mpz_inits(m0, m1, m2, t0, t1, NULL); + d2_ring_t rs, ru; + d2_ring_init(&rs); + d2_ring_init(&ru); + + mpz_mul(out->q, left.q, right.q); + + d2_ring_mul(&out->p, &left.p, &right.p, + h1, h2, m0, m1, m2, t0, t1); + + d2_ring_mul(&rs, &left.p, &right.s, + h1, h2, m0, m1, m2, t0, t1); + mpz_mul(out->s.a, left.s.a, right.q); + mpz_add(out->s.a, out->s.a, rs.a); + mpz_mul(out->s.b, left.s.b, right.q); + mpz_add(out->s.b, out->s.b, rs.b); + + d2_ring_mul(&ru, &left.p, &right.u, + h1, h2, m0, m1, m2, t0, t1); + mpz_mul(out->u.a, left.u.a, right.q); + mpz_add(out->u.a, out->u.a, ru.a); + mpz_mul(out->u.b, left.u.b, right.q); + mpz_add(out->u.b, out->u.b, ru.b); + + d2_ring_clear(&rs); + d2_ring_clear(&ru); + mpz_clears(m0, m1, m2, t0, t1, NULL); + d2_seg_clear(&left); + d2_seg_clear(&right); +} + +static void d2_eval_ring(mpf_t out, + const d2_ring_t *r, + const mpz_t q, + const mpf_t x, + mp_bitcnt_t bits) { + mpf_t a, b, fq; + mpf_init2(a, bits); + mpf_init2(b, bits); + mpf_init2(fq, bits); + mpf_set_z(a, r->a); + mpf_set_z(b, r->b); + mpf_mul(b, b, x); + mpf_add(a, a, b); + mpf_set_z(fq, q); + mpf_div(out, a, fq); + mpf_clears(a, b, fq, NULL); +} + +static int d2_bs_pi(uint64_t digits, + uint64_t guard, + const char *path, + double *seconds, + unsigned long *terms_out, + mp_bitcnt_t extra_bits) { + double t0 = now_seconds(); + mp_bitcnt_t bits = rj_bits_for_decimal(digits, guard) + extra_bits; + + const long double digits_per_term = + 24.95589965765426673091470594098813130578L; + unsigned long terms = + (unsigned long)ceill(((long double)digits + (long double)guard + 30.0L) / + digits_per_term) + 2UL; + + mpz_t h1, h2, disc, tmpz; + mpz_inits(h1, h2, disc, tmpz, NULL); + if (mpz_set_str(h1, D2_H1, 10) != 0 || + mpz_set_str(h2, D2_H2, 10) != 0) { + mpz_clears(h1, h2, disc, tmpz, NULL); + return 0; + } + + d2_seg_t tree; + d2_seg_init(&tree); + d2_bs_build(&tree, 0UL, terms, h1, h2); + + /* Stable small embedding x = H2/J = -2 H2/(H1+sqrt(H1^2-4H2)). */ + mpz_mul(disc, h1, h1); + mpz_mul_ui(tmpz, h2, 4UL); + mpz_sub(disc, disc, tmpz); + + mpf_t fh1, fh2, fdisc, root, x, F, T, y, v, z, alpha, beta, + K, tmp, inv_pi, pi; + mpf_init2(fh1, bits); mpf_init2(fh2, bits); + mpf_init2(fdisc, bits); mpf_init2(root, bits); + mpf_init2(x, bits); mpf_init2(F, bits); + mpf_init2(T, bits); mpf_init2(y, bits); + mpf_init2(v, bits); mpf_init2(z, bits); + mpf_init2(alpha, bits); mpf_init2(beta, bits); + mpf_init2(K, bits); mpf_init2(tmp, bits); + mpf_init2(inv_pi, bits);mpf_init2(pi, bits); + + mpf_set_z(fh1, h1); + mpf_set_z(fh2, h2); + mpf_set_z(fdisc, disc); + mpf_sqrt(root, fdisc); + mpf_add(tmp, fh1, root); + mpf_mul_ui(x, fh2, 2UL); + mpf_neg(x, x); + mpf_div(x, x, tmp); + + d2_eval_ring(F, &tree.s, tree.q, x, bits); + d2_eval_ring(T, &tree.u, tree.q, x, bits); + + mpf_div(y, x, fh2); + if (!tiny_quadratic_root(D2_L1, D2_L2, v, bits)) { + return 0; + } + + mpf_mul_ui(z, y, 1728UL); + mpf_ui_sub(tmp, 1UL, z); + mpf_mul(tmp, tmp, v); + mpf_mul_ui(tmp, tmp, 427UL); + mpf_div(alpha, y, tmp); + mpf_ui_sub(beta, 1UL, alpha); + + mpf_set_ui(K, 427UL); + mpf_sqrt(K, K); + mpf_div_ui(K, K, 6UL); + mpf_ui_sub(tmp, 1UL, z); + mpf_sqrt(tmp, tmp); + mpf_mul(K, K, tmp); + + mpf_mul(inv_pi, beta, F); + mpf_mul_ui(tmp, T, 6UL); + mpf_add(inv_pi, inv_pi, tmp); + mpf_mul(inv_pi, inv_pi, K); + mpf_ui_div(pi, 1UL, inv_pi); + + int ok = write_pi(path, pi, digits); + if (seconds) *seconds = now_seconds() - t0; + if (terms_out) *terms_out = terms; + + mpf_clears(fh1, fh2, fdisc, root, x, F, T, y, v, z, alpha, beta, + K, tmp, inv_pi, pi, NULL); + d2_seg_clear(&tree); + mpz_clears(h1, h2, disc, tmpz, NULL); + return ok; +} + +static int run_d2_bs_case(uint64_t digits) { + const uint64_t guard = 192ULL; + char ref_path[256], d2_path[256]; + snprintf(ref_path, sizeof(ref_path), "build/d2bs_ref_%llu.txt", + (unsigned long long)digits); + snprintf(d2_path, sizeof(d2_path), "build/d2bs_%llu.txt", + (unsigned long long)digits); + + double tref = 0.0, td2 = 0.0; + unsigned long nref = 0UL, nd2 = 0UL; + + if (!chud_bs_pi(digits, guard, ref_path, &tref, &nref)) return 0; + + /* First try the same working precision as the Chudnovsky reference. */ + if (!d2_bs_pi(digits, guard, d2_path, &td2, &nd2, 0UL)) return 0; + int eq = files_equal(ref_path, d2_path); + + printf("%llu,%.9f,%lu,%.9f,%lu,%.6f,%d\n", + (unsigned long long)digits, + tref, nref, + td2, nd2, + td2 / tref, + eq); + fflush(stdout); + return eq; +} + +int main(int argc, char **argv) { + printf("digits,chud_bs_seconds,chud_terms,d2_bs_seconds,d2_terms,d2_over_chud_bs,equal\n"); + if (argc <= 1) { + const uint64_t defaults[] = {1000ULL, 10000ULL, 30000ULL, 100000ULL, 200000ULL}; + for (size_t i = 0; i < sizeof(defaults) / sizeof(defaults[0]); ++i) { + if (!run_d2_bs_case(defaults[i])) return 1; + } + return 0; + } + for (int i = 1; i < argc; ++i) { + uint64_t digits = strtoull(argv[i], NULL, 10); + if (digits == 0ULL || !run_d2_bs_case(digits)) return 1; + } + return 0; +} From 1dfa3dde6521ed0505fbc3838ed73b0d62e1b01d Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 16:44:10 +0800 Subject: [PATCH 02/10] Run accelerated D2 binary splitting benchmark --- .../benchmark-d2-binary-splitting.yml | 29 +++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 .github/workflows/benchmark-d2-binary-splitting.yml diff --git a/.github/workflows/benchmark-d2-binary-splitting.yml b/.github/workflows/benchmark-d2-binary-splitting.yml new file mode 100644 index 0000000..28f114a --- /dev/null +++ b/.github/workflows/benchmark-d2-binary-splitting.yml @@ -0,0 +1,29 @@ +name: D2 binary splitting benchmark + +on: + pull_request: + branches: [ main ] + +jobs: + benchmark: + runs-on: ubuntu-latest + timeout-minutes: 20 + steps: + - uses: actions/checkout@v4 + - name: Install dependencies + run: sudo apt-get update && sudo apt-get install -y build-essential libgmp-dev + - name: Build benchmark + run: | + mkdir -p build + gcc -O3 -march=native -std=c17 -Wall -Wextra -Wpedantic \ + -Isrc -Ivalidation validation/benchmark_d2_bs_vs_classics.c \ + src/ramanujan_c_common.c -lgmp -lm \ + -o build/benchmark_d2_bs_vs_classics + - name: Run benchmark + run: | + build/benchmark_d2_bs_vs_classics 1000 10000 30000 100000 200000 | tee build/d2_bs_vs_chud.csv + - name: Upload benchmark table + uses: actions/upload-artifact@v4 + with: + name: d2-bs-vs-chud + path: build/d2_bs_vs_chud.csv From 85e488993a911249ba9ceccf94def5091f6f245d Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 16:48:06 +0800 Subject: [PATCH 03/10] Fuse D2 linear weight into binary splitting tree --- validation/benchmark_d2_bs_vs_classics.c | 210 ++++++++++------------- 1 file changed, 92 insertions(+), 118 deletions(-) diff --git a/validation/benchmark_d2_bs_vs_classics.c b/validation/benchmark_d2_bs_vs_classics.c index 0b08718..f47baa6 100644 --- a/validation/benchmark_d2_bs_vs_classics.c +++ b/validation/benchmark_d2_bs_vs_classics.c @@ -3,19 +3,26 @@ #undef main /* - * Accelerated comparison only. + * Accelerated D2 only. This is the class-number-2 formula with the same + * product-tree idea used by Chudnovsky binary splitting. * - * D2 uses a balanced exact product tree in the quadratic ring generated by - * x = H2/J, where x^2 + H1*x + H2 = 0. This is the direct analogue of - * Chudnovsky binary splitting: the tree topology is unchanged, while scalar - * integer numerators are replaced by degree-2 algebraic numerators. + * Let x = H2/J. Then x^2 + H1*x + H2 = 0 and z = 1728*x/H2. + * The exact linear weight also lies in Q(x): * - * The tree carries the common denominator Q, the ratio product P, the partial - * hypergeometric sum S, and its theta-weighted sum U = sum n*t_n. No modular - * polynomial, no (z,u,v) transport, and no repeated modular Newton solve is - * used. + * beta = 1-alpha + * = (166325935146432054371646921561600 - 66389*x) + * / 1799585700948647322182745887577600. + * + * Therefore the tree can carry one weighted sum W directly, rather than + * separate F and theta(F) accumulators. Each merge carries only Q, P and W, + * matching the P/Q/T structure of ordinary binary splitting as closely as a + * quadratic coefficient ring permits. */ +static const char *D2_BETA0 = "166325935146432054371646921561600"; +static const char *D2_BETA1 = "-66389"; +static const char *D2_BETAD = "1799585700948647322182745887577600"; + typedef struct { mpz_t a; mpz_t b; @@ -24,8 +31,7 @@ typedef struct { typedef struct { mpz_t q; d2_ring_t p; - d2_ring_t s; - d2_ring_t u; + d2_ring_t w; } d2_seg_t; static void d2_ring_init(d2_ring_t *r) { @@ -41,22 +47,16 @@ static void d2_ring_clear(d2_ring_t *r) { static void d2_seg_init(d2_seg_t *s) { mpz_init(s->q); d2_ring_init(&s->p); - d2_ring_init(&s->s); - d2_ring_init(&s->u); + d2_ring_init(&s->w); } static void d2_seg_clear(d2_seg_t *s) { mpz_clear(s->q); d2_ring_clear(&s->p); - d2_ring_clear(&s->s); - d2_ring_clear(&s->u); + d2_ring_clear(&s->w); } -/* - * (a+b*x)(c+d*x), x^2 = -H1*x-H2. - * Three large products are used for the bilinear part; products by H1/H2 - * have one fixed ~100-bit operand. - */ +/* Three-large-product multiplication in Z[x]/(x^2+H1*x+H2). */ static void d2_ring_mul(d2_ring_t *out, const d2_ring_t *x, const d2_ring_t *y, @@ -69,7 +69,6 @@ static void d2_ring_mul(d2_ring_t *out, mpz_t t1) { mpz_mul(m0, x->a, y->a); mpz_mul(m1, x->b, y->b); - mpz_add(t0, x->a, x->b); mpz_add(t1, y->a, y->b); mpz_mul(m2, t0, t1); @@ -85,11 +84,14 @@ static void d2_ring_mul(d2_ring_t *out, static void d2_leaf(d2_seg_t *out, unsigned long n, - const mpz_t h2) { - mpz_t num, den, g; - mpz_inits(num, den, g, NULL); - - /* r_n = t_{n+1}/t_n = A_n*1728*x / (72*(n+1)^3*H2). */ + const mpz_t h2, + const mpz_t beta0, + const mpz_t beta1, + const mpz_t betad) { + mpz_t num, den, g, weight0; + mpz_inits(num, den, g, weight0, NULL); + + /* r_n = A_n*1728*x / (72*(n+1)^3*H2). */ mpz_set_ui(num, 6UL * n + 1UL); mpz_mul_ui(num, num, 2UL * n + 1UL); mpz_mul_ui(num, num, 6UL * n + 5UL); @@ -101,7 +103,6 @@ static void d2_leaf(d2_seg_t *out, mpz_mul_ui(den, den, 72UL); mpz_mul(den, den, h2); - /* Cancel the scalar part before it enters the tree. */ mpz_gcd(g, num, den); mpz_divexact(num, num, g); mpz_divexact(out->q, den, g); @@ -109,22 +110,29 @@ static void d2_leaf(d2_seg_t *out, mpz_set_ui(out->p.a, 0UL); mpz_set(out->p.b, num); - mpz_set(out->s.a, out->q); - mpz_set_ui(out->s.b, 0UL); - - mpz_mul_ui(out->u.a, out->q, n); - mpz_set_ui(out->u.b, 0UL); - - mpz_clears(num, den, g, NULL); + /* + * W_leaf = beta + 6n. + * The segment invariant is W = w / (BETAD*q), hence multiply the exact + * ring weight by q at the leaf. + */ + mpz_mul_ui(weight0, betad, 6UL * n); + mpz_add(weight0, weight0, beta0); + mpz_mul(out->w.a, weight0, out->q); + mpz_mul(out->w.b, beta1, out->q); + + mpz_clears(num, den, g, weight0, NULL); } static void d2_bs_build(d2_seg_t *out, unsigned long a, unsigned long b, const mpz_t h1, - const mpz_t h2) { + const mpz_t h2, + const mpz_t beta0, + const mpz_t beta1, + const mpz_t betad) { if (b - a == 1UL) { - d2_leaf(out, a, h2); + d2_leaf(out, a, h2, beta0, beta1, betad); return; } @@ -132,67 +140,39 @@ static void d2_bs_build(d2_seg_t *out, d2_seg_t left, right; d2_seg_init(&left); d2_seg_init(&right); - d2_bs_build(&left, a, m, h1, h2); - d2_bs_build(&right, m, b, h1, h2); + d2_bs_build(&left, a, m, h1, h2, beta0, beta1, betad); + d2_bs_build(&right, m, b, h1, h2, beta0, beta1, betad); mpz_t m0, m1, m2, t0, t1; mpz_inits(m0, m1, m2, t0, t1, NULL); - d2_ring_t rs, ru; - d2_ring_init(&rs); - d2_ring_init(&ru); + d2_ring_t rw; + d2_ring_init(&rw); mpz_mul(out->q, left.q, right.q); - d2_ring_mul(&out->p, &left.p, &right.p, h1, h2, m0, m1, m2, t0, t1); - d2_ring_mul(&rs, &left.p, &right.s, + /* W = W_L + P_L*W_R. */ + d2_ring_mul(&rw, &left.p, &right.w, h1, h2, m0, m1, m2, t0, t1); - mpz_mul(out->s.a, left.s.a, right.q); - mpz_add(out->s.a, out->s.a, rs.a); - mpz_mul(out->s.b, left.s.b, right.q); - mpz_add(out->s.b, out->s.b, rs.b); + mpz_mul(out->w.a, left.w.a, right.q); + mpz_add(out->w.a, out->w.a, rw.a); + mpz_mul(out->w.b, left.w.b, right.q); + mpz_add(out->w.b, out->w.b, rw.b); - d2_ring_mul(&ru, &left.p, &right.u, - h1, h2, m0, m1, m2, t0, t1); - mpz_mul(out->u.a, left.u.a, right.q); - mpz_add(out->u.a, out->u.a, ru.a); - mpz_mul(out->u.b, left.u.b, right.q); - mpz_add(out->u.b, out->u.b, ru.b); - - d2_ring_clear(&rs); - d2_ring_clear(&ru); + d2_ring_clear(&rw); mpz_clears(m0, m1, m2, t0, t1, NULL); d2_seg_clear(&left); d2_seg_clear(&right); } -static void d2_eval_ring(mpf_t out, - const d2_ring_t *r, - const mpz_t q, - const mpf_t x, - mp_bitcnt_t bits) { - mpf_t a, b, fq; - mpf_init2(a, bits); - mpf_init2(b, bits); - mpf_init2(fq, bits); - mpf_set_z(a, r->a); - mpf_set_z(b, r->b); - mpf_mul(b, b, x); - mpf_add(a, a, b); - mpf_set_z(fq, q); - mpf_div(out, a, fq); - mpf_clears(a, b, fq, NULL); -} - static int d2_bs_pi(uint64_t digits, uint64_t guard, const char *path, double *seconds, - unsigned long *terms_out, - mp_bitcnt_t extra_bits) { + unsigned long *terms_out) { double t0 = now_seconds(); - mp_bitcnt_t bits = rj_bits_for_decimal(digits, guard) + extra_bits; + mp_bitcnt_t bits = rj_bits_for_decimal(digits, guard); const long double digits_per_term = 24.95589965765426673091470594098813130578L; @@ -200,34 +180,37 @@ static int d2_bs_pi(uint64_t digits, (unsigned long)ceill(((long double)digits + (long double)guard + 30.0L) / digits_per_term) + 2UL; - mpz_t h1, h2, disc, tmpz; - mpz_inits(h1, h2, disc, tmpz, NULL); + mpz_t h1, h2, beta0, beta1, betad, disc, tmpz; + mpz_inits(h1, h2, beta0, beta1, betad, disc, tmpz, NULL); if (mpz_set_str(h1, D2_H1, 10) != 0 || - mpz_set_str(h2, D2_H2, 10) != 0) { - mpz_clears(h1, h2, disc, tmpz, NULL); + mpz_set_str(h2, D2_H2, 10) != 0 || + mpz_set_str(beta0, D2_BETA0, 10) != 0 || + mpz_set_str(beta1, D2_BETA1, 10) != 0 || + mpz_set_str(betad, D2_BETAD, 10) != 0) { + mpz_clears(h1, h2, beta0, beta1, betad, disc, tmpz, NULL); return 0; } d2_seg_t tree; d2_seg_init(&tree); - d2_bs_build(&tree, 0UL, terms, h1, h2); + d2_bs_build(&tree, 0UL, terms, h1, h2, beta0, beta1, betad); - /* Stable small embedding x = H2/J = -2 H2/(H1+sqrt(H1^2-4H2)). */ mpz_mul(disc, h1, h1); mpz_mul_ui(tmpz, h2, 4UL); mpz_sub(disc, disc, tmpz); - mpf_t fh1, fh2, fdisc, root, x, F, T, y, v, z, alpha, beta, - K, tmp, inv_pi, pi; - mpf_init2(fh1, bits); mpf_init2(fh2, bits); + mpf_t fh1, fh2, fdisc, root, x, fw0, fw1, fq, fbd, + W, y, z, K, tmp, inv_pi, pi; + mpf_init2(fh1, bits); mpf_init2(fh2, bits); mpf_init2(fdisc, bits); mpf_init2(root, bits); - mpf_init2(x, bits); mpf_init2(F, bits); - mpf_init2(T, bits); mpf_init2(y, bits); - mpf_init2(v, bits); mpf_init2(z, bits); - mpf_init2(alpha, bits); mpf_init2(beta, bits); - mpf_init2(K, bits); mpf_init2(tmp, bits); - mpf_init2(inv_pi, bits);mpf_init2(pi, bits); - + mpf_init2(x, bits); mpf_init2(fw0, bits); + mpf_init2(fw1, bits); mpf_init2(fq, bits); + mpf_init2(fbd, bits); mpf_init2(W, bits); + mpf_init2(y, bits); mpf_init2(z, bits); + mpf_init2(K, bits); mpf_init2(tmp, bits); + mpf_init2(inv_pi, bits); mpf_init2(pi, bits); + + /* Stable small embedding x = -2H2/(H1+sqrt(H1^2-4H2)). */ mpf_set_z(fh1, h1); mpf_set_z(fh2, h2); mpf_set_z(fdisc, disc); @@ -237,20 +220,18 @@ static int d2_bs_pi(uint64_t digits, mpf_neg(x, x); mpf_div(x, x, tmp); - d2_eval_ring(F, &tree.s, tree.q, x, bits); - d2_eval_ring(T, &tree.u, tree.q, x, bits); + /* W = (w.a + w.b*x)/(BETAD*q). */ + mpf_set_z(fw0, tree.w.a); + mpf_set_z(fw1, tree.w.b); + mpf_mul(fw1, fw1, x); + mpf_add(W, fw0, fw1); + mpf_set_z(fq, tree.q); + mpf_div(W, W, fq); + mpf_set_z(fbd, betad); + mpf_div(W, W, fbd); mpf_div(y, x, fh2); - if (!tiny_quadratic_root(D2_L1, D2_L2, v, bits)) { - return 0; - } - mpf_mul_ui(z, y, 1728UL); - mpf_ui_sub(tmp, 1UL, z); - mpf_mul(tmp, tmp, v); - mpf_mul_ui(tmp, tmp, 427UL); - mpf_div(alpha, y, tmp); - mpf_ui_sub(beta, 1UL, alpha); mpf_set_ui(K, 427UL); mpf_sqrt(K, K); @@ -259,20 +240,17 @@ static int d2_bs_pi(uint64_t digits, mpf_sqrt(tmp, tmp); mpf_mul(K, K, tmp); - mpf_mul(inv_pi, beta, F); - mpf_mul_ui(tmp, T, 6UL); - mpf_add(inv_pi, inv_pi, tmp); - mpf_mul(inv_pi, inv_pi, K); + mpf_mul(inv_pi, K, W); mpf_ui_div(pi, 1UL, inv_pi); int ok = write_pi(path, pi, digits); if (seconds) *seconds = now_seconds() - t0; if (terms_out) *terms_out = terms; - mpf_clears(fh1, fh2, fdisc, root, x, F, T, y, v, z, alpha, beta, - K, tmp, inv_pi, pi, NULL); + mpf_clears(fh1, fh2, fdisc, root, x, fw0, fw1, fq, fbd, + W, y, z, K, tmp, inv_pi, pi, NULL); d2_seg_clear(&tree); - mpz_clears(h1, h2, disc, tmpz, NULL); + mpz_clears(h1, h2, beta0, beta1, betad, disc, tmpz, NULL); return ok; } @@ -286,13 +264,10 @@ static int run_d2_bs_case(uint64_t digits) { double tref = 0.0, td2 = 0.0; unsigned long nref = 0UL, nd2 = 0UL; + if (!chud_bs_pi(digits, guard, ref_path, &tref, &nref) || + !d2_bs_pi(digits, guard, d2_path, &td2, &nd2)) return 0; - if (!chud_bs_pi(digits, guard, ref_path, &tref, &nref)) return 0; - - /* First try the same working precision as the Chudnovsky reference. */ - if (!d2_bs_pi(digits, guard, d2_path, &td2, &nd2, 0UL)) return 0; int eq = files_equal(ref_path, d2_path); - printf("%llu,%.9f,%lu,%.9f,%lu,%.6f,%d\n", (unsigned long long)digits, tref, nref, @@ -307,9 +282,8 @@ int main(int argc, char **argv) { printf("digits,chud_bs_seconds,chud_terms,d2_bs_seconds,d2_terms,d2_over_chud_bs,equal\n"); if (argc <= 1) { const uint64_t defaults[] = {1000ULL, 10000ULL, 30000ULL, 100000ULL, 200000ULL}; - for (size_t i = 0; i < sizeof(defaults) / sizeof(defaults[0]); ++i) { + for (size_t i = 0; i < sizeof(defaults) / sizeof(defaults[0]); ++i) if (!run_d2_bs_case(defaults[i])) return 1; - } return 0; } for (int i = 1; i < argc; ++i) { From dfe1570cdf5c368eb9ecc6de4a0c0058d75f96e1 Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 16:50:54 +0800 Subject: [PATCH 04/10] Use integral quadratic basis for D2 binary splitting --- validation/benchmark_d2_bs_vs_classics.c | 186 +++++++++++------------ 1 file changed, 91 insertions(+), 95 deletions(-) diff --git a/validation/benchmark_d2_bs_vs_classics.c b/validation/benchmark_d2_bs_vs_classics.c index f47baa6..b305244 100644 --- a/validation/benchmark_d2_bs_vs_classics.c +++ b/validation/benchmark_d2_bs_vs_classics.c @@ -3,25 +3,29 @@ #undef main /* - * Accelerated D2 only. This is the class-number-2 formula with the same - * product-tree idea used by Chudnovsky binary splitting. + * Accelerated D2 only. * - * Let x = H2/J. Then x^2 + H1*x + H2 = 0 and z = 1728*x/H2. - * The exact linear weight also lies in Q(x): + * Instead of the large-coefficient x=H2/J basis, use the integral basis + * omega = (1 + sqrt(61))/2, omega^2 = omega + 15. * - * beta = 1-alpha - * = (166325935146432054371646921561600 - 66389*x) - * / 1799585700948647322182745887577600. + * For the selected D=-427 branch the two algebraic quantities needed by the + * explicit formula become * - * Therefore the tree can carry one weighted sum W directly, rather than - * separate F and theta(F) accumulators. Each merge carries only Q, P and W, - * matching the P/Q/T structure of ordinary binary splitting as closely as a - * quadratic coefficient ring permits. + * z = (Z0 + Z1*omega)/ZD + * beta = 1-alpha = (B0 + B1*omega)/BD + * + * with the small exact coefficients below. This keeps the same binary- + * splitting/product-tree skeleton as Chudnovsky while avoiding the large + * H1/H2 coefficients in every quadratic-ring multiplication. */ -static const char *D2_BETA0 = "166325935146432054371646921561600"; -static const char *D2_BETA1 = "-66389"; -static const char *D2_BETAD = "1799585700948647322182745887577600"; +static const char *D2_Z0 = "-59818419102592333"; +static const char *D2_Z1 = "13579278976889262"; +static const char *D2_ZD = "609541351191872000"; + +static const char *D2_B0 = "320004750671"; +static const char *D2_B1 = "-56552805760"; +static const char *D2_BD = "766923569391"; typedef struct { mpz_t a; @@ -56,12 +60,15 @@ static void d2_seg_clear(d2_seg_t *s) { d2_ring_clear(&s->w); } -/* Three-large-product multiplication in Z[x]/(x^2+H1*x+H2). */ +/* + * (a+b*w)(c+d*w), w^2=w+15. + * m0=ac, m1=bd, m2=(a+b)(c+d), hence + * constant = m0 + 15*m1 + * w-coeff = m2 - m0. + */ static void d2_ring_mul(d2_ring_t *out, const d2_ring_t *x, const d2_ring_t *y, - const mpz_t h1, - const mpz_t h2, mpz_t m0, mpz_t m1, mpz_t m2, @@ -73,52 +80,45 @@ static void d2_ring_mul(d2_ring_t *out, mpz_add(t1, y->a, y->b); mpz_mul(m2, t0, t1); - mpz_mul(t0, h2, m1); - mpz_sub(out->a, m0, t0); - - mpz_sub(m2, m2, m0); - mpz_sub(m2, m2, m1); - mpz_mul(t0, h1, m1); - mpz_sub(out->b, m2, t0); + mpz_mul_ui(t0, m1, 15UL); + mpz_add(out->a, m0, t0); + mpz_sub(out->b, m2, m0); } static void d2_leaf(d2_seg_t *out, unsigned long n, - const mpz_t h2, - const mpz_t beta0, - const mpz_t beta1, - const mpz_t betad) { + const mpz_t z0, + const mpz_t z1, + const mpz_t zd, + const mpz_t b0, + const mpz_t b1, + const mpz_t bd) { mpz_t num, den, g, weight0; mpz_inits(num, den, g, weight0, NULL); - /* r_n = A_n*1728*x / (72*(n+1)^3*H2). */ + /* r_n = A_n*z / (72*(n+1)^3). */ mpz_set_ui(num, 6UL * n + 1UL); mpz_mul_ui(num, num, 2UL * n + 1UL); mpz_mul_ui(num, num, 6UL * n + 5UL); - mpz_mul_ui(num, num, 1728UL); mpz_set_ui(den, n + 1UL); mpz_mul_ui(den, den, n + 1UL); mpz_mul_ui(den, den, n + 1UL); mpz_mul_ui(den, den, 72UL); - mpz_mul(den, den, h2); + mpz_mul(den, den, zd); mpz_gcd(g, num, den); mpz_divexact(num, num, g); mpz_divexact(out->q, den, g); - mpz_set_ui(out->p.a, 0UL); - mpz_set(out->p.b, num); + mpz_mul(out->p.a, z0, num); + mpz_mul(out->p.b, z1, num); - /* - * W_leaf = beta + 6n. - * The segment invariant is W = w / (BETAD*q), hence multiply the exact - * ring weight by q at the leaf. - */ - mpz_mul_ui(weight0, betad, 6UL * n); - mpz_add(weight0, weight0, beta0); + /* Segment invariant: W = w/(BD*q), with W_leaf=beta+6n. */ + mpz_mul_ui(weight0, bd, 6UL * n); + mpz_add(weight0, weight0, b0); mpz_mul(out->w.a, weight0, out->q); - mpz_mul(out->w.b, beta1, out->q); + mpz_mul(out->w.b, b1, out->q); mpz_clears(num, den, g, weight0, NULL); } @@ -126,13 +126,14 @@ static void d2_leaf(d2_seg_t *out, static void d2_bs_build(d2_seg_t *out, unsigned long a, unsigned long b, - const mpz_t h1, - const mpz_t h2, - const mpz_t beta0, - const mpz_t beta1, - const mpz_t betad) { + const mpz_t z0, + const mpz_t z1, + const mpz_t zd, + const mpz_t b0, + const mpz_t b1, + const mpz_t bd) { if (b - a == 1UL) { - d2_leaf(out, a, h2, beta0, beta1, betad); + d2_leaf(out, a, z0, z1, zd, b0, b1, bd); return; } @@ -140,8 +141,8 @@ static void d2_bs_build(d2_seg_t *out, d2_seg_t left, right; d2_seg_init(&left); d2_seg_init(&right); - d2_bs_build(&left, a, m, h1, h2, beta0, beta1, betad); - d2_bs_build(&right, m, b, h1, h2, beta0, beta1, betad); + d2_bs_build(&left, a, m, z0, z1, zd, b0, b1, bd); + d2_bs_build(&right, m, b, z0, z1, zd, b0, b1, bd); mpz_t m0, m1, m2, t0, t1; mpz_inits(m0, m1, m2, t0, t1, NULL); @@ -149,12 +150,10 @@ static void d2_bs_build(d2_seg_t *out, d2_ring_init(&rw); mpz_mul(out->q, left.q, right.q); - d2_ring_mul(&out->p, &left.p, &right.p, - h1, h2, m0, m1, m2, t0, t1); + d2_ring_mul(&out->p, &left.p, &right.p, m0, m1, m2, t0, t1); /* W = W_L + P_L*W_R. */ - d2_ring_mul(&rw, &left.p, &right.w, - h1, h2, m0, m1, m2, t0, t1); + d2_ring_mul(&rw, &left.p, &right.w, m0, m1, m2, t0, t1); mpz_mul(out->w.a, left.w.a, right.q); mpz_add(out->w.a, out->w.a, rw.a); mpz_mul(out->w.b, left.w.b, right.q); @@ -180,59 +179,56 @@ static int d2_bs_pi(uint64_t digits, (unsigned long)ceill(((long double)digits + (long double)guard + 30.0L) / digits_per_term) + 2UL; - mpz_t h1, h2, beta0, beta1, betad, disc, tmpz; - mpz_inits(h1, h2, beta0, beta1, betad, disc, tmpz, NULL); - if (mpz_set_str(h1, D2_H1, 10) != 0 || - mpz_set_str(h2, D2_H2, 10) != 0 || - mpz_set_str(beta0, D2_BETA0, 10) != 0 || - mpz_set_str(beta1, D2_BETA1, 10) != 0 || - mpz_set_str(betad, D2_BETAD, 10) != 0) { - mpz_clears(h1, h2, beta0, beta1, betad, disc, tmpz, NULL); + mpz_t z0, z1, zd, b0, b1, bd; + mpz_inits(z0, z1, zd, b0, b1, bd, NULL); + if (mpz_set_str(z0, D2_Z0, 10) != 0 || + mpz_set_str(z1, D2_Z1, 10) != 0 || + mpz_set_str(zd, D2_ZD, 10) != 0 || + mpz_set_str(b0, D2_B0, 10) != 0 || + mpz_set_str(b1, D2_B1, 10) != 0 || + mpz_set_str(bd, D2_BD, 10) != 0) { + mpz_clears(z0, z1, zd, b0, b1, bd, NULL); return 0; } d2_seg_t tree; d2_seg_init(&tree); - d2_bs_build(&tree, 0UL, terms, h1, h2, beta0, beta1, betad); - - mpz_mul(disc, h1, h1); - mpz_mul_ui(tmpz, h2, 4UL); - mpz_sub(disc, disc, tmpz); + d2_bs_build(&tree, 0UL, terms, z0, z1, zd, b0, b1, bd); - mpf_t fh1, fh2, fdisc, root, x, fw0, fw1, fq, fbd, - W, y, z, K, tmp, inv_pi, pi; - mpf_init2(fh1, bits); mpf_init2(fh2, bits); - mpf_init2(fdisc, bits); mpf_init2(root, bits); - mpf_init2(x, bits); mpf_init2(fw0, bits); + mpf_t omega, root61, fz0, fz1, fzd, fw0, fw1, fq, fbd, + z, W, K, tmp, inv_pi, pi; + mpf_init2(omega, bits); mpf_init2(root61, bits); + mpf_init2(fz0, bits); mpf_init2(fz1, bits); + mpf_init2(fzd, bits); mpf_init2(fw0, bits); mpf_init2(fw1, bits); mpf_init2(fq, bits); - mpf_init2(fbd, bits); mpf_init2(W, bits); - mpf_init2(y, bits); mpf_init2(z, bits); - mpf_init2(K, bits); mpf_init2(tmp, bits); - mpf_init2(inv_pi, bits); mpf_init2(pi, bits); - - /* Stable small embedding x = -2H2/(H1+sqrt(H1^2-4H2)). */ - mpf_set_z(fh1, h1); - mpf_set_z(fh2, h2); - mpf_set_z(fdisc, disc); - mpf_sqrt(root, fdisc); - mpf_add(tmp, fh1, root); - mpf_mul_ui(x, fh2, 2UL); - mpf_neg(x, x); - mpf_div(x, x, tmp); - - /* W = (w.a + w.b*x)/(BETAD*q). */ + mpf_init2(fbd, bits); mpf_init2(z, bits); + mpf_init2(W, bits); mpf_init2(K, bits); + mpf_init2(tmp, bits); mpf_init2(inv_pi, bits); + mpf_init2(pi, bits); + + mpf_set_ui(root61, 61UL); + mpf_sqrt(root61, root61); + mpf_add_ui(omega, root61, 1UL); + mpf_div_ui(omega, omega, 2UL); + + /* z = (Z0 + Z1*omega)/ZD. */ + mpf_set_z(fz0, z0); + mpf_set_z(fz1, z1); + mpf_mul(fz1, fz1, omega); + mpf_add(z, fz0, fz1); + mpf_set_z(fzd, zd); + mpf_div(z, z, fzd); + + /* W = (w.a + w.b*omega)/(BD*q). */ mpf_set_z(fw0, tree.w.a); mpf_set_z(fw1, tree.w.b); - mpf_mul(fw1, fw1, x); + mpf_mul(fw1, fw1, omega); mpf_add(W, fw0, fw1); mpf_set_z(fq, tree.q); mpf_div(W, W, fq); - mpf_set_z(fbd, betad); + mpf_set_z(fbd, bd); mpf_div(W, W, fbd); - mpf_div(y, x, fh2); - mpf_mul_ui(z, y, 1728UL); - mpf_set_ui(K, 427UL); mpf_sqrt(K, K); mpf_div_ui(K, K, 6UL); @@ -247,10 +243,10 @@ static int d2_bs_pi(uint64_t digits, if (seconds) *seconds = now_seconds() - t0; if (terms_out) *terms_out = terms; - mpf_clears(fh1, fh2, fdisc, root, x, fw0, fw1, fq, fbd, - W, y, z, K, tmp, inv_pi, pi, NULL); + mpf_clears(omega, root61, fz0, fz1, fzd, fw0, fw1, fq, fbd, + z, W, K, tmp, inv_pi, pi, NULL); d2_seg_clear(&tree); - mpz_clears(h1, h2, beta0, beta1, betad, disc, tmpz, NULL); + mpz_clears(z0, z1, zd, b0, b1, bd, NULL); return ok; } From 89aa952d86a3aead88a630fb7effdc04c00d8f5f Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 16:52:45 +0800 Subject: [PATCH 05/10] Add accelerated classical comparison without bare series runs --- validation/benchmark_agm_borwein_only.c | 120 ++++++++++++++++++++++++ 1 file changed, 120 insertions(+) create mode 100644 validation/benchmark_agm_borwein_only.c diff --git a/validation/benchmark_agm_borwein_only.c b/validation/benchmark_agm_borwein_only.c new file mode 100644 index 0000000..a30ae43 --- /dev/null +++ b/validation/benchmark_agm_borwein_only.c @@ -0,0 +1,120 @@ +#define main original_d2_benchmark_main +#include "benchmark_original_d2_basic_vs_chud.c" +#undef main + +/* Classical iterative references only. The bare direct D2/Chud series are + * compiled only to reuse common helpers and are never executed here. */ + +static unsigned long ceil_log2_u64_fast(uint64_t n) { + unsigned long k = 0UL; + uint64_t x = 1ULL; + while (x < n && k < 63UL) { + x <<= 1; + ++k; + } + return k; +} + +static int agm_pi_only(uint64_t digits, uint64_t guard, const char *path, + double *seconds, unsigned long *iterations_out) { + double t0 = now_seconds(); + mp_bitcnt_t bits = rj_bits_for_decimal(digits, guard); + mpf_t a,b,t,p,an,bn,delta,tmp,sum,pi; + mpf_init2(a,bits); mpf_init2(b,bits); mpf_init2(t,bits); + mpf_init2(p,bits); mpf_init2(an,bits); mpf_init2(bn,bits); + mpf_init2(delta,bits); mpf_init2(tmp,bits); mpf_init2(sum,bits); + mpf_init2(pi,bits); + + mpf_set_ui(a,1UL); + mpf_set_ui(tmp,2UL); mpf_sqrt(tmp,tmp); mpf_ui_div(b,1UL,tmp); + mpf_set_ui(t,1UL); mpf_div_ui(t,t,4UL); + mpf_set_ui(p,1UL); + + unsigned long iterations = ceil_log2_u64_fast(digits + guard + 1ULL) + 4UL; + for (unsigned long i=0UL;i Date: Mon, 24 Aug 2026 16:53:06 +0800 Subject: [PATCH 06/10] Compare accelerated D2 with AGM and Borwein on one runner --- .../benchmark-d2-binary-splitting.yml | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/.github/workflows/benchmark-d2-binary-splitting.yml b/.github/workflows/benchmark-d2-binary-splitting.yml index 28f114a..b3ca1ab 100644 --- a/.github/workflows/benchmark-d2-binary-splitting.yml +++ b/.github/workflows/benchmark-d2-binary-splitting.yml @@ -12,18 +12,27 @@ jobs: - uses: actions/checkout@v4 - name: Install dependencies run: sudo apt-get update && sudo apt-get install -y build-essential libgmp-dev - - name: Build benchmark + - name: Build accelerated benchmarks run: | mkdir -p build gcc -O3 -march=native -std=c17 -Wall -Wextra -Wpedantic \ -Isrc -Ivalidation validation/benchmark_d2_bs_vs_classics.c \ src/ramanujan_c_common.c -lgmp -lm \ -o build/benchmark_d2_bs_vs_classics - - name: Run benchmark + gcc -O3 -march=native -std=c17 -Wall -Wextra -Wpedantic \ + -Isrc -Ivalidation validation/benchmark_agm_borwein_only.c \ + src/ramanujan_c_common.c -lgmp -lm \ + -o build/benchmark_agm_borwein_only + - name: Run accelerated D2 and Chudnovsky run: | build/benchmark_d2_bs_vs_classics 1000 10000 30000 100000 200000 | tee build/d2_bs_vs_chud.csv - - name: Upload benchmark table + - name: Run AGM and Borwein against the same reference files + run: | + build/benchmark_agm_borwein_only 1000 10000 30000 100000 200000 | tee build/agm_borwein.csv + - name: Upload benchmark tables uses: actions/upload-artifact@v4 with: - name: d2-bs-vs-chud - path: build/d2_bs_vs_chud.csv + name: accelerated-pi-comparison + path: | + build/d2_bs_vs_chud.csv + build/agm_borwein.csv From 434e13414bc403447cfcb518116601fe1ef82f64 Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 16:54:24 +0800 Subject: [PATCH 07/10] Extend accelerated D2 comparison to one million digits --- .github/workflows/benchmark-d2-binary-splitting.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/benchmark-d2-binary-splitting.yml b/.github/workflows/benchmark-d2-binary-splitting.yml index b3ca1ab..86d8c6c 100644 --- a/.github/workflows/benchmark-d2-binary-splitting.yml +++ b/.github/workflows/benchmark-d2-binary-splitting.yml @@ -25,10 +25,10 @@ jobs: -o build/benchmark_agm_borwein_only - name: Run accelerated D2 and Chudnovsky run: | - build/benchmark_d2_bs_vs_classics 1000 10000 30000 100000 200000 | tee build/d2_bs_vs_chud.csv + build/benchmark_d2_bs_vs_classics 1000 10000 30000 100000 200000 300000 1000000 | tee build/d2_bs_vs_chud.csv - name: Run AGM and Borwein against the same reference files run: | - build/benchmark_agm_borwein_only 1000 10000 30000 100000 200000 | tee build/agm_borwein.csv + build/benchmark_agm_borwein_only 1000 10000 30000 100000 200000 300000 1000000 | tee build/agm_borwein.csv - name: Upload benchmark tables uses: actions/upload-artifact@v4 with: From 3ad87dea10a22b1095c65f4168e62a105f5dd6ed Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 16:55:47 +0800 Subject: [PATCH 08/10] Add accelerated D2 binary splitting benchmark target --- Makefile | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/Makefile b/Makefile index edcfba4..e3a9f1e 100644 --- a/Makefile +++ b/Makefile @@ -25,7 +25,7 @@ PI1000_EXPLICIT := $(BUILD)/pi1000_explicit.txt PI1000_TRACE := $(BUILD)/pi1000_newton_trace.txt EXPECTED_SHA256 := e898fea26734a6d3af5396b9f4c60ae5dcc88fc40944d835911a9ee8a672ea1b -.PHONY: all clean core certify restore1000 explicit1000 trace1000 dimensionbench loworderbench d2basicbench d2classicbench smoke +.PHONY: all clean core certify restore1000 explicit1000 trace1000 dimensionbench loworderbench d2basicbench d2classicbench d2bsbench smoke all: $(BUILD)/ramanujan_core $(BUILD)/ramanujan_certifier $(BUILD)/ramanujan_restore $(BUILD)/ramanujan_explicit_pi @@ -56,6 +56,12 @@ $(BUILD)/ramanujan_d2_basic_bench: validation/benchmark_original_d2_basic_vs_chu $(BUILD)/ramanujan_d2_classic_bench: validation/benchmark_d2_vs_agm_borwein.c validation/benchmark_original_d2_basic_vs_chud.c $(COMMON) $(HDR) | $(BUILD) $(CC) $(CFLAGS) -I$(SRC) -Ivalidation validation/benchmark_d2_vs_agm_borwein.c $(COMMON) $(LDLIBS_COMMON) -o $@ +$(BUILD)/ramanujan_d2_bs_bench: validation/benchmark_d2_bs_vs_classics.c validation/benchmark_original_d2_basic_vs_chud.c $(COMMON) $(HDR) | $(BUILD) + $(CC) $(CFLAGS) -I$(SRC) -Ivalidation validation/benchmark_d2_bs_vs_classics.c $(COMMON) $(LDLIBS_COMMON) -o $@ + +$(BUILD)/ramanujan_agm_borwein_only_bench: validation/benchmark_agm_borwein_only.c validation/benchmark_original_d2_basic_vs_chud.c $(COMMON) $(HDR) | $(BUILD) + $(CC) $(CFLAGS) -I$(SRC) -Ivalidation validation/benchmark_agm_borwein_only.c $(COMMON) $(LDLIBS_COMMON) -o $@ + core: $(BUILD)/ramanujan_core $(BUILD)/ramanujan_core @@ -83,6 +89,10 @@ d2basicbench: $(BUILD)/ramanujan_d2_basic_bench d2classicbench: $(BUILD)/ramanujan_d2_classic_bench $(BUILD)/ramanujan_d2_classic_bench 1000 10000 30000 100000 200000 +d2bsbench: $(BUILD)/ramanujan_d2_bs_bench $(BUILD)/ramanujan_agm_borwein_only_bench + $(BUILD)/ramanujan_d2_bs_bench 1000 10000 30000 100000 200000 300000 1000000 + $(BUILD)/ramanujan_agm_borwein_only_bench 1000 10000 30000 100000 200000 300000 1000000 + smoke: restore1000 explicit1000 @restored=$$(sha256sum $(PI1000) | awk '{print $$1}'); \ explicit=$$(sha256sum $(PI1000_EXPLICIT) | awk '{print $$1}'); \ From f5fe944ce5cbb12690213ec73e612600ab40601a Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 16:56:10 +0800 Subject: [PATCH 09/10] Record accelerated D2 binary splitting results --- validation/d2_binary_splitting_results.md | 68 +++++++++++++++++++++++ 1 file changed, 68 insertions(+) create mode 100644 validation/d2_binary_splitting_results.md diff --git a/validation/d2_binary_splitting_results.md b/validation/d2_binary_splitting_results.md new file mode 100644 index 0000000..d1a1761 --- /dev/null +++ b/validation/d2_binary_splitting_results.md @@ -0,0 +1,68 @@ +# D2 binary splitting vs classical pi algorithms + +This benchmark removes the bare direct-series comparison from the active timing path. The D2 formula is accelerated with a balanced exact product tree, directly analogous to Chudnovsky binary splitting. + +## Algebraic basis + +For the selected class-number-2 branch (`Delta = -427`), use the integral quadratic basis + +`omega = (1 + sqrt(61))/2`, so `omega^2 = omega + 15`. + +The two algebraic quantities required by the explicit formula reduce to + +`z = (-59818419102592333 + 13579278976889262*omega) / 609541351191872000` + +and + +`beta = 1-alpha = (320004750671 - 56552805760*omega) / 766923569391`. + +This basis is materially cheaper than the earlier large-coefficient `x = H2/J` basis. A quadratic-ring multiplication + +`(a+b*omega)(c+d*omega)` + +needs three large products: + +`m0 = a*c`, `m1 = b*d`, `m2 = (a+b)(c+d)`, + +followed by + +`constant = m0 + 15*m1`, `omega_coeff = m2 - m0`. + +The binary-splitting tree carries only `Q`, ratio product `P`, and the already-weighted partial sum `W`, so the structure is the quadratic-ring analogue of the usual Chudnovsky `P/Q/T` tree. + +## Protocol + +GitHub-hosted Ubuntu 24.04, GMP 6.3.0, GCC C17, `-O3 -march=native`. + +The D2-BS and Chudnovsky-BS timings are end-to-end and include initialization, computation, decimal serialization, and file close. AGM and Borwein quartic are run immediately afterward on the same runner with the same precision budget and serializer. Their outputs are compared against the exact Chudnovsky-BS reference files produced in the first phase. + +The final comparison does not execute the bare direct D2 or bare direct Chudnovsky series. + +AGM and Borwein quartic use full target precision for every iteration; no adaptive precision-doubling schedule is used, so their timings are conservative rather than maximally optimized. + +Every output at every tested precision was byte-identical to the Chudnovsky-BS reference. + +| digits | Chud BS | D2 BS | D2 / Chud | AGM | Borwein-4 | +| ---: | ---: | ---: | ---: | ---: | ---: | +| 1,000 | 0.000127 s | 0.000139 s | 1.09x | 0.000165 s | 0.000146 s | +| 10,000 | 0.001029 s | 0.001718 s | 1.67x | 0.002711 s | 0.003610 s | +| 30,000 | 0.004262 s | 0.007249 s | 1.70x | 0.013150 s | 0.019744 s | +| 100,000 | 0.022866 s | 0.037463 s | 1.64x | 0.076566 s | 0.115266 s | +| 200,000 | 0.059652 s | 0.095458 s | 1.60x | 0.201342 s | 0.280549 s | +| 300,000 | 0.104132 s | 0.163471 s | 1.57x | 0.326274 s | 0.471700 s | +| 1,000,000 | 0.464147 s | 0.736330 s | 1.59x | 1.463959 s | 2.043881 s | + +Term counts for the one-million-digit run were 70,530 for Chudnovsky and 40,082 for D2, preserving the approximately 0.568 term-count ratio implied by the formula-level convergence constants. + +## Result + +After receiving the same class of product-tree acceleration, D2 is no longer hundreds of times slower than Chudnovsky. At one million digits it is about 1.59x slower than the current Chudnovsky binary-splitting implementation while using only about 56.8% as many series terms. + +On this same-run benchmark, D2-BS is faster than the full-precision Gauss-Legendre AGM and Borwein quartic implementations at every tested precision above the tiny startup regime. At one million digits: + +- D2-BS / AGM = about 0.503, so D2-BS is about 1.99x faster. +- D2-BS / Borwein-4 = about 0.360, so D2-BS is about 2.78x faster. + +These AGM/Borwein comparisons are implementation results, not a claim that D2-BS beats every optimized variant of those algorithms. Adaptive-precision implementations can reduce their cost. + +The remaining gap to Chudnovsky is therefore no longer caused by a missing asymptotically efficient summation scheme. It is now primarily the constant-factor cost of performing the product tree in a quadratic coefficient ring instead of the ordinary integer ring. The current measured gap is approximately 1.6x at 100k-1M digits. From 8379d7cd8fe7e6669fea35203b874fba0f7f73cf Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 16:57:48 +0800 Subject: [PATCH 10/10] Remove temporary D2 binary splitting benchmark workflow --- .../benchmark-d2-binary-splitting.yml | 38 ------------------- 1 file changed, 38 deletions(-) delete mode 100644 .github/workflows/benchmark-d2-binary-splitting.yml diff --git a/.github/workflows/benchmark-d2-binary-splitting.yml b/.github/workflows/benchmark-d2-binary-splitting.yml deleted file mode 100644 index 86d8c6c..0000000 --- a/.github/workflows/benchmark-d2-binary-splitting.yml +++ /dev/null @@ -1,38 +0,0 @@ -name: D2 binary splitting benchmark - -on: - pull_request: - branches: [ main ] - -jobs: - benchmark: - runs-on: ubuntu-latest - timeout-minutes: 20 - steps: - - uses: actions/checkout@v4 - - name: Install dependencies - run: sudo apt-get update && sudo apt-get install -y build-essential libgmp-dev - - name: Build accelerated benchmarks - run: | - mkdir -p build - gcc -O3 -march=native -std=c17 -Wall -Wextra -Wpedantic \ - -Isrc -Ivalidation validation/benchmark_d2_bs_vs_classics.c \ - src/ramanujan_c_common.c -lgmp -lm \ - -o build/benchmark_d2_bs_vs_classics - gcc -O3 -march=native -std=c17 -Wall -Wextra -Wpedantic \ - -Isrc -Ivalidation validation/benchmark_agm_borwein_only.c \ - src/ramanujan_c_common.c -lgmp -lm \ - -o build/benchmark_agm_borwein_only - - name: Run accelerated D2 and Chudnovsky - run: | - build/benchmark_d2_bs_vs_classics 1000 10000 30000 100000 200000 300000 1000000 | tee build/d2_bs_vs_chud.csv - - name: Run AGM and Borwein against the same reference files - run: | - build/benchmark_agm_borwein_only 1000 10000 30000 100000 200000 300000 1000000 | tee build/agm_borwein.csv - - name: Upload benchmark tables - uses: actions/upload-artifact@v4 - with: - name: accelerated-pi-comparison - path: | - build/d2_bs_vs_chud.csv - build/agm_borwein.csv