0

首先，我使用的是自然对数的近似值。或查看此处(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，同时 e^1 = 2.71828182846。但不幸的是，自然对数的情况并非如此，我什至不知道自然对数的这个系列是如何派生的（我的意思是来自上面的链接）。而且我找不到关于这个系列的任何来源。

精度误差从何而来？

4

2 回答 2

3

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 求和。

于 2022-01-17T12:52:21.733 回答
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
于 2022-01-17T17:39:08.327 回答