import random, cmath
# Laurent polynomials in t as dict {exponent: complex}
def mul(a,b):
    c={}
    for i,u in a.items():
        for j,v in b.items(): c[i+j]=c.get(i+j,0)+u*v
    return c
def add(a,b,s=1):
    c=dict(a)
    for k,v in b.items(): c[k]=c.get(k,0)+s*v
    return c
def minval(a,tol=1e-9): return min(k for k,v in a.items() if abs(v)>tol)
def matmul(A,B): return [[add(mul(A[i][0],B[0][j]),mul(A[i][1],B[1][j])) for j in range(2)] for i in range(2)]
def det(A): return add(mul(A[0][0],A[1][1]),mul(A[0][1],A[1][0]),-1)
def isone(p,tol=1e-9): return all(abs(v)<tol for k,v in p.items() if k!=0) and abs(p.get(0,0)-1)<tol
I=1j
# --- slice at x_1 = diag(t,1/t): h in L^{<0}SL_2 (h(infty)=I), A = h x_1 symmetric with entries in t^{-2}O.
# claimed: h = [[1+c/t+c^2/t^2+c^3/t^3, e/t],[e/t^3, 1-c/t]], e^2 = -c^4, and A in O_2 for c != 0.
for trial in range(5):
    c=complex(random.uniform(-2,2),random.uniform(-2,2)); 
    for e in (I*c*c, -I*c*c):
        h=[[{0:1,-1:c,-2:c*c,-3:c**3},{-1:e}],[{-3:e},{0:1,-1:-c}]]
        x1=[[{1:1},{}],[{},{-1:1}]]
        A=matmul(h,x1)
        sym = all(abs(v)<1e-9 for v in add(A[0][1],A[1][0],-1).values())
        print('x1-slice: det h = 1:',isone(det(h)),' A symmetric:',sym,' min val of A entries:',min(minval(A[i][j]) for i in range(2) for j in range(2) if any(abs(v)>1e-9 for v in A[i][j].values())))
# --- slice at x_0: Y_2 = {h in L^{<0}SL_2 symmetric, entries in t^{-2}O}: three components claimed.
def h_of(a1,b1,a2,b2,d2): return [[{0:1,-1:a1,-2:a2},{-1:b1,-2:b2}],[{-1:b1,-2:b2},{0:1,-1:-a1,-2:d2}]]
for trial in range(4):
    p,r,mu=[complex(random.uniform(-2,2),random.uniform(-2,2)) for _ in range(3)]
    s=mu*mu*(p*p+r*r)
    h=h_of(2*mu*p*r, mu*(r*r-p*p), s*p*p, s*p*r, s*r*r)
    print('Y_main point: det=1:',isone(det(h)),' minval',minval(h[0][0]) if False else min(minval(h[i][j]) for i in range(2) for j in range(2)))
    a1,nu=[complex(random.uniform(-2,2),random.uniform(-2,2)) for _ in range(2)]
    for sgn in (1,-1):
        h=h_of(a1, sgn*I*a1, nu, sgn*I*nu, -nu)
        print('plane P_%+d point: det=1:'%sgn, isone(det(h)))
# --- stabilizer of x_m: g = [[a, b t^m],[-b t^-m, a]], a^2+b^2=1, b in t^m O: check g x_m g^T = x_m for m=1 numerically with b = t*beta
m=1; beta=complex(0.3,0.7)
b={1:beta}  # b = beta t
# a = sqrt(1-b^2) = 1 - b^2/2 - b^4/8 - ... truncated series
b2=mul(b,b); b4=mul(b2,b2)
a=add(add({0:1},b2,-0.5),b4,-0.125)
g=[[a,mul(b,{m:1})],[mul(b,{-m:-1}),a]]
xm=[[{m:1},{}],[{},{-m:1}]]
gT=[[g[0][0],g[1][0]],[g[0][1],g[1][1]]]
R=matmul(matmul(g,xm),gT)
print('g x_m g^T - x_m, entries with |coef|>1e-3 up to t^4 (truncation error expected at t^5+):')
for i in range(2):
    for j in range(2):
        d=add(R[i][j],xm[i][j],-1); print('  ',i,j,{k:round(abs(v),4) for k,v in d.items() if abs(v)>1e-3 and k<=4})