"""Exact rational sensitivity of two-by-two Ax=b, infinity norm, no rounding test."""
from fractions import Fraction as F
import json
def inverse(A):
    a,b=A[0];c,d=A[1];det=a*d-b*c
    if not det:raise ValueError('singular')
    return [[d/det,-b/det],[-c/det,a/det]]
def mv(A,x):return [sum(a*b for a,b in zip(row,x)) for row in A]
def vn(x):return max(map(abs,x))
def mn(A):return max(sum(map(abs,row)) for row in A)
def case(e):
    A=[[F(1),F(1)],[F(1),1+e]];Ai=inverse(A);x=[F(1),F(1)];b=mv(A,x);db=[F(0),F(1,1000)]
    xp=mv(Ai,[v+d for v,d in zip(b,db)]);dx=[v-u for v,u in zip(xp,x)]
    k=mn(A)*mn(Ai);relb=vn(db)/vn(b);relx=vn(dx)/vn(x);bound=k*relb
    assert mv(A,xp)==[v+d for v,d in zip(b,db)] and relx<=bound
    return {'A':A,'inverse':Ai,'b':b,'delta_b':db,'x':x,'perturbed_x':xp,'condition_infinity':k,'relative_b_change_infinity':relb,'relative_x_change_infinity':relx,'condition_times_relative_b_bound':bound}
def demo():
    difficult=case(F(1,1000));better=case(F(1))
    assert difficult['condition_infinity']==F(4004001,1000) and difficult['perturbed_x']==[0,2]
    assert better['condition_infinity']==9 and better['relative_x_change_infinity']==F(1,1000)
    try:inverse([[F(1),F(1)],[F(1),F(1)]])
    except ValueError:rejected=True
    else:raise AssertionError('singular accepted')
    return {'near_singular_case':difficult,'comparison_case':better,'singular_rejected':rejected,'arithmetic':'exact Fraction; no floating point rounding measured','norm':'infinity, row-sum induced matrix norm','A_perturbation_tested':False}
if __name__=='__main__':print(json.dumps(demo(),default=str,indent=2))