首先，我使用的是自然對數的近似值。或查看此處(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。但不幸的是，自然對數的情況并非如此，我什至不知道自然對數的這個系列是如何派生的（我的意思是來自上面的鏈接）。而且我找不到關于這個系列的任何來源。
精度誤差從何而來？
uj5u.com熱心網友回復：
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 求和。
uj5u.com熱心網友回復：
這條線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
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/414607.html
標籤：