"""Rigorous enclosure of G = 1/2 - P4(R) + P4(R') via the Plackett path integral, in Arb."""
from flint import arb, acb, ctx
import itertools, time
ctx.prec=80
PAIRS=[(0,1),(0,2),(0,3),(1,2),(1,3),(2,3)]
OTH={p:tuple(k for k in range(4) if k not in p) for p in PAIRS}
PI=arb.pi()
R33=2*arb(2).sqrt()/3
R34=2/(arb(6).root(4)+1/arb(6).root(4))
def det3(m):
    return (m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])-m[0][1]*(m[1][0]*m[2][2]-m[1][2]*m[2][0])+m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]))
def adj(M,k,l):  # (k,l) entry of adjugate = cofactor C_{lk}
    rows=[r for r in range(4) if r!=l]; cols=[c for c in range(4) if c!=k]
    sub=[[M[r][c] for c in cols] for r in rows]; return (-1)**(k+l)*det3(sub)
def asin_a(z,a):
    i=acb(0,1); return -i*(i*z+(1-z*z).sqrt(analytic=a)).log(analytic=a)
def integrand(R):
    def f(u,a):
        t=1-u*u
        M=[[ (acb(1) if p==q else acb(0))*(1-t)+t*R[p][q] for q in range(4)] for p in range(4)]
        s=acb(0)
        for (i,j) in PAIRS:
            k,l=OTH[(i,j)]; rho=t*R[i][j]
            akl=adj(M,k,l); akk=adj(M,k,k); all_=adj(M,l,l)
            pc=-akl/(akk*all_).sqrt(analytic=a)
            s+=R[i][j]/(2*PI*(1-rho*rho).sqrt(analytic=a))*(acb(1)/4+asin_a(pc,a)/(2*PI))
        return s*2*u
    return f
def orth4(R):
    v=acb.integral(integrand(R),0,1,rel_tol=arb(2)**-40)
    return acb(1)/16+v
def star(rs,th):
    s=[(1-r*r).sqrt() for r in rs]; C=[[arb(1) if i==j else None for j in range(3)] for i in range(3)]
    for (i,j),a in zip(((0,1),(0,2),(1,2)),th): C[i][j]=C[j][i]=arb(a).cos()
    R=[[arb(1)]*4 for _ in range(4)]
    for i in range(3): R[0][i+1]=R[i+1][0]=rs[i]
    for i in range(3):
        for j in range(3):
            R[i+1][j+1]=arb(1) if i==j else rs[i]*rs[j]+s[i]*s[j]*C[i][j]
    return R
def G(rs,th):
    R=star(rs,th); Rp=[row[:] for row in R]
    for i in range(1,4): Rp[0][i]=-R[0][i]; Rp[i][0]=-R[i][0]
    return (acb(1)/2-orth4(R)+orth4(Rp)).real
if __name__=="__main__":
    import numpy as np, bb
    for th in ([2.0943951,2.0943951,2.0943951-1e-6],[1.9957090050392194,2.1445051414641263,2.1429711606762405-1e-7],[1.0,1.5,2.0]):
        t=time.time(); g=G([R33,R33,R34],th); dt=time.time()-t
        print(g, " float:",bb.G(np.array([bb.r33,bb.r33,bb.r34]),np.array(th)), f"{dt:.2f}s")

# ---- split version: rigorous integral on [0, 1-eta] plus rigorous tail enclosure ----
def integrand_t(R):
    def f(t,a):
        M=[[ (acb(1) if p==q else acb(0))*(1-t)+t*R[p][q] for q in range(4)] for p in range(4)]
        s=acb(0)
        for (i,j) in PAIRS:
            k,l=OTH[(i,j)]; rho=t*R[i][j]
            akl=adj(M,k,l); akk=adj(M,k,k); all_=adj(M,l,l)
            pc=-akl/(akk*all_).sqrt(analytic=a)
            s+=R[i][j]/(2*PI*(1-rho*rho).sqrt(analytic=a))*(acb(1)/4+asin_a(pc,a)/(2*PI))
        return s
    return f
def tail(R,eta):
    # each term R_ij*phi2(t R_ij)*c(t), c in [0,1/2]; enclose integral over [1-eta,1]
    tot=arb(0)
    for (i,j) in PAIRS:
        rho=R[i][j]; ar=abs(rho)
        if ar==0: continue
        a1=ar if ar<1 else arb(1)
        I=((a1).asin()-((1-eta)*a1).asin())/a1          # int dt/sqrt(1-t^2 rho^2)
        term=R[i][j]*I/(2*PI)                            # times c in [0,1/2]
        tot+=term*arb(0.25,0.25)                         # c ranges over [0,1/2]
    return tot
def orth4_split(R,eta=arb(10)**-7,tol=arb(2)**-36):
    v=acb.integral(integrand_t(R),0,1-eta,rel_tol=tol,abs_tol=tol)
    return arb(1)/16+v.real+tail(R,eta)
def G_split(rs,th,eta=arb(10)**-7):
    R=star(rs,th); Rp=[row[:] for row in R]
    for i in range(1,4): Rp[0][i]=-R[0][i]; Rp[i][0]=-R[i][0]
    return arb(1)/2-orth4_split(R,eta)+orth4_split(Rp,eta)