"""Characteristic function of X ~ g(t)^2/m2 on [0,T], g = 1/(c+bt), b = k-1, in Arb (acb):
phi(xi) = (1/m2) int_0^T e^{i xi t} (c+bt)^{-2} dt.
With u = c+bt, a = xi/b:  int_0^T e^{i xi t}(c+bt)^{-2} dt = (e^{-i a c}/b) int_{c}^{c+bT} e^{i a u} u^{-2} du,
int_{u1}^{u2} e^{iau} u^{-2} du = [-e^{iau}/u]_{u1}^{u2} + i a (E1(-i a u1) - E1(-i a u2)).
Arb's E1 is badly conditioned at some arguments (|z| ~ 10^2), so constants are recomputed at each working
precision and the precision is raised until the result is sharp.  E1 has a branch cut on (-inf,0]; for the
'analytic' flag we refuse balls whose argument may touch it."""
from flint import acb, arb, ctx

def make_phi(k, theta, beta):
    cache = {}
    def consts():
        p = ctx.prec
        if p not in cache:
            kk = arb(k); L = kk.log(); c = arb(theta) / L; T = arb(beta) / L; b = kk - 1
            m2 = (1 / b) * (1 / c - 1 / (c + b * T))
            cache[p] = (c, T, b, m2, c, c + b * T)
        return cache[p]
    def phi_raw(xi, analytic=False):
        c, T, b, m2, u1, u2 = consts()
        a = xi / b
        z1 = acb(0, -1) * a * u1
        z2 = acb(0, -1) * a * u2
        if analytic:
            for z in (z1, z2):
                re, im = z.real, z.imag
                if re.lower() <= 0 and im.lower() <= 0 <= im.upper():
                    return acb("nan")
        E1 = lambda z: z.expint(1)
        prim = lambda u: -(acb(0, 1) * a * u).exp() / u
        val = (prim(u2) - prim(u1)) + acb(0, 1) * a * (E1(z1) - E1(z2))
        return (acb(0, -1) * a * c).exp() / b * val / m2
    def phi(xi, analytic=False, tol=1e-40):
        p0 = ctx.prec
        try:
            v = None
            for prec in (p0, 2 * p0, 4 * p0, 8 * p0):
                ctx.prec = prec
                v = phi_raw(acb(xi), analytic)
                if not v.is_finite():
                    return v
                if float(v.real.rad()) < tol and float(v.imag.rad()) < tol:
                    break
            return v
        finally:
            ctx.prec = p0
    c0, T0, b0, m20, _, _ = consts()
    return phi, dict(c=c0, T=T0, b=b0, m2=m20)

if __name__ == "__main__":
    ctx.prec = 256
    phi, P = make_phi(54, 1.0, 3.0)
    for x in [10, 232.2857, 233.639, 1000]:
        v = phi(acb(arb(str(x))))
        print(x, float(v.real.rad()), float(v.abs_upper()))