为什么我的自然对数函数如此不精确？

Question

首先，我使用这个 approximation of a natural log. Or look here (4.1.27) 来更好地表示公式。

这是我的实现：

constexpr double eps = 1e-12;

constexpr double my_exp(const double& power)
{
    double numerator = 1;
    ull denominator = 1;
    size_t count = 1;
    double term = numerator / denominator;
    double sum = 0;
    while (count < 20)
    {
        sum += term;
        numerator *= power;
        #ifdef _DEBUG
            if (denominator > std::numeric_limits<ull>::max() / count)
                throw std::overflow_error("Denominator has overflown at count " + std::to_string(count));
        #endif // _DEBUG
        denominator *= count++;
        term = numerator / denominator;
    }
    return sum;
}

constexpr double E = my_exp(1);

constexpr double my_log(const double& num)
{
    if (num < 1)
        return my_log(num * E) - 1;
    else if (num > E)
        return my_log(num / E) + 1;
    else
    {
        double s = 0;
        size_t tmp_odd = 1;
        double tmp = (num - 1) / (num + 1);
        double mul = tmp * tmp;
        while (tmp >= eps)
        {
            s += tmp;
            tmp_odd += 2;
            tmp *= mul / tmp_odd;
        }
        return 2 * s;
    }
}

你大概明白我为什么要实现这些功能了。基本上，我想实现一个 pow 函数。但是我的方法仍然给出了非常不精确的答案，例如 my_log(10) = 2.30256，但是根据 google (ln 10 ~ 2.30259).

my_exp() 非常精确，因为它的泰勒展开高度收敛。 my_exp(1) = 2.718281828459，同时根据 google，e^1 = 2.71828182846。但不幸的是，自然对数的情况不同，我什至不知道自然对数的这个系列是如何导出的（我的意思是从上面的链接）。而且我找不到关于这个系列的任何资源。

哪里来的精度误差？

Answer 1

if (num < 1) return my_log(num * E) - 1; 乘法不精确。乘以 2 更准确。当然，my_log(num) = my_log(2*num) - ln(2) 所以你需要更改 1.0 常量。

是的，现在您将在 -ln(2) 中出现舍入误差，而不是在 *E 中出现舍入误差。这通常没那么糟糕。

此外，您可以通过先检查 if (num<1/16) 然后使用 my_log(num) = my_log(16*num) - ln(16) 来保存重复的舍入错误。这只是一个舍入误差。

至于你的核心循环中的错误，我怀疑罪魁祸首是s += tmp;。这是一个重复的添加。您可以在那里使用 Kahan 求和。

Answer 2

行 tmp *= mul / tmp_odd; 表示每一项也除以所有先前项的分母，即 1, 1*3, 1*3*5, 1*3*5*7, ... 而不是公式所述的 1, 3, 5, 7, ...。

因此分子和分母应该独立计算：

double sum = 0;
double value = (num - 1) / (num + 1);
double mul = value * value;
size_t denom = 1;
double power = value;
double term = value;
while (term > eps)
{
    sum += term;
    power *= mul;
    denom += 2;
    term = power / denom;
}
return 2 * sum;

...

// Output for num = 1.5, eps = 1e-12
My func:   0.405465108108004513
Cmath log: 0.405465108108164385
             ------------

好多了！

将 epsilon 减少到 1e-18，我们达到了朴素求和的精度限制：

// Output for num = 1.5, eps = 1e-18
My func:   0.40546510810816444
Cmath log: 0.405465108108164385
             ---------------

Kahan-Neumaier 救援：

double sum = 0;
double error = 0;
double value = (num - 1) / (num + 1);
double mul = value * value;
size_t denom = 1;
double power = value;
double term = value;
while (term > eps)
{
    double temp = sum + term;
    if (abs(sum) >= abs(term))
        error += (sum - temp) + term;
    else
        error += (term - temp) + sum;
    sum = temp;
    power *= mul;
    denom += 2;
    term = power / denom;
}
return 2 * (sum + error);

...

// Output for num = 1.5, eps = 1e-18
My func:   0.405465108108164385
Cmath log: 0.405465108108164385

为什么我的自然对数函数如此不精确？

Why my natural log function is so imprecise?

c++

math

exponent

taylor-series

natural-logarithm