#!/usr/bin/env python3
"""Finite exact checks accompanying the compact-existence pages.

Python 3, standard library only. Run with --help for parameters. Rational
intervals enclose Euler node values; no floating-point result is used as proof.
The infinite-dimensional existence arguments are explicitly not run here.
"""
import argparse
from fractions import Fraction as Q
from math import isqrt, factorial
import json

STATUSES = ('passed', 'failed', 'blocked', 'not_run')

def sqrt_interval(x, denominator=10**12):
    x = Q(x)
    if x < 0 or denominator <= 0:
        raise ValueError('sqrt requires nonnegative x and positive denominator')
    k = isqrt((x.numerator * denominator**2) // x.denominator)
    lower = Q(k, denominator)
    upper = lower if lower*lower == x else Q(k+1, denominator)
    if not lower*lower <= x <= upper*upper:
        raise ArithmeticError('square-root enclosure failed')
    return lower, upper

def exp_one_interval(n=12):
    # For the first omitted term 1/(n+1)!, successive ratios <= 1/(n+2).
    lower = sum((Q(1, factorial(k)) for k in range(n+1)), Q(0))
    tail = Q(1, factorial(n+1)) / (1-Q(1,n+2))
    return lower, lower + tail

def record(name, status, **data):
    if status not in STATUSES:
        raise ValueError('unknown status')
    return dict(name=name, status=status, **data)

def calculate(step=Q(1,1000), tolerance=Q(1,100), nodes=False, max_steps=100000):
    if step <= 0 or tolerance <= 0 or max_steps < 1:
        raise ValueError('step, residual tolerance and max steps must be positive')
    endpoint = Q(1,4)
    count = -(-endpoint // step)
    result = []
    if count > max_steps:
        result.append(record('euler_mesh', 'blocked', required_steps=count,
                             max_steps=max_steps, reason='No shortened trajectory substituted'))
    else:
        t, lo, hi, mesh = Q(0), Q(0), Q(0), []
        for j in range(count):
            length = min(step, endpoint-t)
            root_lo, _ = sqrt_interval(lo)
            _, root_hi = sqrt_interval(hi)
            lo, hi = lo+length*(t+root_lo), hi+length*(t+root_hi)
            t += length
            if not Q(0) <= lo <= hi <= 1:
                raise ArithmeticError('Euler enclosure left certified rectangle')
            mesh.append({'time':str(t),'length':str(length),
                         'value_lower':str(lo),'value_upper':str(hi)})
        data = {'steps':count, 'endpoint':str(t), 'last_length':mesh[-1]['length'],
                'last_value_interval':[str(lo),str(hi)],
                'intervals_checked':len(mesh),'arithmetic':'exact rational bounds'}
        if nodes:
            data['nodes'] = mesh
        result.append(record('euler_mesh', 'passed', **data))
    actual_h = min(step, endpoint)
    root_lo, root_hi = sqrt_interval(Q(5,4)*actual_h)
    residual_upper = (actual_h+root_hi)/4
    result.append(record('integral_residual_budget',
                         'passed' if residual_upper <= tolerance else 'failed',
                         maximum_step=str(actual_h),upper_bound=str(residual_upper),
                         requested_tolerance=str(tolerance),
                         meaning='failed means this sufficient bound does not certify the budget; not a lower bound on actual residual',
                         bound_target='integral equation residual, not distance to a chosen solution'))
    centers = [Q(0),Q(1),Q(2)]; y=Q(3,4); epsilon=Q(3,2)
    a = [max(Q(0),epsilon-abs(y-z)) for z in centers]
    denominator=sum(a); w=[x/denominator for x in a]
    point=sum((wi*z for wi,z in zip(w,centers)),Q(0))
    error=abs(point-y)
    good = denominator>0 and all(x>=0 for x in w) and sum(w)==1 and error<epsilon
    result.append(record('finite_net_weights','passed' if good else 'failed',
                         weights=list(map(str,w)),weighted_point=str(point),
                         error=str(error),epsilon=str(epsilon),actual_affine_dimension=1))
    e_lo,e_hi=exp_one_interval()
    selfmap_upper=1+(e_hi-1)/2
    derivative_upper=1+Q(3,4)*e_hi
    result.append(record('sine_integral_constants','passed' if selfmap_upper<2 else 'failed',
                         e_interval=[str(e_lo),str(e_hi)],
                         selfmap_upper=str(selfmap_upper),derivative_upper=str(derivative_upper)))
    z=Q(1,8); ratio_lower=16*(z-z**3/6)
    result.append(record('finite_noncontraction_witness', 'passed' if ratio_lower>1 else 'failed',
                         constant_input='1/64', comparison_input='0',evaluation_point='0',
                         output_to_input_ratio_lower=str(ratio_lower),
                         analytic_ingredient='sin(z) >= z-z^3/6 for 0<=z<=1, proved in companion page'))
    _,sqrt2_hi=sqrt_interval(Q(2))
    migrated_selfmap=Q(1,2)+Q(3,8)*sqrt2_hi
    result.append(record('nonsmooth_migration_constants','passed' if migrated_selfmap<2 else 'failed',
                         kernel_integral_max='3/2',selfmap_upper=str(migrated_selfmap),
                         common_lipschitz_upper=str(Q(1,2)+sqrt2_hi/4)))
    result.append(record('existence_uniqueness_and_iteration','not_run',
                         reason='Finite arithmetic does not prove compactness, extract an infinite subsequence, prove existence, or certify uniqueness/iteration convergence. Read the proofs.'))
    return {'status_vocabulary':list(STATUSES),'checks':result}

def parse_fraction(value):
    try:
        return Q(value)
    except (ValueError, ZeroDivisionError) as err:
        raise argparse.ArgumentTypeError("use a valid fraction with nonzero denominator") from err

def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--step',type=parse_fraction,default=Q(1,1000))
    parser.add_argument('--residual-tol',type=parse_fraction,default=Q(1,100))
    parser.add_argument('--max-steps',type=int,default=100000)
    parser.add_argument('--nodes',action='store_true')
    args=parser.parse_args()
    try:
        report=calculate(args.step,args.residual_tol,args.nodes,args.max_steps)
    except (ValueError, ArithmeticError) as err:
        report={'status_vocabulary':list(STATUSES),'checks':[record('input_or_arithmetic','failed',reason=str(err))]}
    print(json.dumps(report,ensure_ascii=False,indent=2))
    return 1 if any(x['status'] in ('failed','blocked') for x in report['checks']) else 0

if __name__=='__main__':
    raise SystemExit(main())