numpy.log1p处理极小复数时的异常行为及修正方案

我在使用numpy.log1p函数计算极小复数的log(1+x)值时，得到了不符合预期的结果。理论上输出应与输入近似相等，但以下示例并非如此：

np.log1p(1e-14 * (1 + 1j))
Out[75]: (9.992007221626358e-15+9.9999999999999e-15j)
np.log1p(1e-15 * (1 + 1j))
Out[76]: (1.110223024625156e-15+9.999999999999989e-16j)
np.log1p(1e-16 * (1 + 1j))
Out[77]: 1e-16j

scipy.special的log1p函数工作正常，但我因numba需求必须使用numpy版本。当前环境：

• numpy 1.26.4
• Python 3.10.10

版本验证代码：

np.__version__
Out[78]: '1.26.4'

更新 #

正如相关回答指出，numpy的复数log1p实现较为粗糙。我找到了一种修正方法：

def log1p_corrected(x):
    """Precise calculation of log(1 + x) in the complex case"""
    re_x = np.real(x)
    im_x = np.imag(x)

    # |1 + x|**2 - 1
    norm2_min_1 = 2 * re_x + re_x**2 + im_x**2

    # real part of log(1+x)
    re_log1p = np.log1p(norm2_min_1) / 2

    # imaginary part of log(1+x)
    im_log1p = np.arctan2(im_x, 1 + re_x)

    return re_log1p + 1j * im_log1p

与numpy实现的主要区别在于对数实部（log|1+x|）的计算——我重写公式以使用对实数处理良好的np.log1p。

我注意到基于级数展开的替代实现对极小x表现优异且高效，但我希望避免因根据|x|值使用不同方法而引入微小不连续性。

内容的提问来源于stack exchange，提问作者Lorenzo Guerini