Python中高效且防溢出的双曲正弦计算：基于对数空间的优化策略

import numpy as np
import mpmath as mp

# 假设A1是一个包含大数值的NumPy数组
# 示例：
h = 100
b = 1000
N = 10
a = 1
y1 = np.linspace(0, h, 250, False)
y2 = np.linspace(h, b, 800)
y = np.concatenate((y1, y2))
chi_1n = np.arange(1, N + 1) * np.pi / a
y_zone1 = np.reshape(y[y <= h], (1, -1))
chi1_y1 = chi_1n.reshape(-1, 1) * y_zone1
A1 = chi1_y1 # 假设A1是经过计算得到的，可能包含大数值

# 使用NumPy直接计算（可能溢出）
# sinhZ1_numpy = np.sinh(A1)

# 使用mpmath逐元素计算（性能低下）
mp.dps = 16 # 设置精度
sinhZ1_mpmath = np.empty(A1.shape, dtype=object)
for i in range(A1.shape[0]):
    for j in range(A1.shape[1]):
        sinhZ1_mpmath[i, j] = mp.sinh(A1[i, j])

# 此时，sinhZ1_mpmath的计算速度远低于NumPy，即使NumPy可能溢出

import numpy as np
from scipy.special import logsumexp

def logsinh(x):
    """
    计算 log(sinh(x))，避免溢出并支持向量化。
    参数:
        x (np.ndarray): 输入数组。
    返回:
        np.ndarray: log(sinh(x)) 的值。
    """
    ones = np.ones_like(x)
    # logsumexp([x, -x], b=[ones, -ones], axis=0) 计算 log(1*e^x + (-1)*e^-x) = log(e^x - e^-x)
    return logsumexp([x, -x], b=[ones, -ones], axis=0) - np.log(2)

from mpmath import mp
import numpy as np
from scipy.special import logsumexp
import timeit

mp.dps = 16 # 设置mpmath精度

# 定义logsinh函数
def logsinh(x):
    ones = np.ones_like(x)
    return logsumexp([x, -x], b=[ones, -ones], axis=0) - np.log(2)

# 生成测试数据
rng = np.random.default_rng()
x_test = rng.random(size=(100000)) * 10000 # 包含大数值的测试数组

# 向量化mpmath.sinh用于对比
# 注意：mp.sinh本身不支持向量化，这里通过np.vectorize进行包装，
# 但其内部仍是逐元素调用mp.sinh，所以性能依然很低。
mpsinh_vec = np.vectorize(mp.sinh)
mplog_vec = np.vectorize(mp.log)

print("Timing logsinh(x):")
# 运行100次，取平均
time_logsinh = timeit.timeit(lambda: logsinh(x_test), number=100)
print(f"logsinh(x) 平均耗时: {time_logsinh/100:.6f} 秒")

print("\nTiming mpsinh(x) (vectorized wrapper):")
# 运行1次，因为mpmath通常非常慢
time_mpsinh = timeit.timeit(lambda: mpsinh_vec(x_test), number=1)
print(f"mpsinh_vec(x) 平均耗时: {time_mpsinh:.6f} 秒")

# 验证结果的数值准确性
# 将mpmath的结果取对数并转换为float64进行比较
np.testing.assert_allclose(logsinh(x_test), mplog_vec(mpsinh_vec(x_test)).astype(np.float64), rtol=5e-15)
print("\n数值准确性验证通过！")

# 验证负数参数的处理
x_negative = -x_test + 0j # 将x_test转换为复数，确保logsinh输出复数结果
# logsinh(-x) = log(sinh(-x)) = log(-sinh(x)) = log(sinh(x)) + i*pi
np.testing.assert_allclose(logsinh(x_negative), logsinh(x_test) + np.pi * 1j)
print("\n负数参数处理验证通过！")