# 독립 oracle와 재현 수치 검토에 사용한 exact/rational·spectral·Gauss 검산 코드를 보존한다. 코드의 결과는 당시 정식화 검산이며 실제 solver의 전체 검증을 대신하지 않는다. 실제 실행·소스·한계는 [구현 검증](implementation-verification.md)과 [수치 검토](numerical-review.md)를 따른다. # 독립 numerical oracle scratch 이 파일은 numerical-review.md의 사전 계산을 재현하는 문서 부록이다. Solver 또는 production test가 아니며 외부 package·reference·출력 파일을 사용하지 않는다. Python 표준 라이브러리만 사용한다. Fraction은 oracle에서 입력 double의 exact 산술을 확인하는 용도이며 production arbitrary-precision dependency를 요구하지 않는다. T1 root에서 다음 명령은 각 Python block을 순서대로 메모리에서 실행한다. ```powershell @' from pathlib import Path import re p = Path('docs/shell-transient-heat-transfer/numerical-oracle-scratch.md') blocks = re.findall(r'```python\n(.*?)```', p.read_text(encoding='utf-8-sig'), re.S) assert len(blocks) == 5 scope = {} for block in blocks: exec(compile(block, str(p), 'exec'), scope) print('ALL SCRATCH BLOCKS COMPLETED') '@ | python - ``` 실행 exit는 0이어야 한다. Assertion 범위는 exact algebra와 독립 oracle 자체에 한정된다. 수렴/이산 오차의 출력 숫자 자체를 production tolerance로 사용하지 않고 본 보고서의 roundoff·conditioning·오차 전파 기준을 적용한다. NumPy 부재 probe exit1은 실행 실패 이력으로 보고서에 남겼으며 아래 코드에는 NumPy를 쓰지 않는다. ## A. Slab spectral·quadrature·공간/시간 oracle ```python import math,json def zeros(n,m=None): return [[0.0]*(n if m is None else m) for _ in range(n)] def tr(A):return list(map(list,zip(*A))) def mm(A,B):return [[sum(a*b for a,b in zip(row,col)) for col in zip(*B)]for row in A] def mv(A,x):return [sum(a*b for a,b in zip(row,x)) for row in A] def solve(A,b): a=[list(r)+[v]for r,v in zip(A,b)];n=len(b) for k in range(n): p=max(range(k,n),key=lambda i:abs(a[i][k]));a[k],a[p]=a[p],a[k] d=a[k][k] for j in range(k,n+1):a[k][j]/=d for i in range(n): if i!=k: d=a[i][k] for j in range(k,n+1):a[i][j]-=d*a[k][j] return [r[-1]for r in a] def chol(A): n=len(A);L=zeros(n) for i in range(n): for j in range(i+1): v=A[i][j]-sum(L[i][k]*L[j][k]for k in range(j)) L[i][j]=math.sqrt(v)if i==j else v/L[j][j] return L def eig(A): a=[r[:]for r in A];n=len(a);v=[[float(i==j)for j in range(n)]for i in range(n)] for z in range(100*n*n): p,q=max(((i,j)for i in range(n)for j in range(i+1,n)),key=lambda ij:abs(a[ij[0]][ij[1]])) if abs(a[p][q])<1e-15*sum(abs(a[i][i])for i in range(n)):break phi=.5*math.atan2(2*a[p][q],a[q][q]-a[p][p]);c=math.cos(phi);s=math.sin(phi) J=[[float(i==j)for j in range(n)]for i in range(n)] J[p][p]=J[q][q]=c;J[p][q]=s;J[q][p]=-s a=mm(mm(tr(J),a),J);v=mm(v,J) ix=sorted(range(n),key=lambda i:a[i][i]) return [a[i][i]for i in ix],[[v[j][i]for i in ix]for j in range(n)] def slabmat(ne): n=2*ne+1;C=zeros(n);K=zeros(n);l=1/ne for e in range(ne): for i in range(3): for j in range(3): C[2*e+i][2*e+j]+=l/30*[[4,2,-1],[2,16,2],[-1,2,4]][i][j] K[2*e+i][2*e+j]+=[[7,-8,1],[-8,16,-8],[1,-8,7]][i][j]/(3*l) return C,K C,K=slabmat(2);Cf=[r[1:]for r in C[1:]];Kf=[r[1:]for r in K[1:]] L=chol(Cf);Li=tr([solve(L,[float(i==j)for i in range(4)])for j in range(4)]) lam,V=eig(mm(mm(Li,Kf),tr(Li)));W=mm(tr(Li),V) yp=[-1-x for x in solve(Cf,[r[0]for r in C[1:]])] coef=mv(tr(W),mv(Cf,yp)) def spectral(t,n=None): amp=[a*(math.exp(-l*t)if n is None else math.exp(-n*math.log1p(l*t/n)))for a,l in zip(coef,lam)] T=[1]+[1+x for x in mv(W,amp)] if n is None:dT=[0]+mv(W,[-l*x for l,x in zip(lam,amp)]) else: oldamp=[x*(1+l*t/n)for l,x in zip(lam,amp)] old=[1]+[1+x for x in mv(W,oldamp)] dT=[(a-b)/(t/n)for a,b in zip(T,old)] E=sum(mv(C,T));R=mv(C,dT)[0]+mv(K,T)[0] q=[] for e,ss in [(0,-1),(0,0),(0,1),(1,-1),(1,0),(1,1)]: q.append(-4*sum(a*b for a,b in zip([ss-.5,-2*ss,ss+.5],T[2*e:2*e+3]))) return T,q,E,R def exact(t): M=80 T=[1-sum(4/((2*m+1)*math.pi)*math.sin((2*m+1)*math.pi*x/2)*math.exp(-((2*m+1)*math.pi/2)**2*t)for m in range(M))for x in [0,.25,.5,.75,1]] q=[sum(2*math.cos((2*m+1)*math.pi*x/2)*math.exp(-((2*m+1)*math.pi/2)**2*t)for m in range(M))for x in [0,.25,.5,.5,.75,1]] E=1-sum(8/((2*m+1)*math.pi)**2*math.exp(-((2*m+1)*math.pi/2)**2*t)for m in range(M)) return T,q,E,q[0] def errors(a,b):return [max(abs(x-y)for x,y in zip(a[0],b[0])),max(abs(x-y)for x,y in zip(a[1],b[1])),abs(a[2]-b[2]),abs(a[3]-b[3])] out={"lambda":lam,"projected_initial_free":[1+x for x in yp],"slab":[]} for t in [.001,.01,.05,.1,.5,1]: s=spectral(t);x=exact(t) out["slab"].append({"t":t,"penetration":math.sqrt(t),"spatial_error":errors(s,x),"BE4096_temporal":errors(spectral(t,4096),s),"BE8192_temporal":errors(spectral(t,8192),s),"spectral":s}) print(json.dumps(out,indent=2)) def gauss(n): x=[];w=[] for i in range(1,n+1): z=math.cos(math.pi*(i-.25)/(n+.5)) for _ in range(30): p0,p1=1.,z for j in range(2,n+1):p0,p1=p1,((2*j-1)*z*p1-(j-1)*p0)/j dp=n*(z*p1-p0)/(z*z-1);nz=z-p1/dp if abs(nz-z)<2e-16: z=nz;break z=nz x.append(z);w.append(2/((1-z*z)*dp*dp)) return list(zip(x,w)) def dot(x,y):return sum(a*b for a,b in zip(x,y)) def cross(a,b):return [a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]] def geometry(X,xi,eta): signs=[(-1,-1),(1,-1),(1,1),(-1,1)] N=[(1+s*xi)*(1+t*eta)/4 for s,t in signs] D=[[s*(1+t*eta)/4,t*(1+s*xi)/4]for s,t in signs] a=[sum(D[i][0]*X[i][j]for i in range(4))for j in range(3)] b=[sum(D[i][1]*X[i][j]for i in range(4))for j in range(3)] W=cross(a,b);J=math.sqrt(dot(W,W));norm=[z/J for z in W] aa,ab,bb=dot(a,a),dot(a,b),dot(b,b);det=J*J ra=[(bb*x-ab*y)/det for x,y in zip(a,b)];rb=[(aa*y-ab*x)/det for x,y in zip(a,b)] G=[[D[i][0]*ra[j]+D[i][1]*rb[j]for j in range(3)]for i in range(4)] return N,G,J,norm def surface(X,n): M=zeros(4);S=zeros(4) for xi,wx in gauss(n): for eta,wy in gauss(n): N,G,J,no=geometry(X,xi,eta) for i in range(4): for j in range(4): M[i][j]+=N[i]*N[j]*J*wx*wy;S[i][j]+=dot(G[i],G[j])*J*wx*wy return M,S Z0,Z1=slabmat(2) # h=1: Z0 physical mass already; Z1 physical normal stiffness already. def mats(X,n): M,S=surface(X,n);C=zeros(20);K=zeros(20) for i in range(4): for j in range(4): for a in range(5): for b in range(5): C[5*i+a][5*j+b]=M[i][j]*Z0[a][b] K[5*i+a][5*j+b]=S[i][j]*Z0[a][b]+M[i][j]*Z1[a][b] return C,K def fn(A):return math.sqrt(sum(x*x for r in A for x in r)) def rel(A,B):return fn([[a-b for a,b in zip(r,s)]for r,s in zip(A,B)])/fn(B) def response(X,C,K): A=[[c/.1+k for c,k in zip(cr,kr)]for cr,kr in zip(C,K)] T=solve(A,[float(i==0)for i in range(20)]) flux=[] for xi,eta in [(-1/math.sqrt(3),-1/math.sqrt(3)),(1/math.sqrt(3),-1/math.sqrt(3)),(1/math.sqrt(3),1/math.sqrt(3)),(-1/math.sqrt(3),1/math.sqrt(3))]: N,G,J,no=geometry(X,xi,eta) for off,s in [(0,-1),(0,0),(0,1),(2,-1),(2,0),(2,1)]: Ls=[s*(s-1)/2,1-s*s,s*(s+1)/2];Ds=[4*(s-.5),-8*s,4*(s+.5)] flux.append([-sum((G[i][j]*Ls[a]+no[j]*N[i]*Ds[a])*T[5*i+off+a]for i in range(4)for a in range(3))for j in range(3)]) return T,flux def rv(a,b):return math.sqrt(sum((x-y)**2 for x,y in zip(a,b)))/math.sqrt(sum(x*x for x in b)) cases={"affine":[[0,0,0],[1,0,0],[1,1,0],[0,1,0]],"distorted":[[0,0,0],[2,0,0],[1.5,1,0],[0,1,0]]} for w in [.1,.5,2]:cases["saddle"+str(w)]=[[x,y,w*x*y]for x,y in [(-1,-1),(1,-1),(1,1),(-1,1)]] ans={} for name,X in cases.items(): ds=[] for n in [3,5,7]: C,K=mats(X,n);T,q=response(X,C,K);ds.append((C,K,T,q)) rows=[] for a,b in [(ds[0],ds[1]),(ds[1],ds[2])]: rows.append([rel(a[0],b[0]),rel(a[1],b[1]),rv(a[2],b[2]),max(math.sqrt(dot([x-y for x,y in zip(q,r)],[x-y for x,y in zip(q,r)]))for q,r in zip(a[3],b[3]))/max(math.sqrt(dot(r,r))for r in b[3])]) ans[name]=rows print(json.dumps({"quadrature":ans},indent=2)) def direct54(X): C=zeros(20);K=zeros(20) for xi,wx in gauss(3): for eta,wy in gauss(3): N,G,J,no=geometry(X,xi,eta) for off in [0,2]: for s,ws in gauss(3): ls=[s*(s-1)/2,1-s*s,s*(s+1)/2] ds=[4*(s-.5),-8*s,4*(s+.5)] H=[0.]*20;B=zeros(20,3) for i in range(4): for a in range(3): p=5*i+off+a;H[p]=N[i]*ls[a] B[p]=[G[i][j]*ls[a]+no[j]*N[i]*ds[a]for j in range(3)] weight=wx*wy*ws*J/4 for i in range(20): for j in range(20): C[i][j]+=weight*H[i]*H[j];K[i][j]+=weight*dot(B[i],B[j]) return C,K for name,X in cases.items(): C,K=mats(X,3);Cd,Kd=direct54(X) print("54point",name,rel(C,Cd),rel(K,Kd)) assert rel(C,Cd)<1e-12 and rel(K,Kd)<1e-12 # Generalized eigen residual and C orthogonality, independently reassemble. C,K=slabmat(2);Cf=[r[1:]for r in C[1:]];Kf=[r[1:]for r in K[1:]] res=fn([[a-b*l for a,b,l in zip(kr,cr,lam)]for kr,cr in zip(mm(Kf,W),mm(Cf,W))])/fn(mm(Kf,W)) orth=fn([[v-float(i==j)for j,v in enumerate(row)]for i,row in enumerate(mm(mm(tr(W),Cf),W))]) print("slab eigen residual",res,"C orthogonality",orth) assert res<1e-12 and orth<1e-12 def spatial_bounded(n,steps=32,t=.001,ng=12): h=1/n;l=6*(1-math.cos(math.pi/n))/(h*h*(2+math.cos(math.pi/n))) af=math.exp(-steps*math.log1p(l*t/steps));ase=math.exp(-l*t);ae=math.exp(-math.pi**2*t) et=eg=etime=ese=0. for i in range(n): for s,w in gauss(ng): x=(i+(s+1)/2)*h;measure=w*h/2 v=((1-s)*math.cos(math.pi*i/n)+(1+s)*math.cos(math.pi*(i+1)/n))/2 g=(math.cos(math.pi*(i+1)/n)-math.cos(math.pi*i/n))/h et+=(af*v-ae*math.cos(math.pi*x))**2*measure eg+=(af*g+math.pi*ae*math.sin(math.pi*x))**2*measure etime+=((af-ase)*v)**2*measure ese+=(ase*v-ae*math.cos(math.pi*x))**2*measure return math.sqrt(et),math.sqrt(eg),math.sqrt(etime),math.sqrt(ese) rows=[] for n in [8,16,32,64]: a=spatial_bounded(n);b=spatial_bounded(n,ng=16) ratio=a[2]/a[3];quad=max(abs(x-y)/abs(x)for x,y in zip(a,b)) assert ratio<1/32 and quad<1e-8 rows.append(a) fluxbytes=n*n*24*3*32*8 tempbytes=(n+1)**2*5*33*8 print("bounded space",n,"T,grad,time,semi",a,"time/space",ratio,"quad",quad,"historybytes",fluxbytes+tempbytes) for a,b in zip(rows,rows[1:]): pt=math.log(a[0]/b[0],2);pg=math.log(a[1]/b[1],2) assert 1.95<=pt<=2.05 and .98<=pg<=1.03 print("bounded order",pt,pg) time_orders=[] for n in [16,32,64,128]: e1=abs(math.exp(-n*math.log1p(1/n))-math.exp(-1)) e2=abs(math.exp(-2*n*math.log1p(1/(2*n)))-math.exp(-1)) order=math.log(e1/e2,2) assert .97<=order<=1.02 time_orders.append(order) print("fixed spatial time orders",time_orders) ``` ## B. Exact Fraction integration·rank·jump ```python from fractions import Fraction as F def fmat(A,d):return [[F(v,d)for v in r]for r in A] def mul(a,b): c=[F(0)]*(len(a)+len(b)-1) for i,x in enumerate(a): for j,y in enumerate(b):c[i+j]+=x*y return c def der(a):return [i*a[i] for i in range(1,len(a))] def integ(a):return sum(F(x,i+1)for i,x in enumerate(a)) sh=[[1,-3,2],[0,4,-4],[0,-1,2]] assert [[integ(mul(a,b))for b in sh]for a in sh]==fmat([[4,2,-1],[2,16,2],[-1,2,4]],30) assert [[integ(mul(der(a),der(b)))for b in sh]for a in sh]==fmat([[7,-8,1],[-8,16,-8],[1,-8,7]],3) M=fmat([[4,2,1,2],[2,4,2,1],[1,2,4,2],[2,1,2,4]],36) S=fmat([[4,-1,-2,-1],[-1,4,-1,-2],[-2,-1,4,-1],[-1,-2,-1,4]],6) fac=[([1,-1],[1,-1]),([0,1],[1,-1]),([0,1],[0,1]),([1,-1],[0,1])] assert [[integ(mul(a,c))*integ(mul(b,d))for c,d in fac]for a,b in fac]==M assert [[integ(mul(der(a),der(c)))*integ(mul(b,d))+integ(mul(a,c))*integ(mul(der(b),der(d)))for c,d in fac]for a,b in fac]==S Z0=[[F(0)for _ in range(5)]for _ in range(5)] Z1=[[F(0)for _ in range(5)]for _ in range(5)] for off in [0,2]: for i in range(3): for j in range(3): Z0[i+off][j+off]+=F([[4,2,-1],[2,16,2],[-1,2,4]][i][j],30) Z1[i+off][j+off]+=F([[7,-8,1],[-8,16,-8],[1,-8,7]][i][j],3) C=[[M[i//5][j//5]*Z0[i%5][j%5]for j in range(20)]for i in range(20)] K=[[S[i//5][j//5]*Z0[i%5][j%5]+M[i//5][j//5]*Z1[i%5][j%5]for j in range(20)]for i in range(20)] def pivots(A): a=[r[:]for r in A];p=[] for k in range(len(a)): d=a[k][k];p.append(d) for i in range(k+1,len(a)): for j in range(k+1,len(a)):a[i][j]-=a[i][k]*a[k][j]/d return p assert all(x>0 for x in pivots(C)) assert all(x>0 for x in pivots([r[:19]for r in K[:19]])) assert all(sum(r)==0 for r in K) assert sum(map(sum,C))==2 cc=[[F(2),F(1)],[F(1),F(2)]];kk=[[F(1),F(-1)],[F(-1),F(1)]] old=[F(0),F(0)] for expected,react in [(F(-1,5),F(24,5)),(F(1,25),F(36,25))]: new=[(2*(2*old[0]+old[1])-1)/5,F(1)] r=[2*sum(cc[i][j]*(new[j]-old[j])for j in range(2))+sum(kk[i][j]*new[j]for j in range(2))for i in range(2)] assert new[0]==expected and r==[0,react] assert 3*sum(new[j]-old[j]for j in range(2))==react/2 old=new print("exact: C SPD 20/20; K rank19; K1=0; 1TC1=2; jump -1/5,1/25; Rc 24/5,36/25") ``` ## C. Time predicate와 경계 ```python import math,random,struct from fractions import Fraction as F def input_cell(x): v=F(x);prev=F(math.nextafter(x,-math.inf)) nxt=math.nextafter(x,math.inf) hi=(v+F(nxt))/2 if math.isfinite(nxt) else v+(v-prev)/2 even=not(struct.unpack('>Q',struct.pack('>d',x))[0]&1) return (prev+v)/2,hi,even def cell_overlap(t,d,m): a,b,ca=input_cell(t);c,e,cb=input_cell(d);c*=m;e*=m lo=max(a,c);hi=min(b,e) return lo0 and d>0 and type(inc)is int and 1<=inc<=2**32-1):return "input-error" if d>=t:return 1,False,t if F(t)/F(d)>inc+1:return "INC" _,e=math.frexp(t);q=math.ldexp(t,-e)/math.ldexp(d,-e);m=round(q) close=cell_overlap(t,d,m) assert close==dyadic_overlap(t,d,m) if close:n=m else: factor=F(2)**(-e+88);X=F(t)*factor;Y=F(d)*factor assert X.denominator==Y.denominator==1 X=int(X);Y=int(Y);f=math.floor(q) if f*Y>X:f-=1 elif (f+1)*Y<=X:f+=1 assert f*Y<=X<(f+1)*Y n=f if X==f*Y else f+1 if n>inc:return "INC" if not close:assert n==-(-F(t)//F(d)) delta=d if close else fm(-(n-1),d,t) assert 00):continue t=math.nextafter(t,math.inf if i%2 else -math.inf) if not math.isfinite(t):continue q=t/d if not math.isfinite(q)or q>2**32:continue mm=round(q) assert cell_overlap(t,d,mm)==dyadic_overlap(t,d,mm) count+=1 print("input-cell/fixed-dyadic independent comparisons",count) assert grid_cell(5.351455173733616,.17262758624947144)[0]==31 assert grid_cell(math.nextafter(.9,math.inf),.3)[0]==4 assert grid_cell(20*math.ulp(0.),6*math.ulp(0.))[0]==3 assert grid_cell(17*math.ulp(0.),5*math.ulp(0.))[0]==4 print("INPUT-CELL ASSERTIONS PASSED") for t,d in [(1.,.1),(math.ldexp(.9,-1000),math.ldexp(.3,-1000)),(math.ldexp(.9,1000),math.ldexp(.3,1000)),(1e-323,5e-324)]: print("input-cell scale",repr(t),repr(d),grid_cell(t,d)) assert grid_cell(10.,3.,3)=="INC" and grid_cell(10.,3.,4)[0]==4 assert grid_cell(float(2**32),1.,2**32-1)=="INC" for bad in [0,-1,1.5,2**32]:assert grid_cell(1.,.1,bad)=="input-error" for bad in [0.,-1.,math.inf,math.nan]:assert grid_cell(bad,.1)=="input-error" print("CURRENT TIME GRID ASSERTIONS PASSED") ``` ## D. 입력 rounding-cell의 독립 반증 검색 고정 seed의 검색은 폐기된 arithmetic-only predicate와 입력 rounding-cell의 차이를 검출한다. 이 출력은 predicate 실행 성공과 승인 시간 의미의 정합을 구별하기 위한 반증 evidence다. 상세 판정은 numerical-review.md §8을 따른다. ```python import math,random from fractions import Fraction as F def rad(x):return max(x-math.nextafter(x,-math.inf),math.nextafter(x,math.inf)-x)/2 def pred(t,d,m): q=t/d;p=m*d B=math.nextafter(d*rad(q)+rad(p),math.inf) return abs(float(F(t)-m*F(d)))<=B,B rng=random.Random(41);found=[] for z in range(20000): d=rng.uniform(.125,.25);m=rng.randint(2,100);t=m*d for dr in [-1,1]: tt=math.nextafter(t,math.inf if dr>0 else -math.inf) loT=(F(tt)+F(math.nextafter(tt,-math.inf)))/2 hiT=(F(tt)+F(math.nextafter(tt,math.inf)))/2 loD=(F(d)+F(math.nextafter(d,-math.inf)))/2 hiD=(F(d)+F(math.nextafter(d,math.inf)))/2 overlap=max(loT,m*loD)<=min(hiT,m*hiD) ok,B=pred(tt,d,m) if overlap and not ok: found.append((tt,d,m,float(abs(F(tt)-m*F(d))),B));break if found:break print("input-cell overlap rejected by retired arithmetic predicate:",found) t=math.nextafter(.9,math.inf);d=.3;m=3 print("specific .9 upper:",float(F(t)-m*F(d)),rad(t)+m*rad(d),"input-cell overlap",abs(F(t)-m*F(d))<=F(rad(t))+m*F(rad(d))) ``` ## E. Library ULP 가정 없는 Fourier interval oracle 고정 dyadic 320-bit outward enclosure, Machin π, Taylor sin/cos/exp와 analytic remainder를 사용한다. Python Fraction은 검증 generator 전용이다. 출력 interval의 outward binary64 hex endpoints를 C++ fixture에 고정할 수 있다. FESA solver와 C++ std::sin/std::exp의 ULP 보장을 가정하지 않는다. ```python from fractions import Fraction as F import math BITS=320;SCALE=1<0 or b[1]<0 return mul(a,(F(1)/b[1],F(1)/b[0])) def scale(a,c):return mul(a,interval(c)) def atan_bounds(q): a=sum((F(-1)**j)*q**(2*j+1)/(2*j+1) for j in range(64)) tail=q**129/129 return a,a+tail a,b=atan_bounds(F(1,5)),atan_bounds(F(1,239)) PI=interval(16*a[0]-4*b[1],16*a[1]-4*b[0]) assert PI[1]-PI[0]=0 if x[0]>=512:return interval(0,F(2)**-512) assert x[1]<=512 y=scale(x,F(1,1024));term=interval(1);s=term for j in range(1,65): term=scale(mul(term,neg(y)),F(1,j));s=add(s,term) tail=F(1,2)**65/math.factorial(65) s=add(s,interval(-tail,tail)) for _ in range(10):s=mul(s,s) return s def sincos_pi8(k): k%=16;quad=(k+2)//4;rem=k-4*quad x=scale(PI,F(rem,8));x2=mul(x,x) st=x;ss=x;ct=interval(1);cs=ct for j in range(1,32): st=scale(mul(st,neg(x2)),F(1,2*j*(2*j+1)));ss=add(ss,st) ct=scale(mul(ct,neg(x2)),F(1,(2*j-1)*2*j));cs=add(cs,ct) ss=add(ss,interval(-F(1,math.factorial(65)),F(1,math.factorial(65)))) cs=add(cs,interval(-F(1,math.factorial(64)),F(1,math.factorial(64)))) q=quad%4 return [(ss,cs),(cs,neg(ss)),(neg(ss),neg(cs)),(neg(cs),ss)][q] def series_interval(t): Ts=[interval(1)for _ in range(5)];qs=[interval(0)for _ in range(5)];E=interval(1) for m in range(80): odd=2*m+1 x=scale(mul(PI,PI),F(t)*odd*odd/4);dec=exneg(x) for j in range(5): si,co=sincos_pi8(odd*j) Ts[j]=sub(Ts[j],scale(div(mul(si,dec),PI),F(4,odd))) qs[j]=add(qs[j],scale(mul(co,dec),2)) E=sub(E,scale(div(dec,mul(PI,PI)),F(8,odd*odd))) assert max(hi-lo for lo,hi in Ts+qs+[E])