#!/usr/bin/env python3
"""Finite quantum linear-algebra certificates, Python standard library only."""
import cmath,itertools,json,math
from fractions import Fraction
CHECKS=[]
def check(name,ok,details=None):
 if not ok:raise AssertionError(name)
 CHECKS.append(dict(check=name,passed=True,details=details))
def eye(n):return [[complex(i==j)for j in range(n)]for i in range(n)]
def adj(a):return [[v.conjugate()for v in col]for col in zip(*a)]
def mul(a,b):return [[sum(x*y for x,y in zip(row,col))for col in zip(*b)]for row in a]
def mv(a,v):return [sum(x*y for x,y in zip(row,v))for row in a]
def scale(c,a):return [[c*x for x in row]for row in a]
def add(a,b):return [[x+y for x,y in zip(ar,br)]for ar,br in zip(a,b)]
def sub(a,b):return add(a,scale(-1,b))
def kron(a,b):return [[x*y for x in ar for y in br]for ar in a for br in b]
def normv(v):return math.sqrt(sum(abs(x)**2 for x in v))
def normalized(v):return [x/normv(v)for x in v]
def close(a,b,tol=1e-11):return all(abs(x-y)<tol for ar,br in zip(a,b)for x,y in zip(ar,br))
def closev(a,b,tol=1e-11):return all(abs(x-y)<tol for x,y in zip(a,b))
def unitary(a):return close(mul(adj(a),a),eye(len(a)))
def diag(xs):return [[complex(xs[i])if i==j else 0j for j in range(len(xs))]for i in range(len(xs))]
def blocks(a,b,c,d):return [ar+br for ar,br in zip(a,b)]+[cr+dr for cr,dr in zip(c,d)]
def top(a,n=2):return [row[:n]for row in a[:n]]
def power(a,n):
 out=eye(len(a))
 for _ in range(n):out=mul(out,a)
 return out
I=eye(2);X=[[0j,1+0j],[1+0j,0j]];Z=diag([1,-1]);H=scale(1/math.sqrt(2),[[1+0j,1+0j],[1+0j,-1+0j]]);S=diag([1,1j])
def exp_pauli(a,t):return add(scale(math.cos(t),I),scale(-1j*math.sin(t),a))
def opnorm2(a):
 b=mul(adj(a),a);u=b[0][0].real;v=b[1][1].real;w=b[0][1]
 return math.sqrt(max(0,(u+v+math.sqrt((u-v)**2+4*abs(w)**2))/2))
plus=[1/math.sqrt(2)]*2
# Hamiltonian simulation and the exact noncommuting two-term example.
K=add(X,Z);U=add(scale(math.cos(math.sqrt(2)),I),scale(-1j*math.sin(math.sqrt(2))/math.sqrt(2),K))
rows=[]
for r in [1,2,4,10,100]:
 V=power(mul(exp_pauli(Z,1/r),exp_pauli(X,1/r)),r);err=opnorm2(sub(V,U));assert err<=1/r+1e-12
 rows.append(dict(r=r,evolutions=2*r,error=err,bound=1/r))
check('X+Z exact exponential and Trotter resource table',unitary(U)and close(mul(K,K),scale(2,I)),rows)
# State preparation by conditional probability masses, with zero handling.
v=[1,1,2,0];mass=[abs(x)**2 for x in v];N=sum(mass);out=[]
for j in range(4):
 branch=j>>1;parent=sum(mass[branch*2:branch*2+2]);a=math.sqrt(parent/N);a*=math.sqrt(mass[j]/parent)if parent else 1
 out.append(a)
check('Prefix mass tree prepares (1,1,2,0)',closev(out,normalized(v)),{'mass_root':6,'children':[2,4],'amplitudes':out[:],'first_angle':2*math.acos(1/math.sqrt(3))})
out[2]*=1j
check('Separate diagonal phase prepares (1,1,2i,0)',closev(out,normalized([1,1,2j,0])))
# Hadamard-test branches, including complex expectation signs.
def inner(a,b):return sum(x.conjugate()*y for x,y in zip(a,b))
states=[[1+0j,0j],[0j,1+0j],list(map(complex,plus)),[1/math.sqrt(2),1j/math.sqrt(2)],[1/math.sqrt(3),1j*math.sqrt(2/3)]]
case_count=0
for psi in states:
 for u in [I,X,Z,S,H]:
  up=mv(u,psi);z=inner(psi,up)
  for phase,target in [(1,z.real),(-1j,z.imag)]:
   zero=[(a+phase*b)/2 for a,b in zip(psi,up)];one=[(a-phase*b)/2 for a,b in zip(psi,up)];p0=normv(zero)**2;p1=normv(one)**2
   assert abs(p0+p1-1)<1e-12 and abs(p0-p1-target)<1e-12
   case_count+=1
