from fractions import Fraction as F
from math import factorial
# Polynomial coefficients in ascending order; exact rational arithmetic.
R=[F(x) for x in [1,0,2,0,1,0,2]]
def derivative(a,n):
    for _ in range(n):a=[(i+1)*a[i+1] for i in range(len(a)-1)]
    return a
S=[F(0)]*7
for j in range(4):
    for i,v in enumerate(derivative(R[:],2*j)):
        S[i]+=F((-1)**j,8**j*factorial(j))*v
assert S==[F(7,32),0,F(49,8),0,F(-13,2),0,F(2)]
assert F(2)*F(13,8)**2+F(27,32)==F(49,8)
print('PASS: exact S_4 polynomial and positive sum-of-squares decomposition.')
normalization=F(1,5)*(4-F(1,2)+3*F(1,2))
assert normalization==1
assert F(2)*F(15,64)-F(13,2)*F(3,16)+F(49,8)*F(1,4)+F(7,32)==1
print('PASS: exact normalization of the stated Y density and Gaussian-polynomial U density.')
# In s=t^2, multiply both sides of the partial fraction identity by the denominator.
assert F(4,5)+F(1,5)==1
assert F(4,5)-F(4,5)==0
print('PASS: 1/[(1+2t^2)(1+t^4)] partial fractions, all coefficients exact.')
assert 2*(F(4,2)+2)==8
print('PASS: explicit-pair variance from the second-order transform coefficient is 8.')
print('Scope: finite exact identities do not certify global parameter regions, external KPS inputs, or novelty.')