#列 行
def getIndexFromDownMatrix(col, row):
    if row > col:
        row, col = col, row
    return (1+col)*col//2+row


# 平方根法求解
def SqrtSolve(A,b):
    n = len(b)
    L = [0]*(n*(n+1)//2)
    for j in range(n):
        for i in range(j,n):
            if i == j:
                L[getIndexFromDownMatrix(i,j)] = A[getIndexFromDownMatrix(i,j)]
                for k in range(j):
                    L[getIndexFromDownMatrix(i,j)] -= L[getIndexFromDownMatrix(j,k)]**2
                L[getIndexFromDownMatrix(i,j)] = L[getIndexFromDownMatrix(i,j)]**0.5
            else:
                L[getIndexFromDownMatrix(i,j)] = A[getIndexFromDownMatrix(i,j)]
                for k in range(j):
                    L[getIndexFromDownMatrix(i,j)] -= L[getIndexFromDownMatrix(i,k)]*L[getIndexFromDownMatrix(j,k)]
                L[getIndexFromDownMatrix(i,j)] /= L[getIndexFromDownMatrix(j,j)]

    for i in range(n):
        for k in range(i):
            b[i] -= L[getIndexFromDownMatrix(i,k)]*b[k]
        b[i] /= L[getIndexFromDownMatrix(i,i)]

    for i in range(n-1,-1,-1):
        for k in range(i+1,n):
            b[i] -= L[getIndexFromDownMatrix(k,i)]*b[k]
        b[i] /= L[getIndexFromDownMatrix(i,i)]
    return b


if __name__ == "__main__":
    # 储存下三角矩阵 a11, a21, a22, a31, a32, a33 ...
    A = [4,2,2,-2,-3,14]


    b = [10,5,4]

    print("x:")
    print(SqrtSolve(A,b))