check('Hadamard Re/Im sign and probabilities',True,{'cases':case_count,'example_Re':.5,'example_Im':.5,'example_P0':.75})
# LCU exact big unitary and both outcomes.
B=[[1/math.sqrt(3),-math.sqrt(2/3)],[math.sqrt(2/3),1/math.sqrt(3)]]
SELECT=blocks(I,scale(0,I),scale(0,I),Z);W=mul(mul(kron(adj(B),I),SELECT),kron(B,I));A=add(I,scale(2,Z))
check('LCU PREP-SELECT-PREP† upper block',unitary(W)and close(top(W),scale(1/3,A)))
w=mv(W,plus+[0j,0j]);p=normv(w[:2])**2
check('LCU success 5/9 and failure 4/9',abs(p-5/9)<1e-12 and closev(normalized(w[:2]),normalized([3,-1]))and abs(normv(w[2:])**2-4/9)<1e-12,{'success':p,'failure':normv(w[2:])**2})
# Oblivious amplification, reflection sign R=I-2Pi.
R=diag([-1,-1,1,1]);zero=scale(0,I)
for u in [I,X,H,S]:
 w=blocks(scale(.5,u),scale(math.sqrt(3)/2,I),scale(math.sqrt(3)/2,I),scale(-.5,adj(u)))
 a=scale(-1,mul(mul(mul(mul(w,R),adj(w)),R),w))
 assert unitary(w)and close(top(a),u)and all(abs(a[i][j])<1e-12 for i in[2,3]for j in[0,1])
check('One-step oblivious amplification for four unknown-input unitaries',True,{'W_calls':2,'W_dagger_calls':1,'aux_reflections':2})
b=diag([.5,.25]);d=diag([math.sqrt(.75),math.sqrt(15/16)]);w=blocks(b,d,d,scale(-1,b));a=scale(-1,mul(mul(mul(mul(w,R),adj(w)),R),w));v=mv(a,plus+[0j,0j]);p=normv(v[:2])**2
check('Nonunitary OAA distorts different singular values',close(top(a),diag([1,11/16]))and abs(p-377/512)<1e-12 and closev(normalized(v[:2]),normalized([16,11])),{'success':p,'desired_second_coefficient':.5,'actual_second_coefficient':11/16})
# Block multiplication with independent ancillas.
a=diag([1,.5]);d=diag([0,math.sqrt(3)/2]);ua=blocks(a,d,d,scale(-1,a));b=mul(mul(H,a),H);db=mul(mul(H,d),H);ub=blocks(b,db,db,scale(-1,b))
check('Explicit four-dimensional block dilation',unitary(ua)and close(top(ua),a))
check('Reused-ancilla square gives I, not A²',close(top(mul(ua,ua)),I)and not close(I,mul(a,a)))
def embed_two(u,q0,q1,n=3):
 size=1<<n;out=[[0j]*size for _ in range(size)]
 for j in range(size):
  bits=[j>>(n-1-k)&1 for k in range(n)];jlocal=2*bits[q0]+bits[q1]
  for ilocal in range(4):
   new=bits[:];new[q0]=ilocal>>1;new[q1]=ilocal&1;i=sum(b<<(n-1-k)for k,b in enumerate(new));out[i][j]=u[ilocal][jlocal]
 return out
up=mul(embed_two(ua,0,2),embed_two(ub,1,2))
check('Independent ancillas give ordered AB block',unitary(up)and close(top(up),mul(a,b)),{'ancillas':2,'data_qubits':1})
# Qubitization sign conventions, endpoints and general-unitary extension.
a=diag([.5,-.5]);d=scale(math.sqrt(3)/2,I);s=blocks(a,d,d,scale(-1,a));walk=mul(diag([1,1,-1,-1]),s)
for j,lam in enumerate([.5,-.5]):
 subspace=[[walk[j][j],walk[j][j+2]],[walk[j+2][j],walk[j+2][j+2]]]
 assert close(subspace,[[lam,math.sqrt(1-lam*lam)],[-math.sqrt(1-lam*lam),lam]])
check('Qubitization two planes and arccos phases',True,{'phases_positive_lambda':math.pi/3,'phases_negative_lambda':2*math.pi/3})
for k in range(7):
 assert close(top(power(walk,k)),diag([math.cos(k*math.acos(.5)),math.cos(k*math.acos(-.5))]))
