"""fF_enclose.py: rigorous Arb enclosures of 1 - f(6) and F(7) - 1 for the linear-sieve functions f, F
((sF)' = f(s-1), (sf)' = F(s-1), F = 2e^gamma/s, f = 0 on (0, 2], as in Johnston-Starichkova, arXiv:2208.01229v3, Section 2), the two inputs
that their (Fbound) and (fbound) state as F(7) - 1 <= 5e-6 and 1 - f(6) <= 1.05e-4 and that Appendix B of the paper uses.

Closed forms quoted by Johnston-Starichkova (after Cai, Lemma 2):
  f(s) = (2e^g/s) (log(s-1) + int_3^{s-1} (1/t) int_2^{t-1} log(u-1)/u du dt),                       4 <= s <= 6;
  F(s) = (2e^g/s) (1 + int_2^{s-1} log(t-1)/t dt + int_2^{s-3} (log(t-1)/t) int_{t+2}^{s-1} (1/u) log((u-1)/(t+1)) du dt),  5 <= s <= 7.
Inner integrals in closed form: int log(u-1)/u du = Phi(u) := (log u)^2/2 + Li_2(1/u) (d/du Phi = (log u + log(1 - 1/u))/u).
Outer integrals by Arb's rigorous integrator (acb.integral; the integrands are analytic on the paths, and the analytic flag is checked).
Controls: f on [2,4] is 2e^g log(s-1)/s, and the [4,6] formula at s = 4 must equal it; F on [3,5] and the [5,7] formula must agree at s = 5;
an independent floating-point quadrature gave 1.04940e-4 and 4.99452e-6.
"""
from flint import arb, acb, ctx
import sys
ctx.prec = 256
eg = arb.const_euler().exp()
NAN = acb('nan')

def Phi(u, an=False):
    u = acb(u)
    z = 1 / u
    if an and not (z.real < 1):          # Li_2 has its cut on [1, inf); make sure the ball stays left of it
        return NAN
    lu = u.log(analytic=an)
    return lu * lu / 2 + z.polylog(2)

def G(x, an=False):                       # int_2^x log(u-1)/u du
    return Phi(x, an) - Phi(2)

def f_46(s):                               # 4 <= s <= 6
    s = arb(s)
    I = acb.integral(lambda t, an: G(t - 1, an) / t, 3, s - 1) if s > 4 else acb(0)
    return 2 * eg / s * (arb(s - 1).log() + I.real)

def F_57(s):                               # 5 <= s <= 7
    s = arb(s)
    def H(t, an):                          # int_{t+2}^{s-1} (1/u) log((u-1)/(t+1)) du
        return Phi(s - 1) - Phi(t + 2, an) - (t + 1).log(analytic=an) * ((s - 1) / (t + 2)).log(analytic=an)
    I = acb.integral(lambda t, an: (t - 1).log(analytic=an) / t * H(t, an), 2, s - 3) if s > 5 else acb(0)
    return 2 * eg / s * (1 + G(s - 1).real + I.real)

def F_35(s):                               # 3 <= s <= 5: F = (2e^g/s)(1 + int_2^{s-1} log(t-1)/t dt)
    s = arb(s)
    return 2 * eg / s * (1 + G(s - 1).real)

print(f"# python-flint {__import__('flint').__version__}, Arb at {ctx.prec} bits; python {sys.version.split()[0]}")
c1 = f_46(4) - 2 * eg * arb(3).log() / 4
c2 = F_57(5) - F_35(5)
print("control: f_[4,6](4) - 2e^g log 3/4 =", c1.str(10), "; F_[5,7](5) - F_[3,5](5) =", c2.str(10))
assert abs(c1) < arb('1e-60') and abs(c2) < arb('1e-60')
# control 2 (the DDE (sF)' = f(s-1) on [5,7] reduces, by Fubini, to int_2^{s-3} log(t-1)/t log((s-2)/(t+1)) dt = int_3^{s-2} G(t-1)/t dt;
# checked at s = 6.5, which exercises G and the log/Li_2 pieces of H independently)
sx = arb('6.5')
lhs = acb.integral(lambda t, an: (t - 1).log(analytic=an) / t * ((sx - 2) / (t + 1)).log(analytic=an), 2, sx - 3)
rhs = acb.integral(lambda t, an: G(t - 1, an) / t, 3, sx - 2)
print("control: DDE identity at s = 6.5, difference =", (lhs - rhs).real.str(10))
assert abs((lhs - rhs).real) < arb('1e-60')
one_minus_f6 = 1 - f_46(6)
F7_minus_1 = F_57(7) - 1
print("1 - f(6) =", one_minus_f6.str(20))
print("F(7) - 1 =", F7_minus_1.str(20))
ok1 = one_minus_f6 < arb('1.05e-4'); ok2 = F7_minus_1 < arb('5e-6')
print(f"1 - f(6) < 1.05e-4: {ok1};  F(7) - 1 < 5e-6: {ok2};  margins (lower ends): "
      f"{(arb('1.05e-4') - one_minus_f6).lower().str(6)}, {(arb('5e-6') - F7_minus_1).lower().str(6)}")
assert ok1 and ok2