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}'); \ 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;ia); + 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->w); +} + +static void d2_seg_clear(d2_seg_t *s) { + mpz_clear(s->q); + d2_ring_clear(&s->p); + d2_ring_clear(&s->w); +} + +/* + * (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, + 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_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 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*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_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, zd); + + mpz_gcd(g, num, den); + mpz_divexact(num, num, g); + mpz_divexact(out->q, den, g); + + mpz_mul(out->p.a, z0, num); + mpz_mul(out->p.b, z1, num); + + /* 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, b1, 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 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, z0, z1, zd, b0, b1, bd); + 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, 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); + 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, m0, m1, m2, t0, t1); + + /* W = W_L + P_L*W_R. */ + 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); + mpz_add(out->w.b, out->w.b, rw.b); + + d2_ring_clear(&rw); + mpz_clears(m0, m1, m2, t0, t1, NULL); + d2_seg_clear(&left); + d2_seg_clear(&right); +} + +static int d2_bs_pi(uint64_t digits, + uint64_t guard, + const char *path, + double *seconds, + unsigned long *terms_out) { + double t0 = now_seconds(); + mp_bitcnt_t bits = rj_bits_for_decimal(digits, guard); + + 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 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, z0, z1, zd, b0, b1, bd); + + 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(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, omega); + mpf_add(W, fw0, fw1); + mpf_set_z(fq, tree.q); + mpf_div(W, W, fq); + mpf_set_z(fbd, bd); + mpf_div(W, W, fbd); + + 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, 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(omega, root61, fz0, fz1, fzd, fw0, fw1, fq, fbd, + z, W, K, tmp, inv_pi, pi, NULL); + d2_seg_clear(&tree); + mpz_clears(z0, z1, zd, b0, b1, bd, 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) || + !d2_bs_pi(digits, guard, d2_path, &td2, &nd2)) 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; +} 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.