@triton.jit
def fused_recurrent_bwd_kernel(q, k, k_l2, dq, dk, dk_l2, dk_l2_partial_fwd,
    dk_l2_partial_bwd, dq_reflected, dk_reflected, T, D: tl.constexpr, BK:
    tl.constexpr):
    i_b, i_h = tl.program_id(0), tl.program_id(1)
    d_h = tl.zeros([BK, BK], dtype=tl.float32)
    offset = i_b * T * D + i_h * BK + tl.arange(0, BK) + (T - 1) * D
    p_q = q + offset
    p_k = k + offset
    p_k_l2 = k_l2 + offset
    p_dq_reflected = dq_reflected + offset
    p_dk_reflected = dk_reflected + offset
    p_dq = dq + offset
    p_dk = dk + offset
    p_dk_l2 = dk_l2 + offset
    p_dk_l2_partial_fwd = dk_l2_partial_fwd + offset
    p_dk_l2_partial_bwd = dk_l2_partial_bwd + offset
    for _ in range(T):
        b_q = tl.load(p_q).to(tl.float32)
        b_k = tl.load(p_k).to(tl.float32)
        b_k_l2 = tl.load(p_k_l2).to(tl.float32)
        b_dq_reflected = tl.load(p_dq_reflected).to(tl.float32)
        b_dk_reflected = tl.load(p_dk_reflected).to(tl.float32)
        d_h += b_q[None, :] * b_dq_reflected[:, None] + b_k[None, :
            ] * b_dk_reflected[:, None]
        b_dk_l2_partial_fwd = tl.load(p_dk_l2_partial_fwd).to(tl.float32)
        b_dk_l2 = -2 * tl.sum(b_dk_l2_partial_fwd[:, None] * d_h, axis=0)
        tl.store(p_dk_l2, b_dk_l2.to(p_dk_l2.dtype.element_ty))
        b_dk_l2_partial = tl.sum(d_h * b_k_l2[None, :], axis=1)
        d_h -= 2 * b_k_l2[None, :] * b_dk_l2_partial[:, None]
        tl.store(p_dk_l2_partial_bwd, b_dk_l2_partial.to(
            p_dk_l2_partial_bwd.dtype.element_ty))
        p_dq_reflected -= D
        p_dk_reflected -= D
        p_q -= D
        p_k -= D
        p_k_l2 -= D
        p_dk_l2_partial_fwd -= D
        p_dk_l2_partial_bwd -= D
        p_dk_l2 -= D
    tl.debug_barrier()
    offset = i_b * T * D + i_h * BK + tl.arange(0, BK)
    p_q = q + offset
    p_k = k + offset
    p_k_l2 = k_l2 + offset
    p_dq_reflected = dq_reflected + offset
    p_dk_reflected = dk_reflected + offset
    p_dq = dq + offset
    p_dk = dk + offset
    p_dk_l2 = dk_l2 + offset
    p_dk_l2_partial_bwd = dk_l2_partial_bwd + offset
    h = tl.zeros([BK, BK], dtype=tl.float32) + (tl.arange(0, BK)[:, None] ==
        tl.arange(0, BK)[None, :])
    for _ in range(T):
        b_k_l2 = tl.load(p_k_l2).to(tl.float32)
        b_dk_l2_partial = tl.load(p_dk_l2_partial_bwd).to(tl.float32)
        b_dk_l2 = -2 * tl.sum(h * b_dk_l2_partial[:, None], axis=0)
        b_dk_l2 += tl.load(p_dk_l2)
        tl.store(p_dk_l2, b_dk_l2.to(p_dk_l2.dtype.element_ty))
        tmp = tl.sum(h * b_k_l2[None, :], axis=1)
        h -= 2 * b_k_l2[None, :] * tmp[:, None]
        b_dq_reflected = tl.load(p_dq_reflected).to(tl.float32)
        b_dk_reflected = tl.load(p_dk_reflected).to(tl.float32)
        b_dq = tl.sum(b_dq_reflected[:, None] * h, axis=0)
        b_dk = tl.sum(b_dk_reflected[:, None] * h, axis=0)
        tl.store(p_dq, b_dq.to(p_dq.dtype.element_ty))
        tl.store(p_dk, b_dk.to(p_dk.dtype.element_ty))
        p_q += D
        p_k += D
        p_k_l2 += D
        p_dq_reflected += D
        p_dk_reflected += D
        p_dk_l2 += D
        p_dk_l2_partial_bwd += D
        p_dq += D
        p_dk += D
