Listing: Listing

#!/usr/bin/env python3
"""Em03-H-v8 rational defect checker. New; unexecuted artifact."""
import argparse
from fractions import Fraction as F
import json
import math
from pathlib import Path
import shutil
from em03h_v8_tools import need, sha, atomic_json, execute

BITS = 48
SCALE = 1 << BITS
TERMS = 40

def down(x):
    x = F(x)
    return F((x.numerator*SCALE)//x.denominator,SCALE)

def up(x):
    return -down(-F(x))

class I:
    def __init__(self,lo,hi=None):
        lo,hi = F(lo),F(lo if hi is None else hi)
        need(lo <= hi,"Reversed interval")
        self.lo,self.hi = down(lo),up(hi)

    def __add__(self,y):
        if isinstance(y,C):
            return NotImplemented
        y = ri(y)
        return I(self.lo+y.lo,self.hi+y.hi)

    __radd__ = __add__

    def __neg__(self):
        return I(-self.hi,-self.lo)

    def __sub__(self,y):
        if isinstance(y,C):
            return NotImplemented
        return self+(-ri(y))

    def __rsub__(self,y):
        return ri(y)+(-self)

    def __mul__(self,y):
        # Mixed real/complex arithmetic is handled by C.__rmul__.
        if isinstance(y,C):
            return NotImplemented
        y = ri(y)
        vals = [a*b for a in (self.lo,self.hi) for b in (y.lo,y.hi)]
        return I(min(vals),max(vals))

    __rmul__ = __mul__

    def __truediv__(self,y):
        y = ri(y)
        need(not y.lo <= 0 <= y.hi,"Division through zero")
        return self*I(1/y.hi,1/y.lo)

    def bound(self):
        return max(abs(self.lo),abs(self.hi))

    def width(self):
        return self.hi-self.lo

def ri(x):
    return x if isinstance(x,I) else I(x)

class C:
    def __init__(self,re=0,im=0):
        self.re,self.im = ri(re),ri(im)

    def __add__(self,y):
        y = ci(y)
        return C(self.re+y.re,self.im+y.im)

    __radd__ = __add__

    def __neg__(self):
        return C(-self.re,-self.im)

    def __sub__(self,y):
        return self+(-ci(y))

    def __rsub__(self,y):
        return ci(y)+(-self)

    def __mul__(self,y):
        y = ci(y)
        return C(self.re*y.re-self.im*y.im,
                 self.re*y.im+self.im*y.re)

    __rmul__ = __mul__

    def bound(self):
        return self.re.bound()+self.im.bound()

    def width(self):
        return self.re.width()+self.im.width()

def ci(x):
    return x if isinstance(x,C) else C(x)

def exp_minus(x):
    x = F(x)
    need(x >= 0,"Negative exp-minus argument")
    if x >= BITS:
        return I(0,F(1,SCALE))
    y,s = x,0
    while y > 1:
        y /= 2
        s += 1
    total,term = F(1),F(1)
    for k in range(1,TERMS+1):
        term *= -y/k
        total += term
    error = 3*y**(TERMS+1)/math.factorial(TERMS+1)
    out = I(max(F(0),total-error),min(F(1),total+error))
    for _ in range(s):
        out = out*out
        out = I(max(F(0),out.lo),min(F(1),out.hi))
    return out

def arctan(x):
    x = F(x)
    total = sum((-1)**k*x**(2*k+1)/F(2*k+1) for k in range(TERMS))
    error = x**(2*TERMS+1)/F(2*TERMS+1)
    return I(total-error,total+error)

def pi_ball():
    return 16*arctan(F(1,5))-4*arctan(F(1,239))

def sine(x):
    x = F(x)
    total = sum((-1)**k*x**(2*k+1)/math.factorial(2*k+1)
                for k in range(TERMS))
    error = x**(2*TERMS+1)/math.factorial(2*TERMS+1)
    return I(total-error,total+error)

def cosine(x):
    x = F(x)
    total = sum((-1)**k*x**(2*k)/math.factorial(2*k)
                for k in range(TERMS))
    error = x**(2*TERMS)/math.factorial(2*TERMS)
    return I(total-error,total+error)

def phase(t):
    need(0 <= t.lo <= t.hi <= 1,"Phase domain")
    sa,sb = sine(t.lo),sine(t.hi)
    ca,cb = cosine(t.lo),cosine(t.hi)
    return C(I(cb.lo,ca.hi),I(-sb.hi,-sa.lo))

def kernel(t):
    need(t.lo >= 0,"Kernel domain")
    den = I(4)+t*t
    return C(I(4)/den,(-2*t)/den)

def bump_interval(a,b):
    a,b = max(F(a),F(1,4)),min(F(b),F(3,4))
    if a >= b:
        return I(0)
    def d(x):
        return (x-F(1,4))*(F(3,4)-x)
    dmin = min(d(a),d(b))
    dmax = F(1,16) if a <= F(1,2) <= b else max(d(a),d(b))
    lo = F(0) if dmin == 0 else exp_minus(1/dmin).lo
    hi = exp_minus(1/dmax).hi
    return I(lo,hi)

def rational(pair):
    need(len(pair) == 2 and int(pair[1]) > 0,"Rational representation")
    return F(int(pair[0]),int(pair[1]))

def nodes(values):
    result = []
    for value in values:
        need(len(value) == 2,"Complex node representation")
        result.append((rational(value[0]),rational(value[1])))
    return result

def interpolate(values,t):
    m = len(values)-1
    if t == 1:
        return values[-1]
    x = t*m
    j = x.numerator//x.denominator
    f = x-j
    return tuple((1-f)*values[j][k]+f*values[j+1][k] for k in (0,1))

def enclosure(values,a,b):
    va,vb = interpolate(values,a),interpolate(values,b)
    return C(I(min(va[0],vb[0]),max(va[0],vb[0])),
             I(min(va[1],vb[1]),max(va[1],vb[1])))

def text(x):
    x = F(x)
    return [str(x.numerator),str(x.denominator)]

def interval_text(x):
    return {"lo":text(x.lo),"hi":text(x.hi)}

def complex_text(x):
    return {"real":interval_text(x.re),"imag":interval_text(x.im)}

def sqrt_bounds(x):
    x = F(x)
    need(x >= 0,"Square-root domain")
    k = math.isqrt((x.numerator*SCALE*SCALE)//x.denominator)
    lo = F(k,SCALE)
    return lo,lo if lo*lo == x else F(k+1,SCALE)

def contains(interval,value):
    return interval.lo <= value <= interval.hi

def selftests(r):
    for name,value in (("I_times_C",I(2)*C(3,4)),
                       ("C_times_I",C(3,4)*I(2))):
        r.exact(name,contains(value.re,6) and contains(value.im,8))
    r.exact("mixed_addition",contains((I(2)+C(3,4)).re,5))
    r.exact("mixed_subtraction",contains((I(2)-C(3,4)).im,-4))
    r.exact("exp_zero",contains(exp_minus(0),1))
    r.exact("sin_zero",contains(sine(0),0))
    r.exact("cos_zero",contains(cosine(0),1))
    r.exact("kernel_zero_real",contains(kernel(I(0)).re,1))
    r.exact("kernel_zero_imag",contains(kernel(I(0)).im,0))
    lo,hi = sqrt_bounds(F(9,16))
    r.exact("rational_square_root",lo <= F(3,4) <= hi)

def calculate(args,root,r):
    global BITS,SCALE,TERMS
    need(args.bits >= 32 and args.refine >= 1 and args.pulse_cells >= 8,
         "Parameter ranges")
    need(args.refine & (args.refine-1) == 0,"Power-of-two refinement")
    BITS,SCALE = args.bits,1 << args.bits
    TERMS = max(32,(BITS+32)//2)
    selftests(r)
    path = Path(args.nodes)
    shutil.copyfile(path,root/"nodes-input.json")
    helper = Path(__file__).with_name("em03h_v8_tools.py")
    shutil.copyfile(helper,root/helper.name)
    r.meta["helper_sha256"] = sha(helper)
    r.meta["nodes_sha256"] = sha(path)
    data = json.loads((root/"nodes-input.json").read_text())
    contract = {"t0":1,"T":1,"Omega":1,"v":1,"w":0,"packet":"fixed_f"}
    r.exact("primary_contract",data["contract"] == contract)
    r.exact("piecewise_linear",
            data["interpolation"] == "continuous_piecewise_linear")
    a,b = nodes(data["A"]),nodes(data["B"])
    need(len(a) == len(b) and len(a) >= 9,"Node lengths")
    m = len(a)-1
    need(m & (m-1) == 0,"Original power-of-two grid")
    n = m*args.refine
    need(n <= SCALE,"Grid must be exactly representable at arithmetic bits")
    h = F(1,n)
    t = [F(j,n) for j in range(n+1)]
    pi = pi_ball()
    b0 = I(0)
    for j in range(args.pulse_cells):
        bump = bump_interval(F(j,args.pulse_cells),F(j+1,args.pulse_cells))
        b0 = b0+F(1,args.pulse_cells)*bump
    r.record("pulse_normalization",interval_text(b0))
    need(b0.lo > 0,"Pulse lower bound unresolved")
    kap = [pi*bump_interval(t[j],t[j+1])/(2*b0) for j in range(n)]
    av = [enclosure(a,t[j],t[j+1]) for j in range(n)]
    bv = [enclosure(b,t[j],t[j+1]) for j in range(n)]
    cv = [av[j]+bv[j]*kap[j] for j in range(n)]
    jv = [av[j]*kap[j] for j in range(n)]
    kg = {d:kernel(I(F(d-1,n),F(d+1,n))) for d in range(1,n)}
    kp = {d:phase(I(F(d-1,n),F(d+1,n))) for d in range(1,n)}
    partial_g,partial_p = kernel(I(0,h)),phase(I(0,h))
    epsilon = F(0)
    records = []
    for j in range(n):
        ia,ib = C(),C()
        for k in range(j):
            ia = ia+(kg[j-k]*cv[k])*h
            ib = ib+(kp[j-k]*jv[k])*h
        ia = ia+(partial_g*cv[j])*I(0,h)
        ib = ib+(partial_p*jv[j])*I(0,h)
        da = av[j]-kernel(I(t[j]+1,t[j+1]+1))+C(0,1)*ia
        db = bv[j]+C(0,1)*ib
        local = max(da.bound(),db.bound())
        epsilon = max(epsilon,local)
        record = {"cell":j,"dA":complex_text(da),"dB":complex_text(db),
                  "modulus_upper":text(local),
                  "rectangle_width":text(max(da.width(),db.width()))}
        records.append(record)
        r.record("cell_%d"%j,record)
    atomic_json(root/"cell_defects.json",records)
    denom = exp_minus(1+pi.hi/2).lo
    need(denom > 0,"Amplification enclosure unresolved")
    amplification = 1/denom
    radius = epsilon*amplification
    qr,qi = b[-1]
    modlo,modhi = sqrt_bounds(qr*qr+qi*qi)
    etalo = max(F(0),modlo-radius)**2
    etahi = min(F(1),(modhi+radius)**2)
    need(0 <= etalo <= etahi <= 1,"Inconsistent amplitude enclosure")
    verdict = ("meets_quarter_error" if etalo >= F(3,4) else
               "does_not_meet_quarter_error" if etahi < F(3,4) else
               "unresolved")
    sphase = min(F(1),radius/modlo) if modlo > radius else F(1)
    certificate = {
        "arithmetic":"exact_rational_with_outward_dyadic_rounding",
        "bits":BITS,"Taylor_terms":TERMS,
        "original_cells":m,"validation_cells":n,
        "pulse_cells":args.pulse_cells,
        "rounding_grid":text(F(1,SCALE)),
        "B0":interval_text(b0),
        "B0_relative_width_upper":text(b0.width()/b0.lo),
        "pi":interval_text(pi),
        "uniform_defect_upper":text(epsilon),
        "amplification_upper":text(amplification),
        "q_center":[text(qr),text(qi)],
        "center_modulus":{"lo":text(modlo),"hi":text(modhi)},
        "amplitude_radius":text(radius),
        "eta":{"lo":text(etalo),"hi":text(etahi)},
        "specified_error":{"lo":text(1-etahi),"hi":text(1-etalo)},
        "ideal_phase_verdict":verdict,
        "estimated_phase_error_upper":text(1-etalo*(1-sphase)),
        "JSON_NPZ_link":"requires separate successful inspection",
        "independently_inspected":False,
        "physical_implementation_certified":False}
    atomic_json(root/"certificate.json",certificate)
    return certificate

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--nodes",required=True)
    parser.add_argument("--out",required=True)
    parser.add_argument("--bits",type=int,default=48)
    parser.add_argument("--refine",type=int,default=1)
    parser.add_argument("--pulse-cells",type=int,default=512)
    args = parser.parse_args()
    execute(args,calculate,__file__,numerical=False,
            success="BOUND_COMPUTED_INSPECTION_REQUIRED")

if __name__ == "__main__":
    main()