from sage.all import PolynomialRing, ZZ
n = 5
R = PolynomialRing(ZZ, ["a%s" % i for i in range(n + 1)])
a = R.gens()
S = PolynomialRing(R, "x")
x = S.gen()
f = sum(a[i] * x**i for i in range(n + 1))
sign = 1 if (n * (n - 1) // 2) % 2 == 0 else -1
numerator = sign * f.resultant(f.derivative())
disc, rem = numerator.quo_rem(a[n])
assert rem == 0
discThe generator computes exact integer coefficients in $\mathbb Z[a_0,\dots,a_n]$ from the resultant identity in Formula (1) and all arithmetic is exact.
It checks Formula (2) by substituting the elementary symmetric polynomials of the roots, $a_i=a_n(-1)^{n-i}e_{n-i}(\alpha_1,\dots,\alpha_n)$, for every stored degree, compares the quadratic and cubic rows from Formula (3) and Formula (4) against the stored values, and compares integer specialisations with Sage's exact univariate discriminant over $\mathbb Z[x]$.