check('Qubitization powers give Chebyshev response',True,{'powers':list(range(7))})
send=blocks(diag([1,-1]),zero,zero,diag([-1,1]));wend=mul(diag([1,1,-1,-1]),send)
check('Qubitization endpoints are one-dimensional',close(top(wend),diag([1,-1])))
u=mul(diag([1,1,1j,1j]),s);sp=blocks(scale(0,eye(4)),u,adj(u),scale(0,eye(4)))
g=[[0j]*2 for _ in range(8)]
for j in range(2):g[j][j]=g[j+4][j]=1/math.sqrt(2)
check('Extra flag Hermitianizes a non-Hermitian signal oracle',not close(u,adj(u))and unitary(sp)and close(sp,adj(sp))and close(mul(mul(adj(g),sp),g),a))
# Complete QSP pairs and exact phase words.
def zphase(phi):return diag([cmath.exp(1j*phi),cmath.exp(-1j*phi)])
def signal(x):return [[complex(x),1j*math.sqrt(1-x*x)],[1j*math.sqrt(1-x*x),complex(x)]]
def reflect(x):return [[complex(x),complex(math.sqrt(1-x*x))],[complex(math.sqrt(1-x*x)),complex(-x)]]
def qsp(x,phases):
 m=zphase(phases[0])
 for phi in phases[1:]:m=mul(mul(m,signal(x)),zphase(phi))
 return m
qsp_cases=0
for x in [k/50 for k in range(-50,51)]:
 for phases,p,q in [([0,0],x,1),([0,0,0,0],4*x**3-3*x,4*x*x-1),([0,math.pi/3,-math.pi/3,0],x**3,x*x+cmath.exp(1j*math.pi/3))]:
  v=qsp(x,phases);expected=[[p,1j*q*math.sqrt(1-x*x)],[1j*complex(q).conjugate()*math.sqrt(1-x*x),complex(p).conjugate()]]
  assert close(v,expected)and abs(abs(p)**2+(1-x*x)*abs(q)**2-1)<1e-11
  qsp_cases+=1
check('QSP complete polynomial pairs and three explicit sequences',True,{'grid_sequence_pairs':qsp_cases})
for x in [-1,-.8,-.2,0,.3,.5,1]:
 v=qsp(x,[math.pi/3,0]);vn=qsp(x,[-math.pi/3,0]);assert abs((v[0][0]+vn[0][0])/2-x/2)<1e-12
check('QSP real part x/2 versus impossible full entry',abs(.5)**2!=1,{'full_entry_endpoint_failure':.5,'real_part_realizable':True})
# Non-Hermitian QSVT, explicit Julia dilation and phase conventions.
a=[[0j,.3+0j],[.8+0j,0j]];dl=diag([math.sqrt(.91),.6]);dr=diag([.6,math.sqrt(.91)]);u=blocks(a,dl,dr,scale(-1,adj(a)))
def projphase(phi):return diag([cmath.exp(1j*phi)]*2+[cmath.exp(-1j*phi)]*2)
v=mul(mul(mul(mul(mul(projphase(math.pi),u),projphase(-math.pi/6)),adj(u)),projphase(-5*math.pi/6)),u)
sv=mul(mul(a,adj(a)),a);a3=power(a,3)
check('QSVT odd x³ three-query sequence',unitary(u)and close(top(v),sv)and close(sv,[[0,.027],[.512,0]]),{'singular_value_cube':[[0,.027],[.512,0]],'matrix_cube':[[0,.072],[.192,0]]})
check('Singular-value cube differs from matrix cube',not close(sv,a3)and close(a3,[[0,.072],[.192,0]]))
vp=mul(mul(mul(projphase(math.pi/4),adj(u)),projphase(-math.pi/4)),u)
vn=mul(mul(mul(projphase(-math.pi/4),adj(u)),projphase(math.pi/4)),u)
realblock=scale(.5,add(top(vp),top(vn)))
check('QSVT even x² requires right space and real-part auxiliary',close(realblock,mul(adj(a),a))and close(realblock,diag([.64,.09])),{'query_count':2,'real_part_ancilla':1,'right_space_diagonal':[.64,.09]})
pert=scale(.99,a);delta=opnorm2(sub(a,pert));err=opnorm2(sub(mul(mul(a,adj(a)),a),mul(mul(pert,adj(pert)),pert)))
check('QSVT coefficient robustness bound for x³',err<=3*delta+1e-12,{'block_difference':delta,'transformed_difference':err,'bound':3*delta})
# HHL exact example, uncomputation and signed eigenvalues.
for vals in [[1,.5],[1,-.5]]:
 raw=[.5*plus[j]/vals[j]for j in range(2)];prob=normv(raw)**2;out=normalized(raw)
 assert abs(prob-5/8)<1e-12 and closev(out,normalized([1,2 if vals[1]>0 else -2]))
