From 1199404c8838fbad7203477e504979c3fe9e357f Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 21:09:27 +0800 Subject: [PATCH 01/10] Add convergence-window direct recurrence experiment --- .../benchmark_direct_convergence_window.c | 138 ++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100644 validation/benchmark_direct_convergence_window.c diff --git a/validation/benchmark_direct_convergence_window.c b/validation/benchmark_direct_convergence_window.c new file mode 100644 index 0000000..a33c246 --- /dev/null +++ b/validation/benchmark_direct_convergence_window.c @@ -0,0 +1,138 @@ +#define main original_d2_direct_benchmark_main +#include "benchmark_original_d2_basic_vs_chud.c" +#undef main + +/* Convergence-Windowed Direct Recurrence (CWDR). + * Keeps the original scalar/direct evaluator, adds tapered precision and + * short magnitude-window accumulation. No binary splitting is used here. + */ +static mp_bitcnt_t cw_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 + 192.0L; + if (keep < 320.0L) keep = 320.0L; + if (keep > (long double)max_bits) keep = (long double)max_bits; + return (mp_bitcnt_t)ceill(keep); +} + +static mp_bitcnt_t cw_window_bits(long double digits_per_term, + unsigned long block_len, + mp_bitcnt_t max_bits) { + long double span = digits_per_term * (long double)(block_len + 2UL) * + 3.32192809488736234787L + 320.0L; + if (span > (long double)max_bits) span = (long double)max_bits; + if (span < 512.0L) span = 512.0L; + return (mp_bitcnt_t)ceill(span); +} + +static void cw_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, bF, bT, nt; + mpf_init2(term, bits); mpf_init2(zwork, bits); + mpf_init2(bF, bits); mpf_init2(bT, bits); mpf_init2(nt, bits); + mpf_set(zwork, z); + 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 wbits = cw_window_bits(digits_per_term, b - a, bits); + mpf_set_prec_raw(bF, wbits); mpf_set_prec_raw(bT, wbits); + mpf_set_prec_raw(nt, wbits); + mpf_set_ui(bF, 0UL); mpf_set_ui(bT, 0UL); + + for (unsigned long n = a; n < b; ++n) { + mp_bitcnt_t abits = cw_active_bits(bits, digits_per_term, n); + mpf_set_prec_raw(term, abits); + mpf_set_prec_raw(zwork, abits); + hyper_term_advance(term, zwork, n); + mpf_add(bF, bF, term); + mpf_mul_ui(nt, term, n); + mpf_add(bT, bT, nt); + } + mpf_add(F, F, bF); + mpf_add(T, T, bT); + } + + mpf_set_prec_raw(term, bits); mpf_set_prec_raw(zwork, bits); + mpf_set_prec_raw(bF, bits); mpf_set_prec_raw(bT, bits); mpf_set_prec_raw(nt, bits); + mpf_clears(term, zwork, bF, bT, nt, NULL); +} + +static int d2_cw_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; + cw_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_cw_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; + cw_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_cw(uint64_t digits,unsigned long block_terms){ + const uint64_t guard=192ULL; char ref[256],db[256],dcw[256],cb[256],ccw[256]; + snprintf(ref,sizeof(ref),"build/cw_ref_%llu_%lu.txt",(unsigned long long)digits,block_terms); + snprintf(db,sizeof(db),"build/cw_d2_base_%llu_%lu.txt",(unsigned long long)digits,block_terms); + snprintf(dcw,sizeof(dcw),"build/cw_d2_%llu_%lu.txt",(unsigned long long)digits,block_terms); + snprintf(cb,sizeof(cb),"build/cw_chud_base_%llu_%lu.txt",(unsigned long long)digits,block_terms); + snprintf(ccw,sizeof(ccw),"build/cw_chud_%llu_%lu.txt",(unsigned long long)digits,block_terms); + double tr=0,tdb=0,tdcw=0,tcb=0,tccw=0; unsigned long nr=0,ndb=0,ndcw=0,ncb=0,nccw=0; + if(!chud_bs_pi(digits,guard,ref,&tr,&nr)||!d2_basic_pi(digits,guard,db,&tdb,&ndb)|| + !d2_cw_pi(digits,guard,block_terms,dcw,&tdcw,&ndcw)|| + !chud_direct_pi(digits,guard,cb,&tcb,&ncb)||!chud_cw_pi(digits,guard,block_terms,ccw,&tccw,&nccw))return 0; + int e1=files_equal(ref,db),e2=files_equal(ref,dcw),e3=files_equal(ref,cb),e4=files_equal(ref,ccw); + printf("%llu,%lu,%.9f,%.9f,%.9f,%.9f,%.9f,%.6f,%.6f,%d,%d,%d,%d\n", + (unsigned long long)digits,block_terms,tr,tdb,tdcw,tcb,tccw,tdcw/tdb,tdcw/tccw,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 b=strtoul(argv[2],NULL,10); + printf("digits,block_terms,chud_bs_s,d2_direct_s,d2_cw_s,chud_direct_s,chud_cw_s,d2_cw_over_d2_direct,d2_cw_over_chud_cw,eq_d2_base,eq_d2_cw,eq_chud_base,eq_chud_cw\n"); + return (digits&&b&&run_cw(digits,b))?0:1; +} From 43b6fdd863b3f7ce2278b1f9c79eb4d39a33920b Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 21:09:52 +0800 Subject: [PATCH 02/10] Remove placeholder from experiment branch --- noop | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 noop diff --git a/noop b/noop deleted file mode 100644 index e69de29..0000000 From a199e0257d058b96375d54d538e517849f2f8253 Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 21:10:10 +0800 Subject: [PATCH 03/10] Add CWDR benchmark workflow --- .../benchmark-direct-convergence-window.yml | 40 +++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 .github/workflows/benchmark-direct-convergence-window.yml diff --git a/.github/workflows/benchmark-direct-convergence-window.yml b/.github/workflows/benchmark-direct-convergence-window.yml new file mode 100644 index 0000000..0df073e --- /dev/null +++ b/.github/workflows/benchmark-direct-convergence-window.yml @@ -0,0 +1,40 @@ +name: Direct convergence-window 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 CWDR benchmark + run: | + mkdir -p build + gcc -O3 -march=native -std=c17 -Wall -Wextra -Wpedantic \ + -Isrc -Ivalidation validation/benchmark_direct_convergence_window.c \ + src/ramanujan_c_common.c -lgmp -lm \ + -o build/benchmark_direct_convergence_window + - name: Sweep convergence-window size + run: | + set -euo pipefail + rm -rf /dev/shm/rj_cw + mkdir -p /dev/shm/rj_cw/build + cp build/benchmark_direct_convergence_window /dev/shm/rj_cw/ + cd /dev/shm/rj_cw + ./benchmark_direct_convergence_window 10000 64 >/dev/null + echo 'block,digits,chud_bs_s,d2_direct_s,d2_cw_s,chud_direct_s,chud_cw_s,d2_cw_over_d2_direct,d2_cw_over_chud_cw,equal' > summary.csv + for block in 16 32 64 128 256; do + row=$(./benchmark_direct_convergence_window 100000 "$block" | tail -n 1) + echo "$block,$(echo "$row" | cut -d, -f1),$(echo "$row" | cut -d, -f3-9),$(echo "$row" | awk -F, '{print ($10&&$11&&$12&&$13)?1:0}')" >> summary.csv + done + cat summary.csv + cp summary.csv "$GITHUB_WORKSPACE/build/" + - uses: actions/upload-artifact@v4 + with: + name: direct-convergence-window + path: build/summary.csv From 5b2190d0b9b9ba8190825468e1d003c98b652359 Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 21:11:30 +0800 Subject: [PATCH 04/10] Capture CWDR correctness failures during sweep --- .github/workflows/benchmark-direct-convergence-window.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/benchmark-direct-convergence-window.yml b/.github/workflows/benchmark-direct-convergence-window.yml index 0df073e..e6932a6 100644 --- a/.github/workflows/benchmark-direct-convergence-window.yml +++ b/.github/workflows/benchmark-direct-convergence-window.yml @@ -26,10 +26,10 @@ jobs: mkdir -p /dev/shm/rj_cw/build cp build/benchmark_direct_convergence_window /dev/shm/rj_cw/ cd /dev/shm/rj_cw - ./benchmark_direct_convergence_window 10000 64 >/dev/null + ./benchmark_direct_convergence_window 10000 64 >/dev/null || true echo 'block,digits,chud_bs_s,d2_direct_s,d2_cw_s,chud_direct_s,chud_cw_s,d2_cw_over_d2_direct,d2_cw_over_chud_cw,equal' > summary.csv for block in 16 32 64 128 256; do - row=$(./benchmark_direct_convergence_window 100000 "$block" | tail -n 1) + row=$( (./benchmark_direct_convergence_window 100000 "$block" || true) | tail -n 1) echo "$block,$(echo "$row" | cut -d, -f1),$(echo "$row" | cut -d, -f3-9),$(echo "$row" | awk -F, '{print ($10&&$11&&$12&&$13)?1:0}')" >> summary.csv done cat summary.csv From c171670cbb05703b9028e821f2510aebec58afe0 Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 21:14:26 +0800 Subject: [PATCH 05/10] Expose CWDR correctness flags --- .../benchmark-direct-convergence-window.yml | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/.github/workflows/benchmark-direct-convergence-window.yml b/.github/workflows/benchmark-direct-convergence-window.yml index e6932a6..a14ba7f 100644 --- a/.github/workflows/benchmark-direct-convergence-window.yml +++ b/.github/workflows/benchmark-direct-convergence-window.yml @@ -19,18 +19,19 @@ jobs: -Isrc -Ivalidation validation/benchmark_direct_convergence_window.c \ src/ramanujan_c_common.c -lgmp -lm \ -o build/benchmark_direct_convergence_window - - name: Sweep convergence-window size + - name: Diagnose correctness run: | set -euo pipefail rm -rf /dev/shm/rj_cw mkdir -p /dev/shm/rj_cw/build cp build/benchmark_direct_convergence_window /dev/shm/rj_cw/ cd /dev/shm/rj_cw - ./benchmark_direct_convergence_window 10000 64 >/dev/null || true - echo 'block,digits,chud_bs_s,d2_direct_s,d2_cw_s,chud_direct_s,chud_cw_s,d2_cw_over_d2_direct,d2_cw_over_chud_cw,equal' > summary.csv - for block in 16 32 64 128 256; do - row=$( (./benchmark_direct_convergence_window 100000 "$block" || true) | tail -n 1) - echo "$block,$(echo "$row" | cut -d, -f1),$(echo "$row" | cut -d, -f3-9),$(echo "$row" | awk -F, '{print ($10&&$11&&$12&&$13)?1:0}')" >> summary.csv + echo 'digits,block,chud_bs_s,d2_direct_s,d2_cw_s,chud_direct_s,chud_cw_s,d2_cw_over_d2_direct,d2_cw_over_chud_cw,eq_d2_base,eq_d2_cw,eq_chud_base,eq_chud_cw' > summary.csv + for digits in 10000 30000; do + for block in 16 64 256; do + row=$( (./benchmark_direct_convergence_window "$digits" "$block" || true) | tail -n 1) + echo "$row" >> summary.csv + done done cat summary.csv cp summary.csv "$GITHUB_WORKSPACE/build/" From ec821989d6afab7e97a104c6a7ca2053a62da388 Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 21:16:12 +0800 Subject: [PATCH 06/10] Replace lossy window sums with safe tapered direct recurrence --- .../benchmark_direct_convergence_window.c | 236 +++++++++++------- 1 file changed, 144 insertions(+), 92 deletions(-) diff --git a/validation/benchmark_direct_convergence_window.c b/validation/benchmark_direct_convergence_window.c index a33c246..14ced72 100644 --- a/validation/benchmark_direct_convergence_window.c +++ b/validation/benchmark_direct_convergence_window.c @@ -2,137 +2,189 @@ #include "benchmark_original_d2_basic_vs_chud.c" #undef main -/* Convergence-Windowed Direct Recurrence (CWDR). - * Keeps the original scalar/direct evaluator, adds tapered precision and - * short magnitude-window accumulation. No binary splitting is used here. +/* 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 cw_active_bits(mp_bitcnt_t max_bits, +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 + 192.0L; - if (keep < 320.0L) keep = 320.0L; + 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 mp_bitcnt_t cw_window_bits(long double digits_per_term, - unsigned long block_len, - mp_bitcnt_t max_bits) { - long double span = digits_per_term * (long double)(block_len + 2UL) * - 3.32192809488736234787L + 320.0L; - if (span > (long double)max_bits) span = (long double)max_bits; - if (span < 512.0L) span = 512.0L; - return (mp_bitcnt_t)ceill(span); -} - -static void cw_series(const mpf_t z, +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, bF, bT, nt; - mpf_init2(term, bits); mpf_init2(zwork, bits); - mpf_init2(bF, bits); mpf_init2(bT, bits); mpf_init2(nt, bits); - mpf_set(zwork, z); - mpf_set_ui(term, 1UL); mpf_set_ui(F, 1UL); mpf_set_ui(T, 0UL); + 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 wbits = cw_window_bits(digits_per_term, b - a, bits); - mpf_set_prec_raw(bF, wbits); mpf_set_prec_raw(bT, wbits); - mpf_set_prec_raw(nt, wbits); - mpf_set_ui(bF, 0UL); mpf_set_ui(bT, 0UL); + 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) { - mp_bitcnt_t abits = cw_active_bits(bits, digits_per_term, n); - mpf_set_prec_raw(term, abits); - mpf_set_prec_raw(zwork, abits); hyper_term_advance(term, zwork, n); - mpf_add(bF, bF, term); + mpf_add(F, F, term); mpf_mul_ui(nt, term, n); - mpf_add(bT, bT, nt); + mpf_add(T, T, nt); } - mpf_add(F, F, bF); - mpf_add(T, T, bT); } - mpf_set_prec_raw(term, bits); mpf_set_prec_raw(zwork, bits); - mpf_set_prec_raw(bF, bits); mpf_set_prec_raw(bT, bits); mpf_set_prec_raw(nt, bits); - mpf_clears(term, zwork, bF, bT, nt, NULL); + mpf_clears(term, zwork, nt, NULL); } -static int d2_cw_pi(uint64_t digits, uint64_t guard, unsigned long block_terms, +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); + 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; - cw_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); + 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_cw_pi(uint64_t digits, uint64_t guard, unsigned long block_terms, +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; - cw_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); + 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_cw(uint64_t digits,unsigned long block_terms){ - const uint64_t guard=192ULL; char ref[256],db[256],dcw[256],cb[256],ccw[256]; - snprintf(ref,sizeof(ref),"build/cw_ref_%llu_%lu.txt",(unsigned long long)digits,block_terms); - snprintf(db,sizeof(db),"build/cw_d2_base_%llu_%lu.txt",(unsigned long long)digits,block_terms); - snprintf(dcw,sizeof(dcw),"build/cw_d2_%llu_%lu.txt",(unsigned long long)digits,block_terms); - snprintf(cb,sizeof(cb),"build/cw_chud_base_%llu_%lu.txt",(unsigned long long)digits,block_terms); - snprintf(ccw,sizeof(ccw),"build/cw_chud_%llu_%lu.txt",(unsigned long long)digits,block_terms); - double tr=0,tdb=0,tdcw=0,tcb=0,tccw=0; unsigned long nr=0,ndb=0,ndcw=0,ncb=0,nccw=0; - if(!chud_bs_pi(digits,guard,ref,&tr,&nr)||!d2_basic_pi(digits,guard,db,&tdb,&ndb)|| - !d2_cw_pi(digits,guard,block_terms,dcw,&tdcw,&ndcw)|| - !chud_direct_pi(digits,guard,cb,&tcb,&ncb)||!chud_cw_pi(digits,guard,block_terms,ccw,&tccw,&nccw))return 0; - int e1=files_equal(ref,db),e2=files_equal(ref,dcw),e3=files_equal(ref,cb),e4=files_equal(ref,ccw); +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,tdcw,tcb,tccw,tdcw/tdb,tdcw/tccw,e1,e2,e3,e4); - fflush(stdout); return e1&&e2&&e3&&e4; + (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 b=strtoul(argv[2],NULL,10); - printf("digits,block_terms,chud_bs_s,d2_direct_s,d2_cw_s,chud_direct_s,chud_cw_s,d2_cw_over_d2_direct,d2_cw_over_chud_cw,eq_d2_base,eq_d2_cw,eq_chud_base,eq_chud_cw\n"); - return (digits&&b&&run_cw(digits,b))?0:1; +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; } From 98a8e393cbd60b026e9c9523a9cc6dea8bcf225f Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 21:17:53 +0800 Subject: [PATCH 07/10] Run high-precision CTDR validation --- .../benchmark-direct-convergence-window.yml | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/.github/workflows/benchmark-direct-convergence-window.yml b/.github/workflows/benchmark-direct-convergence-window.yml index a14ba7f..a84a018 100644 --- a/.github/workflows/benchmark-direct-convergence-window.yml +++ b/.github/workflows/benchmark-direct-convergence-window.yml @@ -12,28 +12,28 @@ jobs: - uses: actions/checkout@v4 - name: Install dependencies run: sudo apt-get update && sudo apt-get install -y build-essential libgmp-dev - - name: Build CWDR benchmark + - name: Build CTDR benchmark run: | mkdir -p build gcc -O3 -march=native -std=c17 -Wall -Wextra -Wpedantic \ -Isrc -Ivalidation validation/benchmark_direct_convergence_window.c \ src/ramanujan_c_common.c -lgmp -lm \ -o build/benchmark_direct_convergence_window - - name: Diagnose correctness + - name: High precision validation run: | set -euo pipefail - rm -rf /dev/shm/rj_cw - mkdir -p /dev/shm/rj_cw/build - cp build/benchmark_direct_convergence_window /dev/shm/rj_cw/ - cd /dev/shm/rj_cw - echo 'digits,block,chud_bs_s,d2_direct_s,d2_cw_s,chud_direct_s,chud_cw_s,d2_cw_over_d2_direct,d2_cw_over_chud_cw,eq_d2_base,eq_d2_cw,eq_chud_base,eq_chud_cw' > summary.csv - for digits in 10000 30000; do - for block in 16 64 256; do - row=$( (./benchmark_direct_convergence_window "$digits" "$block" || true) | tail -n 1) - echo "$row" >> summary.csv - done + rm -rf /dev/shm/rj_ct + mkdir -p /dev/shm/rj_ct/build + cp build/benchmark_direct_convergence_window /dev/shm/rj_ct/ + cd /dev/shm/rj_ct + echo 'digits,block,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' > summary.csv + for spec in '100000 16' '100000 64' '200000 16' '200000 64'; do + set -- $spec + row=$(./benchmark_direct_convergence_window "$1" "$2" | tail -n 1) + echo "$row" >> summary.csv done cat summary.csv + awk -F, 'NR>1 { if (!($10 && $11 && $12 && $13)) exit 1 }' summary.csv cp summary.csv "$GITHUB_WORKSPACE/build/" - uses: actions/upload-artifact@v4 with: From 26e90a4434c0f2c30924fb7472790e0d4ddd6122 Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 21:24:01 +0800 Subject: [PATCH 08/10] Add convergence-tapered direct benchmark target --- Makefile | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) 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}'); \ From e8b05d94c0cd99d35d84d982b2daaa04ed0c4332 Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 21:24:28 +0800 Subject: [PATCH 09/10] Document convergence-tapered direct recurrence results --- .../direct_convergence_taper_results.md | 67 +++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 validation/direct_convergence_taper_results.md 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. From 81f413389883dca9c504db3af428d7f3ad5a8865 Mon Sep 17 00:00:00 2001 From: Kestis Date: Mon, 24 Aug 2026 21:24:49 +0800 Subject: [PATCH 10/10] Remove temporary CTDR benchmark workflow --- .../benchmark-direct-convergence-window.yml | 41 ------------------- 1 file changed, 41 deletions(-) delete mode 100644 .github/workflows/benchmark-direct-convergence-window.yml diff --git a/.github/workflows/benchmark-direct-convergence-window.yml b/.github/workflows/benchmark-direct-convergence-window.yml deleted file mode 100644 index a84a018..0000000 --- a/.github/workflows/benchmark-direct-convergence-window.yml +++ /dev/null @@ -1,41 +0,0 @@ -name: Direct convergence-window 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 CTDR benchmark - run: | - mkdir -p build - gcc -O3 -march=native -std=c17 -Wall -Wextra -Wpedantic \ - -Isrc -Ivalidation validation/benchmark_direct_convergence_window.c \ - src/ramanujan_c_common.c -lgmp -lm \ - -o build/benchmark_direct_convergence_window - - name: High precision validation - run: | - set -euo pipefail - rm -rf /dev/shm/rj_ct - mkdir -p /dev/shm/rj_ct/build - cp build/benchmark_direct_convergence_window /dev/shm/rj_ct/ - cd /dev/shm/rj_ct - echo 'digits,block,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' > summary.csv - for spec in '100000 16' '100000 64' '200000 16' '200000 64'; do - set -- $spec - row=$(./benchmark_direct_convergence_window "$1" "$2" | tail -n 1) - echo "$row" >> summary.csv - done - cat summary.csv - awk -F, 'NR>1 { if (!($10 && $11 && $12 && $13)) exit 1 }' summary.csv - cp summary.csv "$GITHUB_WORKSPACE/build/" - - uses: actions/upload-artifact@v4 - with: - name: direct-convergence-window - path: build/summary.csv