Exact SymPy verification

Exact SymPy verification#

This notebook checks the displayed polynomial map over the exact rational numbers. Every decisive line is an assert: there is no floating-point arithmetic, random testing, or network access.

import sympy as sp

sp.__version__

'1.13.3'

The map and its Jacobian#

The source variables are ordered as \((x,y,z)\) and the output coordinates as \((F_1,F_2,F_3)\).

x, y, z = sp.symbols("x y z")

F1 = (1 + x*y)**3*z + y**2*(1 + x*y)*(4 + 3*x*y)
F2 = y + 3*x*(1 + x*y)**2*z + 3*x*y**2*(4 + 3*x*y)
F3 = 2*x - 3*x**2*y - x**3*z
F = sp.Matrix([F1, F2, F3])
variables = (x, y, z)

F

\[\begin{split}\displaystyle \left[\begin{matrix}y^{2} \left(x y + 1\right) \left(3 x y + 4\right) + z \left(x y + 1\right)^{3}\\3 x y^{2} \left(3 x y + 4\right) + 3 x z \left(x y + 1\right)^{2} + y\\- x^{3} z - 3 x^{2} y + 2 x\end{matrix}\right]\end{split}\]

JF = F.jacobian(variables)
det_JF = sp.factor(JF.det())
assert det_JF == -2
det_JF

\[\displaystyle -2\]

Three exact points in one fiber#

All coordinates are Integer or Rational objects.

Q = sp.Rational
points = [
    (sp.Integer(0), sp.Integer(0), -Q(1, 4)),
    (sp.Integer(1), -Q(3, 2), Q(13, 2)),
    (-sp.Integer(1), Q(3, 2), Q(13, 2)),
]
target = (-Q(1, 4), sp.Integer(0), sp.Integer(0))

images = [
    tuple(sp.expand(f.subs(dict(zip(variables, point)))) for f in F)
    for point in points
]
assert len(set(points)) == 3
assert images == [target, target, target]
list(zip(points, images))

[((0, 0, -1/4), (-1/4, 0, 0)),
 ((1, -3/2, 13/2), (-1/4, 0, 0)),
 ((-1, 3/2, 13/2), (-1/4, 0, 0))]

Determinant-one normalizations#

Scaling the first target coordinate by \(-1/2\) gives the immediate normalization \(\widetilde F=(-F_1/2,F_2,F_3)\). A second check uses an integral normalization whose differential at the origin is the identity.

F_tilde = sp.Matrix([-F1 / 2, F2, F3])
det_tilde = sp.factor(F_tilde.jacobian(variables).det())
assert det_tilde == 1
det_tilde

\[\displaystyle 1\]

X, Y, Z = sp.symbols("X Y Z")
T = 1 + 2*X*Y
H = sp.Matrix([
    X - 3*X**2*Y - 2*X**3*Z,
    Y + 6*X*T**2*Z + 6*X*Y**2*(4 + 6*X*Y),
    T**3*Z + Y**2*T*(4 + 6*X*Y),
])
JH = H.jacobian((X, Y, Z))
det_JH = sp.factor(JH.det())
linear_part = JH.subs({X: 0, Y: 0, Z: 0})
normalized_points = [
    (sp.Integer(0), sp.Integer(0), -Q(1, 4)),
    (Q(1, 2), -Q(3, 2), Q(13, 2)),
    (-Q(1, 2), Q(3, 2), Q(13, 2)),
]
normalized_target = (sp.Integer(0), sp.Integer(0), -Q(1, 4))
normalized_images = [
    tuple(sp.expand(h.subs(dict(zip((X, Y, Z), point)))) for h in H)
    for point in normalized_points
]
assert det_JH == 1
assert linear_part == sp.eye(3)
assert normalized_images == [normalized_target] * 3
(det_JH, linear_part, normalized_images)

(1,
 Matrix([
 [1, 0, 0],
 [0, 1, 0],
 [0, 0, 1]]),
 [(0, 0, -1/4), (0, 0, -1/4), (0, 0, -1/4)])

Binary-cubic identities and the compact certificate#

For abstract output coordinates \((a,b,c)\) define the binary cubic

\[Q_{a,b,c}(U,V)=cU^3-2U^2V+bUV^2-2aV^3.\]

At \((a,b,c)=F(x,y,z)\), \(U=1+xy\), and \(V=x\), the cubic and its two first derivatives satisfy three unusually rigid identities. We verify them directly, then use a second weighted-coordinate form to compress the determinant calculation.

a, b, c, Uc, Vc = sp.symbols("a b c U V")
binary_cubic = c*Uc**3 - 2*Uc**2*Vc + b*Uc*Vc**2 - 2*a*Vc**3
cubic_substitution = {
    a: F1, b: F2, c: F3,
    Uc: 1 + x*y, Vc: x,
}
binary_cubic_identities = [
    sp.factor(binary_cubic.subs(cubic_substitution)),
    sp.factor(sp.diff(binary_cubic, Uc).subs(cubic_substitution) - 2*x),
    sp.factor(sp.diff(binary_cubic, Vc).subs(cubic_substitution) + 2*(1 + x*y)),
]
assert binary_cubic_identities == [0, 0, 0]
binary_cubic_identities

[0, 0, 0]

tau = sp.symbols("tau")
affine_cubic = c*tau**3 - 2*tau**2 + b*tau - 2*a
discriminant = sp.factor(sp.discriminant(affine_cubic, tau))
discriminant_factor = (
    27*a**2*c**2 - 18*a*b*c + 16*a + b**3*c - b**2
)
assert sp.expand(discriminant + 4*discriminant_factor) == 0
discriminant

\[\displaystyle - 4 \left(27 a^{2} c^{2} - 18 a b c + 16 a + b^{3} c - b^{2}\right)\]

t, s = sp.symbols("t s")
A = t**3*s + t*(t - 1)**2*(1 + 3*t)
B = t - 1 + 3*t**2*s + 3*(t - 1)**2*(1 + 3*t)
C = 5 - 3*t - s
weighted_substitution = {t: 1 + x*y, s: x**2*z}

weighted_identities = [
    sp.expand(x**2*F1 - A.subs(weighted_substitution)),
    sp.expand(x*F2 - B.subs(weighted_substitution)),
    sp.expand(F3 - x*C.subs(weighted_substitution)),
]
assert weighted_identities == [0, 0, 0]
weighted_identities

[0, 0, 0]

For weighted target coordinates \((A/x^2,B/x,xC)\), the following \(3\times3\) determinant is the chain-rule numerator. Its exact collapse to \(-2\) is a short certificate independent in presentation from expanding the original Jacobian.

W = sp.Matrix([
    [-2*A, sp.diff(A, t), sp.diff(A, s)],
    [-B,   sp.diff(B, t), sp.diff(B, s)],
    [C,    sp.diff(C, t), sp.diff(C, s)],
])
det_W = sp.factor(W.det())
assert det_W == -2
det_W