check('HHL signed reciprocal rotation and success probability',True,{'success':5/8,'Z_expectation':-3/5})
rho_target=[[.2,.4],[.4,.8]];rho_dephased=diag([.2,.8]);fidelity=sum(rho_target[i][j]*rho_dephased[j][i]for i in range(2)for j in range(2)).real
check('Omitting HHL uncomputation dephases the solution',abs(fidelity-.68)<1e-12,{'fidelity_with_target':fidelity,'trace_distance':.4})
# No-fast-forwarding family. Taylor action is only a finite numeric checker.
def path_h(bits):
 n=len(bits);size=2*(n+1);h=[[0j]*size for _ in range(size)]
 for j,x in enumerate(bits):
  w=math.sqrt((n-j)*(j+1))/n
  for b in range(2):h[2*(j+1)+(b^x)][2*j+b]=h[2*j+b][2*(j+1)+(b^x)]=w
 return h
def exp_action(h,t,v):
 term=list(map(complex,v));out=term[:]
 for k in range(1,140):
  term=[(-1j*t/k)*x for x in mv(h,term)];out=[a+b for a,b in zip(out,term)]
  if k>30 and normv(term)<1e-16:break
 return out
nf_cases=0;maxerr=0
for n in range(1,7):
 for bits in itertools.product([0,1],repeat=n):
  h=path_h(bits);size=len(h);assert close(h,adj(h));assert all(sum(abs(x)>0 for x in row)<=2 for row in h)
  initial=[1]+[0]*(size-1);out=exp_action(h,math.pi*n/2,initial);target=[0j]*size;target[2*n+(sum(bits)&1)]=(-1j)**n
  err=normv([a-b for a,b in zip(out,target)]);assert err<1e-9;maxerr=max(maxerr,err);nf_cases+=1
  for j in range(n+1):
   for b in range(2):
    for direction in [-1,1]:
     k=j+direction
     if 0<=k<=n:
      bit=bits[max(j,k)-1];dest=2*k+(b^bit);weight=math.sqrt((n-min(j,k))*(min(j,k)+1))/n
      assert abs(h[dest][2*j+b]-weight)<1e-12
check('Every hidden parity path for N1..6 reaches the correct endpoint',True,{'instances':nf_cases,'maximum_numeric_vector_error':maxerr,'bit_queries_per_sparse_oracle':2})
# Capstone off-diagonal matrix, LCU alpha1, inverse state and readout budget.
a=add(scale(.75,I),scale(.25,X));ainv=[[1.5+0j,-.5+0j],[-.5+0j,1.5+0j]];raw=mv(ainv,[1,0]);state=normalized(raw)
check('Capstone inverse and normalized state',close(mul(a,ainv),I)and closev(state,normalized([3,-1])),{'success_at_C_half':.25*normv(raw)**2,'Z_expectation':inner(state,mv(Z,state)).real,'X_expectation':inner(state,mv(X,state)).real})
prep=[[math.sqrt(.75),-.5],[.5,math.sqrt(.75)]];select=blocks(I,zero,zero,X);w=mul(mul(kron(adj(prep),I),select),kron(prep,I))
check('Capstone LCU encodes (3I+X)/4 at alpha1',unitary(w)and close(top(w),a))
check('Capstone multiplication, fresh-power target and OAA sign',close(mul(w,w),eye(4))and close(sub(scale(3,a),scale(4,power(a,3))),scale(-1,X))and close(power(a,3),add(scale(9/16,I),scale(7/16,X))),{'LCU_success':5/8,'inverse_success':5/8,'cubic_success':65/128})
m=math.ceil(2/.04**2*math.log(2/.01));eta=.001;pmin=(math.sqrt(5/8)-eta)**2;state_bound=2*eta/(math.sqrt(5/8)-eta)
check('Capstone state/readout and conditional success budgets',m==6623 and state_bound<.005 and 2*.005+.04<=.05+1e-12,{'independent_samples':m,'state_error_bound':state_bound,'ideal_expected_attempts':m/(5/8),'approx_success_lower_bound':pmin,'approx_expected_attempts_upper':m/pmin})
RESULT=dict(unit='e-quantum-linear-algebra',passed=True,check_count=len(CHECKS),checks=CHECKS)
if __name__=='__main__':print(json.dumps(RESULT,ensure_ascii=False,indent=2))