#抛物线法解方程
def MullerSolve(fx,x0,x1,x2,err1,err2,N):
    count = 0
    f0 = fx(x0)
    f1 = fx(x1)
    f2 = fx(x2)
    q = (x2 - x1) / (x1 - x0)
    p = 0
    a = 0
    b = 0
    c = 0
    while True:
        p = (x2 - x0) / (x1 - x0)

        a = q**2 * f0 - q*p*f1 + q*f2
        b = q**2 *f0 - p**2 *f1 + (p + q)*f2
        c = p*f2
        h1 = 0
        
        if  b.real < 0:
            h1 = -2 * c / (b - (b**2 - 4*a*c)**0.5)
        else:
            h1 = -2 * c / (b + (b**2 - 4*a*c)**0.5)
        x3 = x2 + h1 * (x2 - x1)
        f3 = fx(x3)
        k = err1 + 1
        if abs(f3) < 1:
            k = abs(x3 - x2)
        else:
            k = abs(x3 - x2) / abs(f3)

        if abs(f3) < err2 or k < err1:
            return x3, 1
        count += 1
        if count > N:
            return None, 0
        x0 = x1
        x1 = x2
        x2 = x3
        f0 = f1
        f1 = f2
        f2 = f3
        q = h1
#计算矩阵的行列式
def Det(A):
    if len(A) == 2:
        return A[0][0] * A[1][1] - A[0][1] * A[1][0]
    det = 0
    for c in range(len(A)):
        sub_matrix = [row[:c] + row[c+1:] for row in A[1:]]
        det += ((-1) ** c) * A[0][c] * Det(sub_matrix)
    return det

if __name__ == "__main__":
    A =[
        [1,0,1],
        [2,2,1],
        [-1,0,0]
    ]
    
    lam = []

    count = 0
    k = -100
    fx = lambda x: Det([[A[i][j] - x * (1 if i == j else 0) for j in range(len(A))] for i in range(len(A))])
    k = 10
    while len(lam) < len(A):
        for i in range(-k,k):
            re,sta = MullerSolve(fx, i, i + 1, i + 2, 1e-10, 1e-10, 100)
            if sta == 1:
                a = round(re.real, 9)
                b = round(re.imag, 9)
                re_t = complex(a, b)
                if re_t not in lam:
                    if re_t.imag != 0:
                        lam.append(re_t)
                        lam.append(re_t.conjugate())
                    else:
                        lam.append(a)
                if len(lam) == len(A):
                    break
        k *= 10
    p = abs(lam[0])
    for i in range(len(lam)):
        print(f"λ{i+1} = {lam[i]}")
        if abs(lam[i]) > p:
            p = abs(lam[i])
    
    print(f"谱半径 = {p:.3f}")