"""Checks accompanying the analytic proofs. Python 3, NumPy and SciPy.

Finite computations supplement the proofs. The normalization audit uses exact
Jacobi endpoint identities and an independently evaluated terminating series.
"""
from fractions import Fraction as F
from math import comb


def rising(x, k):
    z = F(1)
    for i in range(k):
        z *= x+i
    return z


def Q(d,k):
    return (rising(2,k)*rising(d,k)*rising(d+1,k) /
            (rising(d+2,k)*rising(F(d+1,2),k)*rising(F(d+2,2),k)))


def step(d,x):
    return F(d+2,d)*(x+d)*(x+d+1)/(x+d+2)/(2*x+d+1)


def transform(d,x):
    return F(d+2,2*d)+F(d+2,4*d*(d+3))*(
        F(d*d-1)/(x+F(d+1,2))-F(8)/(x+d+2))


def jacobi(n,a,b,z):
    return sum(comb(n+a,k)*comb(n+b,n-k)*
               ((z-1)/2)**(n-k)*((z+1)/2)**k for k in range(n+1))


def cap_moments():
    for d in range(4,31):
        r=F(d*(d+1)*(d-3),16*(d-1))
        assert r*step(d,-2)==F((d+1)*(d+2)*(d-2),16*d)
        for k in range(-2,25):
            assert step(d,k)==transform(d,k)
            if k>=0:
                assert Q(d,k)*step(d,k)==Q(d+1,k)
    print('Positive-factor Mellin identity: d = 4..30, k = -2..24 passed.')


def cancellation():
    for d in range(4,17):
        for j in range(1,13):
            b=F(1)
            mixed=b
            z=F(2,5)
            product=b
            for k in range(j):
                b*=(F(k-j)*(k+j+d+1)*(k+F(d+1,2))*(k+F(d+2,2))/
                    ((k+1)*(k+d)*(k+2)*(k+d+1)))
                product+=b*(1-z*z)**(k+1)
                mixed+=b*Q(d,k+1)
            p=jacobi(j,d-1,1,z)*jacobi(j,d-1,1,-z)
            assert p==product*comb(j+d-1,j)*(-1)**j*(j+1)
            assert mixed==0
    print('Direct Jacobi products and mixed cancellation: d = 4..16, j = 1..12 passed.')


def endpoint_sum():
    for a in range(3,19):
        for n in (1,2,3,8,16,255,256,1025):
            weight=sum(F(2)**(2*j*(a-2)) for j in range(n.bit_length()) if 2**j<n)
            assert weight/n**(2*(a-2))<=F(1)/(1-F(2)**(-2*(a-2)))
    print('Endpoint geometric sum including dyadic boundaries: passed.')


def degree_one_check():
    for d in range(4,31):
        for t in (F(0),F(1,101),F(1,2),F(1)):
            value=jacobi(1,d-2,0,2*t-1)/(d-1)
            assert value==(d*t-1)/(d-1)
            if t:
                assert (value+F(1,d-1))/t==F(d,d-1)
    print('Degree-one fidelity multiplier and quotient: exact checks passed.')


def multiplier_checks():
    import numpy as np
    from scipy.special import eval_jacobi, gammaln
    from math import factorial
    grid=np.r_[np.geomspace(1e-7,0.01,50),np.linspace(0.01,0.5,60)]
    highest_ratio=0.0
    for d in range(4,13):
        a=d-2
        constant=4*2**a*factorial(a)
        for n in range(1,401):
            norm=np.exp(gammaln(n+a+1)-gammaln(n+1)-gammaln(a+1))
            quot=(eval_jacobi(n,a,0,2*grid-1)-(-1)**n)/(norm*grid)
            slope=-n*(n+d-1)*(-1)**n/norm
            maximum=max(abs(slope),np.max(abs(quot)))
            ratio=maximum*n**(a-2)/constant
            highest_ratio=max(highest_ratio,ratio)
            assert np.isfinite(ratio) and ratio<=1+1e-10
    print('Multiplier grid: d = 4..12, n = 1..400, 110 fidelities plus endpoint passed.')
    print(f'Largest ratio to the proved bound on this grid: {highest_ratio:.12f}')


def normalization_checks():
    """Exact endpoints and reproducing-kernel dimension, beyond degree one."""
    for d in range(4, 13):
        a = d-2
        for n in range(25):
            norm = comb(n+a, a)
            assert jacobi(n, a, 0, F(1)) == norm
            assert jacobi(n, a, 0, F(-1)) == (-1)**n
            hn = comb(n+d-1, d-1)**2 - comb(n+d-2, d-1)**2
            assert hn == F(2*n+a+1, a+1)*norm**2
            assert F(jacobi(n, a, 0, F(1)), norm) == 1
    for n in range(1, 401):
        slope = F(n*(n+3), comb(n+2, 2))
        assert slope == 2-F(4, (n+1)*(n+2)) < 2
    print('Normalization: exact m_1(n)=1 and kernel dimensions, d=4..12, n=0..24 passed.')
    print('d=4 endpoint: exact |m_0_prime(n)|=2-4/((n+1)(n+2)), n=1..400 passed.')
    n, a = 2, 3
    value = jacobi(n, a, 0, F(1))
    correct, incorrect = F(value, comb(n+a, a)), F(value, n+a)
    assert (correct, incorrect) == (1, 2)
    print(f'Wrong denominator diagnostic, d=5,n=2: binomial gives {correct}; n+a gives {incorrect}.')


