"""Selected attributed manuscript enclosure routines; filesystem output removed."""
from fractions import Fraction as F
from math import factorial
import numpy as np
from source import exact_source
from validated import I,S,ceildiv

def scalar_floor():
    A,B,C=map(F,['.0008','.3228','.099']);dt=F(1,10);p=12
    def jets(z):
        cs=[z]
        for n in range(p):cs.append(((I(A) if n==0 else I())-I(B)*cs[n]-I(C)*sum((cs[j]*cs[n-j] for j in range(n+1)),I()))/(n+1))
        return cs
    cube=I(F(0),None);cube=I(0,ceildiv(4*S,5))
    cp=jets(cube)[p];M=max(abs(cp.l),abs(cp.u));rem=ceildiv(M,10**p)
    center=0;err=0
    for _ in range(200):
        cs=jets(I(center,center));val=cs[p-1]
        for k in range(p-2,-1,-1):val=val*I(dt)+cs[k]
        nxt=min(4*S//5,max(0,(val.l+val.u)//2))
        err+=max(abs(nxt-val.l),abs(nxt-val.u))+rem;center=nxt
    bounds=[F(center-err,S),F(center+err,S)]
    assert F('.0024725782')<bounds[0]<bounds[1]<F('.0024725784')
    # Independent floating closed-form check is diagnostic only.
    D=float(B*B+4*A*C)**.5;ap=2*float(A)/(float(B)+D);am=-(float(B)+D)/(2*float(C))
    diagnostic=ap*(1-np.exp(-20*D))/(1-ap/am*np.exp(-20*D))
    out=dict(bounds=list(map(str,bounds)),display=list(map(float,bounds)),closed_form_N=diagnostic,dt=str(dt),order=p,bits=110,local_remainder=str(F(rem,S)),error=str(F(err,S)))
    return out

def mean_matrix(e,c,eps):
    d,B,C=exact_source(e,c,eps);A=[row[:6] for row in B[:6]]
    for i in range(6):
        for j,k,p in C[i]:
            if j<6:A[i][j]+=p
            if k<6:A[i][k]+=p
    return A

def linear_flow(phases,v):
    zs=[I(x) for x in v];center=[(x.l+x.u)//2 for x in zs];err=max(max(c-x.l,x.u-c) for c,x in zip(center,zs))
    dt=F(1,20);p=16;growth=1/(1-F('.09')*dt)
    # All exact trajectories bounded by 9*max(v) <1000 over total20.
    assert max(v)*9<1000 and min(v)>=0
    remainder=F(1000)*3**p*dt**p/factorial(p);rem=ceildiv(remainder.numerator*S,remainder.denominator)
    for e,c,t in reversed(phases):
        A=mean_matrix(e,c,F('.1'));assert all(sum(abs(a) for a in row)<=3 and sum(row)<=F('.09') for row in A)
        mat=[[(j,I(a)) for j,a in enumerate(row) if a] for row in A]
        for _ in range(int(t/dt)):
            assert max(center)<=950*S
            cs=[[I(x,x) for x in center]]
            for k in range(p-1):cs.append([sum((a*cs[k][j] for j,a in row),I())/(k+1) for row in mat])
            val=cs[-1]
            for k in range(p-2,-1,-1):val=[a*I(dt)+b for a,b in zip(val,cs[k])]
            nxt=[min(1000*S,max(0,(a.l+a.u)//2)) for a in val]
            local=max(max(abs(c-a.l),abs(c-a.u)) for c,a in zip(nxt,val))+rem
            err=ceildiv(err*growth.numerator,growth.denominator)+local;center=nxt
    return [F(center[5]-err,S),F(center[5]+err,S)]