diff --git a/Makefile b/Makefile index 60ea33d..25cd36f 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 d2bsbench d2parallelbench smoke +.PHONY: all clean core certify restore1000 explicit1000 trace1000 dimensionbench loworderbench d2basicbench d2classicbench d2bsbench d2parallelbench d2directtaperbench smoke all: $(BUILD)/ramanujan_core $(BUILD)/ramanujan_certifier $(BUILD)/ramanujan_restore $(BUILD)/ramanujan_explicit_pi @@ -62,6 +62,9 @@ $(BUILD)/ramanujan_d2_bs_bench: validation/benchmark_d2_bs_vs_classics.c validat $(BUILD)/ramanujan_d2_parallel_bench: validation/benchmark_d2_parallel_lanes.c validation/benchmark_original_d2_basic_vs_chud.c $(COMMON) $(HDR) | $(BUILD) $(CC) $(CFLAGS) -fopenmp -I$(SRC) -Ivalidation validation/benchmark_d2_parallel_lanes.c $(COMMON) $(LDLIBS_COMMON) -o $@ +$(BUILD)/ramanujan_d2_direct_taper_bench: validation/benchmark_direct_convergence_window.c validation/benchmark_original_d2_basic_vs_chud.c $(COMMON) $(HDR) | $(BUILD) + $(CC) $(CFLAGS) -I$(SRC) -Ivalidation validation/benchmark_direct_convergence_window.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 $@ @@ -99,6 +102,10 @@ d2bsbench: $(BUILD)/ramanujan_d2_bs_bench $(BUILD)/ramanujan_agm_borwein_only_be d2parallelbench: $(BUILD)/ramanujan_d2_parallel_bench OMP_NUM_THREADS=2 OMP_DYNAMIC=false OMP_PROC_BIND=true OMP_PLACES=cores D2_PARALLEL_THRESHOLD_BITS=4096 $(BUILD)/ramanujan_d2_parallel_bench 100000 300000 1000000 +d2directtaperbench: $(BUILD)/ramanujan_d2_direct_taper_bench + $(BUILD)/ramanujan_d2_direct_taper_bench 100000 16 + $(BUILD)/ramanujan_d2_direct_taper_bench 200000 16 + smoke: restore1000 explicit1000 @restored=$$(sha256sum $(PI1000) | awk '{print $$1}'); \ explicit=$$(sha256sum $(PI1000_EXPLICIT) | awk '{print $$1}'); \ diff --git a/noop b/noop deleted file mode 100644 index e69de29..0000000 diff --git a/validation/benchmark_direct_convergence_window.c b/validation/benchmark_direct_convergence_window.c new file mode 100644 index 0000000..14ced72 --- /dev/null +++ b/validation/benchmark_direct_convergence_window.c @@ -0,0 +1,190 @@ +#define main original_d2_direct_benchmark_main +#include "benchmark_original_d2_basic_vs_chud.c" +#undef main + +/* Convergence-Tapered Direct Recurrence (CTDR). + * + * This stays on the original scalar/direct evaluator. No binary splitting + * and no quadratic-field product tree are used. The only acceleration is + * precision tapering: term n is about 10^(-q n), so its relative precision + * can fall with n while preserving a fixed absolute error target. Precision + * is changed only at coarse block boundaries to avoid allocator churn. + */ +static mp_bitcnt_t ct_active_bits(mp_bitcnt_t max_bits, + long double digits_per_term, + unsigned long n) { + long double dropped = digits_per_term * (long double)n * 3.32192809488736234787L; + long double keep = (long double)max_bits - dropped + 1024.0L; + if (keep < 2048.0L) keep = 2048.0L; + if (keep > (long double)max_bits) keep = (long double)max_bits; + return (mp_bitcnt_t)ceill(keep); +} + +static void ct_series(const mpf_t z, + long double digits_per_term, + unsigned long terms, + mp_bitcnt_t bits, + unsigned long block_terms, + mpf_t F, + mpf_t T) { + mpf_t term, zwork, nt; + mpf_init2(term, bits); + mpf_init2(zwork, bits); + mpf_init2(nt, bits); + mpf_set_ui(term, 1UL); + mpf_set_ui(F, 1UL); + mpf_set_ui(T, 0UL); + + for (unsigned long a = 1UL; a < terms; a += block_terms) { + unsigned long b = a + block_terms; + if (b > terms) b = terms; + mp_bitcnt_t abits = ct_active_bits(bits, digits_per_term, a); + + /* Safe precision changes at block boundaries. */ + mpf_set_prec(term, abits); + mpf_set_prec(zwork, abits); + mpf_set_prec(nt, abits); + mpf_set(zwork, z); + + for (unsigned long n = a; n < b; ++n) { + hyper_term_advance(term, zwork, n); + mpf_add(F, F, term); + mpf_mul_ui(nt, term, n); + mpf_add(T, T, nt); + } + } + + mpf_clears(term, zwork, nt, NULL); +} + +static int d2_ct_pi(uint64_t digits, uint64_t guard, unsigned long block_terms, + const char *path, double *seconds, unsigned long *terms_out) { + double t0 = now_seconds(); + mp_bitcnt_t bits = rj_bits_for_decimal(digits, guard); + mpf_t y, v, z, alpha, beta, K, F, T, tmp, inv_pi, pi; + 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(F, bits); mpf_init2(T, bits); + mpf_init2(tmp, bits); mpf_init2(inv_pi, bits); + mpf_init2(pi, bits); + + if (!tiny_quadratic_root(D2_H1, D2_H2, y, bits) || + !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); + + const long double q = 24.95589965765426673091470594098813130578L; + unsigned long terms = + (unsigned long)ceill(((long double)digits + (long double)guard + 20.0L) / q) + 2UL; + ct_series(z, q, terms, bits, block_terms, F, T); + + 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(y, v, z, alpha, beta, K, F, T, tmp, inv_pi, pi, NULL); + return ok; +} + +static int chud_ct_pi(uint64_t digits, uint64_t guard, unsigned long block_terms, + const char *path, double *seconds, unsigned long *terms_out) { + double t0 = now_seconds(); + mp_bitcnt_t bits = rj_bits_for_decimal(digits, guard); + mpf_t z, alpha, beta, K, F, T, tmp, den, inv_pi, pi; + mpf_init2(z, bits); mpf_init2(alpha, bits); + mpf_init2(beta, bits); mpf_init2(K, bits); + mpf_init2(F, bits); mpf_init2(T, bits); + mpf_init2(tmp, bits); mpf_init2(den, bits); + mpf_init2(inv_pi, bits); mpf_init2(pi, bits); + + mpf_set_si(z, -1); + mpf_set_ui(den, 53360UL); + mpf_pow_ui(den, den, 3UL); + mpf_div(z, z, den); + mpf_set_ui(alpha, 77265280UL); + mpf_div_ui(alpha, alpha, 90856689UL); + mpf_ui_sub(beta, 1UL, alpha); + + mpf_set_ui(K, 163UL); + 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); + + const long double q = 14.1816474627254776555255216782L; + unsigned long terms = + (unsigned long)ceill(((long double)digits + (long double)guard + 20.0L) / q) + 2UL; + ct_series(z, q, terms, bits, block_terms, F, T); + + 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(z, alpha, beta, K, F, T, tmp, den, inv_pi, pi, NULL); + return ok; +} + +static int run_ct(uint64_t digits, unsigned long block_terms) { + const uint64_t guard = 192ULL; + char ref[256], db[256], dct[256], cb[256], cct[256]; + snprintf(ref, sizeof(ref), "build/ct_ref_%llu_%lu.txt", (unsigned long long)digits, block_terms); + snprintf(db, sizeof(db), "build/ct_d2_base_%llu_%lu.txt", (unsigned long long)digits, block_terms); + snprintf(dct, sizeof(dct), "build/ct_d2_%llu_%lu.txt", (unsigned long long)digits, block_terms); + snprintf(cb, sizeof(cb), "build/ct_chud_base_%llu_%lu.txt", (unsigned long long)digits, block_terms); + snprintf(cct, sizeof(cct), "build/ct_chud_%llu_%lu.txt", (unsigned long long)digits, block_terms); + + double tr = 0.0, tdb = 0.0, tdct = 0.0, tcb = 0.0, tcct = 0.0; + unsigned long nr = 0, ndb = 0, ndct = 0, ncb = 0, ncct = 0; + if (!chud_bs_pi(digits, guard, ref, &tr, &nr) || + !d2_basic_pi(digits, guard, db, &tdb, &ndb) || + !d2_ct_pi(digits, guard, block_terms, dct, &tdct, &ndct) || + !chud_direct_pi(digits, guard, cb, &tcb, &ncb) || + !chud_ct_pi(digits, guard, block_terms, cct, &tcct, &ncct)) return 0; + + int e1 = files_equal(ref, db); + int e2 = files_equal(ref, dct); + int e3 = files_equal(ref, cb); + int e4 = files_equal(ref, cct); + printf("%llu,%lu,%.9f,%.9f,%.9f,%.9f,%.9f,%.6f,%.6f,%d,%d,%d,%d\n", + (unsigned long long)digits, block_terms, + tr, tdb, tdct, tcb, tcct, + tdct / tdb, tdct / tcct, + e1, e2, e3, e4); + fflush(stdout); + return e1 && e2 && e3 && e4; +} + +int main(int argc, char **argv) { + if (argc != 3) { + fprintf(stderr, "usage: %s DIGITS BLOCK_TERMS\n", argv[0]); + return 2; + } + uint64_t digits = strtoull(argv[1], NULL, 10); + unsigned long block = strtoul(argv[2], NULL, 10); + printf("digits,block_terms,chud_bs_s,d2_direct_s,d2_ct_s,chud_direct_s,chud_ct_s,d2_ct_over_d2_direct,d2_ct_over_chud_ct,eq_d2_base,eq_d2_ct,eq_chud_base,eq_chud_ct\n"); + return (digits && block && run_ct(digits, block)) ? 0 : 1; +} diff --git a/validation/direct_convergence_taper_results.md b/validation/direct_convergence_taper_results.md new file mode 100644 index 0000000..d50b63d --- /dev/null +++ b/validation/direct_convergence_taper_results.md @@ -0,0 +1,67 @@ +# Convergence-tapered direct recurrence + +This experiment deliberately stays on the original scalar direct-recurrence path for the class-number-2 member `Delta = -427`. It does not replace the series by binary splitting and it does not carry an exact quadratic-field product tree through the summation. + +## Accelerator + +For the direct hypergeometric term recurrence, the magnitude of term `n` falls approximately as + +`|t_n| ~ 10^(-q n)`. + +For D2, + +`q = 24.9558996576542667...` decimal digits per term. + +For the D=163 / Chudnovsky direct recurrence, + +`q = 14.1816474627254777...` decimal digits per term. + +The Convergence-Tapered Direct Recurrence (CTDR) therefore reduces the relative working precision of the recurrence as the terms become smaller. In the implementation, the active bit precision is approximately + +`P_bits - q*n*log2(10) + safety_bits`, + +with a conservative safety allowance and a minimum precision floor. Precision changes are made only at coarse term-block boundaries to reduce precision-management overhead. The full-precision `F` and `T = sum(n*t_n)` accumulators are retained. + +The formula, term count, term recurrence, and final observation are otherwise unchanged. The same accelerator is applied to D2 and to the one-dimensional D=163 / Chudnovsky direct recurrence. + +## Rejected variant + +An earlier magnitude-window prototype also reduced the precision of local block accumulators. It was faster, but byte validation failed. The reason is structural: although the magnitude of a later term is small, its mantissa still contributes information across the remaining requested digits. Limiting a local accumulator to only the inter-term magnitude span discards deep digits. That variant is rejected and is not the reported CTDR method. + +## Validation + +GitHub-hosted Ubuntu 24.04 runners, GMP, GCC C17, `-O3 -march=native`. + +Every reported D2 CTDR and Chud CTDR output is byte-identical to the Chudnovsky binary-splitting reference at the requested decimal precision. + +The high-precision validation run used GitHub Actions run `32731870165`, job `97445661238`. + +| digits | taper block | Chud BS | D2 direct | D2 CTDR | Chud direct | Chud CTDR | D2 CTDR / D2 direct | D2 CTDR / Chud CTDR | +| ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | +| 100,000 | 16 | 0.024451971 s | 4.376951456 s | 1.952798843 s | 7.641611338 s | 3.396393776 s | 0.446155x | 0.574962x | +| 100,000 | 64 | 0.024074316 s | 4.366041899 s | 1.973523617 s | 7.648643017 s | 3.419784546 s | 0.452017x | 0.577090x | +| 200,000 | 16 | 0.062802076 s | 22.369212151 s | 9.841870785 s | 39.258158445 s | 17.253968477 s | 0.439974x | 0.570412x | +| 200,000 | 64 | 0.062806129 s | 22.368253231 s | 9.981220245 s | 39.240109205 s | 17.294273138 s | 0.446223x | 0.577140x | + +At 200,000 digits with a 16-term taper block: + +- CTDR reduces D2 direct time from `22.369212151 s` to `9.841870785 s`, about `2.27x` faster. +- CTDR reduces Chudnovsky direct time from `39.258158445 s` to `17.253968477 s`. +- Under the same CTDR accelerator, `D2 / Chud = 0.570412x`, so D2 is about `1.75x` faster than the equally accelerated Chudnovsky direct recurrence. + +The measured `0.5704x` ratio is close to the formula-level asymptotic term-count ratio + +`14.1816474627 / 24.9558996577 ~= 0.568`. + +This shows that the original D2 formula-level convergence advantage continues to translate almost one-for-one into wall-clock advantage when both direct recurrences receive the same convergence-aware precision accelerator. + +## Boundary + +CTDR is an accelerator for the direct recurrence state. It is not competitive with optimized Chudnovsky binary splitting at these precisions. At 200,000 digits in the same run, Chudnovsky BS took about `0.0628 s`, while D2 CTDR took about `9.84 s`. + +Therefore the result supports two narrower claims: + +1. Precision tapering is a substantial valid acceleration of the original direct D2 recurrence. +2. D2 retains approximately its full formula-level advantage over Chudnovsky when the same direct-state accelerator is applied to both. + +It does not establish superiority over binary-splitting Chudnovsky.