优化NumPy中的动态衰减累加和：Numba、Cython与纯NumPy方案

import numpy as np

def f_python(x, d):
    """
    纯Python循环实现动态衰减累加和。
    """
    result = np.empty_like(x)
    result[0] = x[0]
    for i in range(1, x.shape[0]):
        result[i] = result[i-1] * d[i] + x[i]
    return result

import numba
import numpy as np

@numba.jit
def f_numba(x, d):
    """
    使用Numba JIT编译优化的动态衰减累加和。
    """
    result = np.empty_like(x)
    result[0] = x[0]
    for i in range(1, x.shape[0]):
        result[i] = result[i-1] * d[i] + x[i]
    return result

# 以下代码需在Jupyter Notebook或IPython环境中运行，或保存为.pyx文件编译
# %%cython
import numpy as np
cimport numpy as np

cpdef np.ndarray[np.float64_t, ndim=1] f_cython(np.ndarray[np.float64_t, ndim=1] x, np.ndarray[np.float64_t, ndim=1] d):
    """
    使用Cython预编译优化的动态衰减累加和。
    """
    cdef:
        int i = 0
        int N = x.shape[0]
        np.ndarray[np.float64_t, ndim=1] result = np.empty_like(x)
    result[0] = x[0]
    for i in range(1, N):
        result[i] = result[i-1] * d[i] + x[i]
    return result

import numpy as np

def f_numpy(x, d):
    """
    纯NumPy分解实现动态衰减累加和（可能存在数值不稳定性）。
    """
    # 确保d[0]不为0，或者根据实际业务逻辑处理
    # 这里为了简化，假设d[0] = 1，并从d[1:]开始累乘
    # 为了与原始循环行为保持一致，需要调整d的累积乘积
    # 一个更准确的累积乘积P应为 P_0=1, P_i = d_i * P_{i-1}
    # 或者 P_i = d_1 * d_2 * ... * d_i

    # 构造一个包含1的d_prime，使得cumprod从1开始
    d_prime = np.concatenate(([1.], d[1:])) 

    # 计算累积乘积 P_i = d_1 * ... * d_i
    # 这里的result实际上是累积乘积 P_i
    # 如果d[0]是有效衰减因子，则需要更复杂的处理
    # 假设d[0] = 1，使得P[0] = 1

    # 修正：为了匹配 c[i] = P[i] * sum(x[k]/P[k])
    # P[0] = 1
    # P[i] = d[1] * d[2] * ... * d[i] for i > 0
    # 这里的d数组是原始的d，d[0]可能不是1
    # 假设d[0]是有效衰减因子，那么P[0] = d[0]
    # P[i] = d[0] * d[1] * ... * d[i]

    # 实际上，如果按照 c[i] = c[i-1] * d[i] + x[i]
    # 那么 P[i] = d[i] * d[i-1] * ... * d[1]
    # 而 P[0] = 1

    # 更直接的分解方式是：
    # 设 p_i = d_i * d_{i-1} * ... * d_1
    # c_i = x_i + d_i x_{i-1} + d_i d_{i-1} x_{i-2} + ... + d_i d_{i-1} ... d_1 x_0
    # c_i = p_i * (x_i/p_i + x_{i-1}/p_{i-1} + ... + x_0/p_0)
    # 其中 p_0 = 1, p_i = d_i * p_{i-1}

    # 重新构建累积乘积 P
    P = np.cumprod(d)

    # 原始答案中的 f_numpy 实现
    # 假设 d[0] 应该为 1
    # 如果 d[0] 为 1，则 P[0] = 1, P[1] = d[1], P[2] = d[1]*d[2], ...
    # 那么 f_numpy 的实现是：
    # result = np.cumprod(d)
    # return result * np.cumsum(x / result)
    # 这假设了 d 数组的第一个元素用于累积乘积的起始，
    # 且 x[0] / P[0] + x[1] / P[1] + ...
    # 这种形式需要 d[0] != 0。

    # 鉴于原始问题中的 d[0] 可能不是1，且循环是 c[i] = c[i-1] * d[i] + x[i]
    # 这里的分解式应为：
    # 令 P_k = d_1 * d_2 * ... * d_k (P_0 = 1)
    # 那么 c_i = P_i * sum_{k=0 to i} (x_k / P_k)
    # 这需要一个辅助数组 P，其中 P[0]=1，P[k]=d[1]*...*d[k]

    # 考虑到原始答案中的 f_numpy 实现
    # result = np.cumprod(d)
    # return result * np.cumsum(x / result)
    # 这个实现是基于 P[k] = d[0] * d[1] * ... * d[k] 的
    # 当 d[0] 参与累积乘积时，这与原始循环 c[0] = x[0] 的语义可能不完全一致
    # 例如，如果 d[0]=0.5, d[1]=0.6, x[0]=10, x[1]=5
    # c[0] = 10
    # c[1] = c[0] * d[1] + x[1] = 10 * 0.6 + 5 = 11
    # f_numpy:
    # P = [0.5, 0.3]
    # x/P = [10/0.5, 5/0.3] = [20, 16.66]
    # cumsum(x/P) = [20, 36.66]
    # P * cumsum(x/P) = [0.5*20, 0.3*36.66] = [10, 11]
    # 这种情况下，结果是匹配的。
    # 因此，原始答案中的 f_numpy 实现是正确的，但它可能在数值上不稳定。

    result = np.cumprod(d)
    return result * np.cumsum(x / result)

import numpy as np

def f_numpy_stable(x, d):
    """
    数值稳定的纯NumPy实现动态衰减累加和（对数域计算）。
    假设 d 中的所有元素都大于0。
    """
    # 计算 log(P_i)
    p_log = np.cumsum(np.log(d))

    # 计算 log(x_k / P_k) = log(x_k) - log(P_k)
    term_log = np.log(x) - p_log

    # 计算 log(sum(exp(log(x_k) - log(P_k))))
    sum_exp_log = np.logaddexp.accumulate(term_log)

    # 最终结果 c_i = exp(log(P_i) + log(sum_exp_log))
    return np.exp(p_log + sum_exp_log)