"""Exact rational endpoints. No binary floating point enters a certificate."""
const Rat = Rational{BigInt}
rational(n::Integer, d::Integer=1) = BigInt(n) // BigInt(d)
rational(x::Rational) = rational(numerator(x), denominator(x))
const Q = BigInt(339949771344778)
const S = BigInt(1000000000000000)
const DOMAIN_SIZE = 32768
const q = Q // S
const LOWER = rational(4196,12343)
const UPPER = rational(7039,20706)
const PERIOD = div(S,2)
const MODULAR_INVERSE = BigInt(281419755786509)

struct Interval
    lo::Rat
    hi::Rat
    function Interval(lo::Rat, hi::Rat)
        lo <= hi || throw(ArgumentError("reversed interval"))
        new(lo,hi)
    end
end
Interval(lo::Union{Integer,Rational},hi::Union{Integer,Rational}) = Interval(rational(lo),rational(hi))
Interval(x::Union{Integer,Rational}) = Interval(x,x)
Base.:+(a::Interval,b::Interval) = Interval(a.lo+b.lo,a.hi+b.hi)
Base.:-(a::Interval) = Interval(-a.hi,-a.lo)
Base.:-(a::Interval,b::Interval) = a+(-b)
function Base.:*(a::Interval,b::Interval)
    v=(a.lo*b.lo,a.lo*b.hi,a.hi*b.lo,a.hi*b.hi)
    Interval(minimum(v),maximum(v))
end
function Base.inv(a::Interval)
    a.lo <= 0 <= a.hi && throw(DomainError(a,"interval contains zero"))
    Interval(inv(a.hi),inv(a.lo))
end
Base.:/(a::Interval,b::Interval) = a*inv(b)
for op in (:+,:-,:*,:/)
    @eval Base.$op(a::Interval,b::Union{Integer,Rational}) = Base.$op(a,Interval(b))
    @eval Base.$op(a::Union{Integer,Rational},b::Interval) = Base.$op(Interval(a),b)
end
function Base.:^(a::Interval,n::Integer)
    n >= 0 || return inv(a)^(-n)
    n == 0 && return Interval(1)
    if iseven(n) && a.lo <= 0 <= a.hi
        return Interval(0,max(abs(a.lo),abs(a.hi))^n)
    end
    v=(a.lo^n,a.hi^n)
    Interval(minimum(v),maximum(v))
end
Base.in(x::Union{Integer,Rational},a::Interval) = a.lo <= x <= a.hi
floor_rat(x::Rational) = fld(numerator(x),denominator(x))
ceil_rat(x::Rational) = cld(numerator(x),denominator(x))

"""For x >= 0, enclose sqrt(x) using an integer square root on a decimal grid."""
function sqrt_enclosure(x::Rat; digits::Int=100)
    x >= 0 || throw(DomainError(x,"square root of negative rational"))
    digits >= 1 || throw(ArgumentError("digits must be positive"))
    scale=big(10)^digits
    n=isqrt(fld(numerator(x)*scale^2,denominator(x)))
    # n² <= x*scale² < (n+1)² follows directly from the floor and isqrt.
    Interval(n//scale,(n+1)//scale)
end

"""Strict arctan(1/b) bounds from an alternating series, b >= 2."""
function atan_inverse_enclosure(b::Integer; digits::Int=100)
    b >= 2 || throw(ArgumentError("b must be at least two"))
    digits >= 1 || throw(ArgumentError("digits must be positive"))
    x=rational(1,b); x2=x*x; power=x; total=rational(0)
    target=rational(1,big(10)^(digits+3))
    j=0
    while true
        total += (iseven(j) ? 1 : -1)*power/(2*j+1)
        j += 1; power *= x2
        next=(iseven(j) ? 1 : -1)*power/(2*j+1)
        if abs(next) < target
            return (interval=Interval(min(total,total+next),max(total,total+next)),terms=j)
        end
    end
end

"""Machin's identity pi = 16 atan(1/5) - 4 atan(1/239); see proof P01."""
function pi_enclosure(; digits::Int=100)
    a=atan_inverse_enclosure(5;digits=digits+2)
    b=atan_inverse_enclosure(239;digits=digits+2)
    (interval=16*a.interval-4*b.interval,terms=(a.terms,b.terms))
end

function trig_point_enclosure(x::Rat, kind::Symbol; terms::Int=8)
    0 <= x <= 1 || throw(DomainError(x,"Taylor certificate requires 0 <= x <= 1"))
    terms >= 1 || throw(ArgumentError("positive term count required"))
    kind in (:sin,:cos) || throw(ArgumentError("expected sin or cos"))
    first=kind == :sin ? 1 : 0
    total=rational(0)
    for j in 0:terms-1
        m=first+2*j
        total += (iseven(j) ? 1 : -1)*x^m/factorial(big(m))
    end
    m=first+2*terms
    next=(iseven(terms) ? 1 : -1)*x^m/factorial(big(m))
    Interval(min(total,total+next),max(total,total+next))
end

function constants_enclosure(; digits::Int=100)
    root=sqrt_enclosure(rational(5);digits)
    phi=(1+root)/2
    pc=pi_enclosure(;digits); piI=pc.interval
    lam=Interval(phi.lo/(phi.lo+piI.hi),phi.hi/(phi.hi+piI.lo))
    eps=piI/648000; delta=eps/2
    # sin is increasing and cos decreasing on [0,1]. The proof is in P01.
    sinI=Interval(trig_point_enclosure(delta.lo,:sin).lo,
                  trig_point_enclosure(delta.hi,:sin).hi)
    cosI=Interval(trig_point_enclosure(delta.hi,:cos).lo,
                  trig_point_enclosure(delta.lo,:cos).hi)
    (root=root,phi=phi,pi=piI,lam=lam,mu=1-lam,epsilon=eps,delta=delta,
     s=sinI,c=cosI,atan_terms=pc.terms,digits=digits)
end

"""Certified decimal truncation. Refine rational bounds; never guess a floor.
The universal p >= 8 result is analytic (P08), not an infinite computation.
"""
function certified_truncation(p::Int; max_digits::Int=max(256,p+128))
    p >= 0 || throw(ArgumentError("precision must be nonnegative"))
    scale=big(10)^p; digits=max(32,p+8)
    while digits <= max_digits
        lam=constants_enclosure(;digits).lam
        lo=floor_rat(scale*lam.lo); hi=floor_rat(scale*lam.hi)
        lo == hi && return lo//scale
        digits *= 2
    end
    error("enclosures did not certify this decimal truncation; no value returned")
end