"""Shared rigorous helpers.  Transcendental functions: Arb (python-flint) balls, converted to float
intervals with one extra ulp of outward rounding.  Bulk arithmetic: numpy float64 with +,-,* only, each
result widened outward by one ulp (IEEE-754 round-to-nearest => the exact result lies in the widened interval)."""
import os, sys, math
import numpy as np
if os.environ.get('CERT_MOCK')=='1':
    print('*** WARNING: running with MOCK flint (float, NOT rigorous) - for debugging only ***',flush=True)
    sys.path.insert(0,os.path.dirname(os.path.abspath(__file__)))
    from mockflint import arb, acb, ctx
else:
    from flint import arb, acb, ctx
INF=np.inf
def dn(x): return np.nextafter(x,-INF)
def up(x): return np.nextafter(x,INF)
def lo(a): return float(dn(float(a.lower())))
def hi(a): return float(up(float(a.upper())))
U=2.0**-53
def sum_up(x):
    """rigorous upper bound for the exact sum of a float64 array x with x>=0 (any summation order)"""
    x=np.asarray(x,dtype=np.float64); n=x.size
    if n==0: return 0.0
    s=float(np.sum(x)); g=n*U/(1-n*U)
    return float(up(up(s)*(1+2*g)+1e-300))
def A(s): return arb(s)            # exact decimal string -> ball
PI=None
def init(prec):
    global PI
    ctx.prec=prec; PI=arb.pi(); return PI
def table_b(N,t,start=1):
    """[lo,hi] float bounds for b_m = exp((t/4) log^2 m), m=start..N"""
    n=N-start+1; blo=np.empty(n); bhi=np.empty(n); t4=t/4
    for m in range(start,N+1):
        L=arb(m).log(); b=(t4*L*L).exp()
        blo[m-start]=lo(b); bhi[m-start]=hi(b)
    return blo,bhi
def geom_blocks(nmax,r=1+2.0**-13,exact_below=4096):
    """integer block starts 1=B0<B1<...; blocks [B_k, B_{k+1}); exact (width 1) below exact_below"""
    B=list(range(1,min(exact_below,nmax)+1)); v=float(B[-1])
    while B[-1]<nmax:
        v=max(v*r,B[-1]+1); nb=int(math.ceil(v))
        if nb>nmax: break
        if nb>B[-1]: B.append(nb)
    B.append(nmax+1)
    return np.array(B,dtype=np.int64)