"""Exact certificates for the integral-point application in Lesson 22.

All intervals use integer outward rounding. No floating point is used.
The example concerns the subgroup generated by (3,5) on y^2=x^3-2.
"""
from fractions import Fraction as F
from math import isqrt
from pathlib import Path
import json

ROOT = Path(__file__).resolve().parent
BITS = 1024
S = 1 << BITS
TERMS = 12000
checks = 0

def check(value):
    global checks
    assert value
    checks += 1

def ceildiv(a, b):
    return -((-a) // b)

def mul(a, b):
    return (a[0]*b[0]//S, ceildiv(a[1]*b[1], S))

def integral(x):
    # I(x)=x^(-1/2) sum c_k (2/x^3)^k/(6k+1), x>cuberoot(2).
    q = F(2)/x**3
    check(0 < q < 1)
    lo = hi = S
    sumlo = sumhi = 0
    for k in range(TERMS):
        sumlo += lo//(6*k+1)
        sumhi += ceildiv(hi, 6*k+1)
        a = (2*k+1)*q.numerator
        b = (2*k+2)*q.denominator
        lo = lo*a//b
        hi = ceildiv(hi*a, b)
    # For k>=TERMS the ratio of consecutive c_k q^k is <q.
    tail = ceildiv(hi*q.denominator,
                   (6*TERMS+1)*(q.denominator-q.numerator))
    z = isqrt((x.denominator*S*S)//x.numerator)
    check(z*z*x.numerator <= x.denominator*S*S < (z+1)**2*x.numerator)
    ans = mul((sumlo, sumhi+tail), (z, z+1))
    check(ans[1]-ans[0] < 1 << 32)
    return ans

def dot(a,b):
    return sum(x*y for x,y in zip(a,b))

def gauss(floor_theta, C):
    basis = [[1,floor_theta],[0,C]]
    change = [[1,0],[0,1]]
    iterations = 0
    while True:
        if dot(basis[1],basis[1]) < dot(basis[0],basis[0]):
            basis.reverse();change.reverse()
        den=dot(basis[0],basis[0]);num=dot(basis[0],basis[1])
        m=(2*num+den)//(2*den)
        if m==0:break
        basis[1]=[b-m*a for a,b in zip(basis[0],basis[1])]
        change[1]=[b-m*a for a,b in zip(change[0],change[1])]
        iterations += 1
    check(abs(change[0][0]*change[1][1]-change[0][1]*change[1][0])==1)
    check(basis==[[u,u*floor_theta+v*C] for u,v in change])
    d1=F(dot(basis[0],basis[0]))
    d2=F(dot(basis[1],basis[1]))-F(dot(basis[0],basis[1]))**2/d1
    check(d1>0 and d2>0)
    check(abs(F(dot(basis[0],basis[1]))/d1)<=F(1,2))
    check(d2 >= (F(3,4)-(F(dot(basis[0],basis[1]))/d1)**2)*d1)
    return basis,change,min(d1,d2),iterations

def add(P,Q):
    if P is None:return Q
    if Q is None:return P
    x,y=P;u,v=Q
    if x==u and y==-v:return None
    m=(3*x*x)/(2*y) if P==Q else (v-y)/(u-x)
    xx=m*m-x-u
    return xx,m*(x-xx)-y

def point_json(P):
    if P is None:return 'O'
    return [str(x) for x in P]

def main():
    I1=integral(F(3));I2=integral(F(129,100))
    omega=(2*I1[0]+I2[0],2*I1[1]+I2[1])
    check(S < omega[0] < omega[1] < 3*S)
    theta=(I2[0]*S//omega[1],ceildiv(I2[1]*S,omega[0]))
    check(0<theta[0]<theta[1]<S//2)
    check(129**3 > 912*3**7)  # q(P)>7/8, since e<3.
    check(F(8,3)**20 > 129**3*17)  # q(2P)<10.
    check(F(3*4*9*729*5,4)<100000)  # period-area budget.
    check(F(2**2*2*3)<25)  # 2 sqrt(2e)<5, since e<3.
    g=2;V=100000;b0=100000;h=3;k=4
    Cstar=3*10**6*10**(6*g)*4**(2*g*g)*(g+1)**(2*g*g+9*g+13)
    H=Cstar*V**g*(b0+1)*(b0+h+1)**(g+1)+2
    N0=ceildiv(24*H*k**k,7)
    N=N0;stages=[]
    for _ in range(5):
        C=10**(2*len(str(N))+2)
        flo=theta[0]*C//S;fhi=theta[1]*C//S
        check(flo==fhi)
        basis,change,l2,it=gauss(flo,C)
        check(l2>2*N*N)
        rad=l2-N*N
        root=isqrt(rad.numerator//rad.denominator)
        check(F(root*root)<=rad<F((root+1)**2))
        delta=root-N
        check(delta>0)
        ratio=F(5*C,delta)
        D=0
        while 10**D < ratio:D+=1
        check(D>0 and ratio<=10**D)
        new=isqrt((72*D)//7)
        stages.append({'old_N':N,'C':C,'floor_C_theta':flo,'basis_rows':basis,
                       'unimodular_rows':change,'min_gram_schmidt_squared':str(l2),
                       'delta_lower':delta,'decimal_exponent_upper':D,
                       'new_N':new,'gauss_iterations':it})
        if new>=N:break
        N=new
    P=(F(3),F(5));Q=None;positive=[];integral_points=[]
    for n in range(N+1):
        if n:Q=add(Q,P)
        if Q is not None:
            x,y=Q;check(y*y==x*x*x-2)
            positive.append({'n':n,'point':point_json(Q),
                             'integral':x.denominator==y.denominator==1})
            if x.denominator==y.denominator==1:
                integral_points.extend([[n,int(x),int(y)],[-n,int(x),-int(y)]])
    # H=3: all remaining integral x lie in [-3,2].
    small=[]
    for x in range(-3,3):
        f=x**3-2
        if f>=0 and isqrt(f)**2==f:small.append([x,isqrt(f)])
        check(f<0 or isqrt(f)**2!=f)
    check(small==[])
    check(integral_points==[[1,3,5],[-1,3,-5]])
    data={'state':'passed','scope':'Integral points in Z(3,5) on y^2=x^3-2; no full Mordell-Weil basis claim',
          'integer_interval_bits':BITS,'series_terms':TERMS,
          'I_3_interval_numerators':I1,'I_129_100_interval_numerators':I2,
          'omega_interval_numerators':omega,'theta_interval_numerators':theta,
          'common_interval_denominator':S,'initial_coefficient_bound':N0,
          'reductions':stages,'final_coefficient_bound':N,
          'positive_multiples':positive,'small_x_integral_points':small,
          'integral_points_with_coefficients':integral_points,'exact_checks':checks}
    (ROOT/'elliptic-integral-points-check.json').write_text(json.dumps(data,indent=2)+'\n',encoding='utf-8')
    print(json.dumps({'state':'passed','initial_bound_digits':len(str(N0)),
                      'bounds':[N0]+[s['new_N'] for s in stages],
                      'integral_points':integral_points,'exact_checks':checks}))

if __name__=='__main__':main()