"""eval_closed_form_offslice.py -- closed form of V(X;P), conformally coupled one-loop triangle, general kinematics.
   V = (sqrt2 pi/8) * Re sum_s c_s LS_s(X,P) F_s(X,P),  s = 4 vertex classes + 6 tube radicands,
   F_s = combinations of D(x) = Li2(x) - Li2(1/x) (vertex classes) and D(x,xbar) = Li2(x) - Li2(xbar) (tubes) read from li2_representation.json.
Branch prescription: principal Li2 of mpmath; the REAL PART of LS_s F_s is taken (validated numerically on the region X_v > P_v; see note).
Sector signs sigma_s (orientation conventions of the symbol normalisation) are stored in li2_representation.json['sigma'].
Region: X_v >= P_v > 0 with P a triangle. Outside it the evaluator raises OutOfRegion (for X_v < P_v use eval_closed_form_V_all.py).
usage: python eval_closed_form_offslice.py X1 X2 X3 P1 P2 P3 [dps]"""
import os, sys, json, itertools
import mpmath as mp
_here = os.path.dirname(os.path.abspath(__file__))
SPEC = json.load(open(os.path.join(_here, 'li2_representation.json')))
SECTORS = ['S0', 'S12', 'S23', 'S31', 't12m', 't12p', 't13m', 't13p', 't23m', 't23p']
COEF = {'S0': mp.mpf(-1)/4, 'S12': mp.mpf(-1)/2, 'S23': mp.mpf(-1)/2, 'S31': mp.mpf(-1)/2, 't12m': -1, 't12p': 1, 't13m': -1, 't13p': 1, 't23m': -1, 't23p': -1}
VERT = {'S0': ((1, 0, 1, -1, 0, 0), (1, 1, 0, 0, -1, 0), (0, 1, 1, 0, 0, -1)), 'S12': None}
_TUBE = {k: [compile(e, 'p', 'eval') for e in SPEC['tube']['ps'][k]] for k in 'pm'}
def _fr(s):
    n, d = (str(s).split('/') + ['1'])[:2]; return mp.mpf(int(n))/int(d)
def _det(M):
    return mp.det(mp.matrix(M))
def _vertex_y(X, cl):
    X1, X2, X3 = X; Q = X1 + X2 + X3
    if cl == 'S0': return [(X3 - X1 - X2)/2, (X1 - X2 - X3)/2, (X2 - X1 - X3)/2]
    if cl == 'S12': return [-Q/2, -X2 + Q/2, -X1 + Q/2]          # planes g1, g2, G12: a = -Q/2, b = -X2 - a, c = -X1 - a
    if cl == 'S23': return [-X2 + Q/2, -Q/2, -X3 + Q/2]          # g2, g3, G23: b = -Q/2, a = -X2 - b, c = -X3 - b
    if cl == 'S31': return [-X1 + Q/2, -X3 + Q/2, -Q/2]          # g1, g3, G31: c = -Q/2, a = -X1 - c, b = -X3 - c
def F_vertex(X, P, cl):
    y = _vertex_y(X, cl); P1, P2, P3 = P; u = [t*t for t in y]
    M = [[0, 1, 1, 1, 1], [1, 0, P2*P2, P1*P1, u[0]], [1, P2*P2, 0, P3*P3, u[1]], [1, P1*P1, P3*P3, 0, u[2]], [1, u[0], u[1], u[2], 0]]
    rad = -2*_det(M); w = mp.sqrt(mp.mpc(rad))
    dist = {(1, 2): P2, (1, 3): P1, (2, 3): P3, (1, 4): y[0], (2, 4): y[1], (3, 4): y[2]}
    def cof(i, j):
        return _det([[M[r][c] for c in range(5) if c != j] for r in range(5) if r != i])*(-1)**(i + j)
    C = {}; ell = {}
    for (i, j) in itertools.combinations(range(1, 5), 2):
        rem = tuple(k for k in range(1, 5) if k not in (i, j)); C[(i, j)] = cof(i, j); ell[(i, j)] = dist[rem]
    Cd = {i: cof(i, i) for i in range(1, 5)}
    tot = mp.mpc(0)
    for c, cyc, eps in SPEC['vertex'][cl]:
        x = mp.mpc(1)
        for e, s in zip(cyc, eps): x *= C[tuple(e)] + s*ell[tuple(e)]*w
        vs = cyc[0] if (len(cyc) == 2 and cyc[0] == cyc[1]) else sorted({i for e in cyc for i in e})
        for i in vs: x /= Cd[i]
        tot += _fr(c)*(mp.polylog(2, x) - mp.polylog(2, 1/x))
    return tot, (X[0] + X[1] + X[2])/w