def series_quotient(d, n, t):
    """Reflected Jacobi 2F1(-n,n+a+1;1;t), evaluated independently.

    P_n^(a,0)(2t-1)=(-1)^n sum_k (-n)_k(n+a+1)_k t^k/(k!)^2.
    The constant is excluded before summation to avoid subtractive loss as t->0.
    """
    from decimal import Decimal
    a = d-2
    term = Decimal(1)
    terms = []
    for k in range(n):
        term *= Decimal((-n+k)*(n+a+1+k))*t/(k+1)**2
        terms.append(term)
    return (-1)**n*sum(terms, Decimal(0))/(t*comb(n+a, a))


def bessel_j0_two():
    """Compute J_0(2) by its absolutely convergent defining series."""
    from decimal import Decimal
    term = Decimal(1)
    total = term
    for k in range(1, 100):
        term *= -Decimal(1)/(k*k)
        total += term
    return total


def diagonal_scaling_checks():
    """80-digit computations on the proposed t=1/n^2 sequence."""
    from decimal import Decimal, localcontext
    from scipy.special import eval_jacobi
    with localcontext() as ctx:
        ctx.prec = 80
        constant = 6*(1-bessel_j0_two())
        print('t=1/n^2: independent Jacobi series; d=5 target n*|q| =',
              f'{constant:.13f}')
        print('     n        |q| d=4        |q| d=5        n*|q| d=5')
        for n in (16, 64, 256, 1024, 4096):
            t = Decimal(1)/n**2
            q4 = abs(series_quotient(4, n, t))
            q5 = abs(series_quotient(5, n, t))
            for d, q in ((4, q4), (5, q5)):
                qs = abs((eval_jacobi(n, d-2, 0, 2*float(t)-1)-(-1)**n)
                         / (float(t)*comb(n+d-2, d-2)))
                assert abs(qs/float(q)-1) < 2e-8
            print(f'{n:6d} {float(q4):16.10e} {float(q5):16.10e} {float(n*q5):16.10e}')
        print('Independent series agrees with SciPy: relative error < 2e-8.')


def harmonic_endpoint_test():
    """Evaluate BOTH sides up to the fixed analytic constant C_5.

    f_n([z])=(2*z_1*conj(z_2))^n, ||z||=1, lies in H_n and has sup norm 1.
    Dirichlet(1,...,1) integration gives
      ||f_n||_r^r = 2^(r*n) Gamma(d) Gamma(r*n/2+1)^2/Gamma(r*n+d).
    We print LHS, the RHS without C_5, and their ratio. C_5 is not replaced
    by a fabricated numerical filter norm. An upper bound for this ratio is
    |q|^3 for EVERY nonzero harmonic in H_n, since ||f||_3^3<=||f||_inf||f||_2^2.
    """
    from decimal import Decimal, localcontext
    from math import factorial
    with localcontext() as ctx:
        ctx.prec = 80
        d = 5
        print('d=5, f_n=(2*z1*conj(z2))^n, ||f_n||_inf=1, t=1/n^2:')
        print('     n         LHS=||T f||_3^3        RHS/C_5=||f||_2^2       LHS/(RHS/C_5)')
        for n in (16, 64, 256, 1024, 4096):
            t = Decimal(1)/n**2
            q = abs(series_quotient(d, n, t))
            def moment(r):
                assert r*n % 2 == 0
                return (Decimal(2)**(r*n)*factorial(d-1)
                        *Decimal(factorial(r*n//2))**2
                        /Decimal(factorial(r*n+d-1)))
            second = moment(2)
            third = moment(3)
            lhs = q**3*third
            ratio = lhs/second
            assert 0 < third <= second and ratio <= q**3
            print(f'{n:6d} {float(lhs):24.12e} {float(second):24.12e} {float(ratio):24.12e}')
        limit = (6*(1-bessel_j0_two()))**3*(Decimal(2)/3)**Decimal('3.5')
        print('Predicted limit n^3*LHS/(RHS/C_5) =', f'{limit:.12f}')
        print('These tests show decay, not the proposed numerical divergence.')


if __name__=='__main__':
    if not __debug__:
        raise SystemExit("These checks require assertions. Run Python without -O.")
    cap_moments()
    cancellation()
    endpoint_sum()
    degree_one_check()
    multiplier_checks()
    normalization_checks()
    diagonal_scaling_checks()
    harmonic_endpoint_test()