复现Shapiro & Watson(1988)时Givens Rotation实现问题求助

我正在复现Shapiro & Watson（1988）的研究，需要通过Givens旋转施加长期约束。目标恒等式为：

C_inf@B=C_inf(P@Q_total)=(C_inf@P)@Q_total=A@Q_total

但我的代码运行后得到的C(1)@B与A@Q_total结果不一致，怀疑自己对Givens旋转的定义或应用存在错误，以下是相关代码：

import numpy as np

def givens_rotation(a, b):
    r = np.hypot(a, b)
    if r == 0:
        return 1.0, 0.0
    c = a / r
    s = b / r
    return c, s
   
def apply_givens(M, i, j, k):
    
    x = M[i, j]
    y = M[i, k]
    c, s = givens_rotation(x, y)

    n = M.shape[0]
    G = np.eye(n)
    G[j, j] = c
    G[j, k] = s
    G[k, j] = -s
    G[k, k] = c
       
    M[:] = M @ G
    return G
   
   
zero_constraints = [
    (0, 4),  # step 1: zero A[0,4]
    (0, 3),  # step 2: zero A[0,3]
    (0, 2),  # step 3: zero A[0,2]
    (0, 1),  # step 4: zero A[0,1]
    #
    (1, 4),  # step 5: zero A[1,4]
    (1, 3),  # step 6: zero A[1,3]
    (1, 2),  # step 7: zero A[1,2]
    #
    (2, 4),  # step 8: zero A[2,4]
    (2, 3)   # step 9: zero A[2,3]
]
   
   
for (i, j) in zero_constraints:
    k = n - 1  
    k = 4 if j < 4 else (3 if j < 3 else (2 if j < 2 else (1 if j < 1 else None)))
    
    G_jk = apply_givens(A, i, j, k)    

内容的提问来源于stack exchange，提问作者samuele ridolfi