def F_tube(X, P, i, j, typ):
    k = 3 - i - j; env = {'X1': X[i], 'X2': X[j], 'X3': X[k], 'P1': P[i], 'P2': P[j], 'P3': P[k]}
    sg = 1 if typ == 'p' else -1; m = P[i]**2 + P[j]**2 - P[k]**2
    T = m*m - 4*P[i]**2*P[j]**2 + 4*P[i]**2*X[j]**2 + 4*P[j]**2*X[i]**2 + 4*sg*X[i]*X[j]*m; w = mp.sqrt(mp.mpc(T))
    ps = [eval(cd, {}, env) for cd in _TUBE[typ]]; tot = mp.mpc(0)
    for a, b, c in SPEC['tube'][typ]['D']:
        tot += _fr(c)*(mp.polylog(2, (ps[a] - w)/(ps[b] - w)) - mp.polylog(2, (ps[a] + w)/(ps[b] + w)))
    for a, R, c in SPEC['tube'][typ]['LL']:
        tot += _fr(c)*mp.log(abs((ps[a] + w)/(ps[a] - w)))*mp.log(abs(eval(R, {}, env)))
    Q = X[0] + X[1] + X[2]
    LS = X[k]*Q/((X[i] - X[j] - X[k])*(X[i] - X[j] + X[k]))/w if typ == 'm' else X[k]/(X[i] + X[j] - X[k])/w
    return tot, LS
class OutOfRegion(ValueError):
    """The input lies outside the region in which this evaluator is valid."""
def _g(v):
    if isinstance(v, (int, str)) or hasattr(v, 'numerator'):
        from fractions import Fraction as _F
        f = _F(v); return mp.mpf(f.numerator)/f.denominator
    return mp.mpf(v)
def check_region(X, P):
    X = [_g(x) for x in X]; P = [_g(p) for p in P]
    if len(X) != 3 or len(P) != 3 or min(P) <= 0: raise OutOfRegion('need three site energies and three positive momenta')
    if 2*max(P) > sum(P): raise OutOfRegion('P = %s violates the triangle inequality' % [mp.nstr(p, 8) for p in P])
    low = [v + 1 for v in range(3) if X[v] < P[v]]
    if low: raise OutOfRegion('X_v < P_v at site(s) %s: this evaluator is valid for X_v >= P_v only; use eval_closed_form_V_all.py' % low)
def sector_terms(X, P):
    check_region(X, P)
    X = [mp.mpf(x) for x in X]; P = [mp.mpf(p) for p in P]; out = {}
    for s in SECTORS:
        F, LS = F_vertex(X, P, s) if s.startswith('S') else F_tube(X, P, int(s[1]) - 1, int(s[2]) - 1, s[3])
        out[s] = mp.re(LS*F)
    return out
def V_offslice(X, P):
    t = sector_terms(X, P); sig = SPEC.get('sigma', {s: 1 for s in SECTORS})
    return mp.sqrt(2)*mp.pi/8*sum(COEF[s]*sig[s]*t[s] for s in SECTORS)
if __name__ == '__main__':
    mp.mp.dps = int(sys.argv[7]) if len(sys.argv) > 7 else 40
    from fractions import Fraction as Fr
    q = lambda s: mp.mpf(Fr(s).numerator)/Fr(s).denominator
    print(mp.nstr(V_offslice([q(s) for s in sys.argv[1:4]], [q(s) for s in sys.argv[4:7]]), mp.mp.dps - 5))