前段时间我写了一个代码，它试图在不使用库函数的情况下昨天，我正在审查旧代码，并试图使其尽可能快（并且正确）。到目前为止，这是我的尝试：
const double ee = exp(1);
double series_ln_taylor(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 )
n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(1 - x) = -x - x**2/2 - x**3/3... */
n = 1 - n;
now = term = n;
for ( i = 1 ; ; ){
lgVal -= now;
term *= n;
now = term / ++i;
if ( now < 1e-17 ) break;
}
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
在这里，我试图找到使得的对数值，它小于 1。此时，的泰勒展开式可以放心使用。
我最近对数值分析产生了兴趣，这就是为什么我不禁要问这个问题，这个代码段在实践中可以运行多快，同时又足够正确？我是否需要切换到其他方法，例如使用连分数，像这样？
C 标准库提供的函数几乎比这个实现快 5.1 倍。
更新 1 ：使用Wikipedia中提到的双曲 arctan 系列，计算似乎比 C 标准库 log 函数慢 2.2 倍。虽然，我没有广泛检查性能，而且对于更大的数字，我当前的实现似乎真的很慢。如果可以的话，我想检查我的两个实现的错误界限和平均时间范围广泛的数字。这是我的第二次努力。
double series_ln_arctanh(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
for ( i = 3 ; ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
任何建议或批评表示赞赏。
更新 2：根据下面的建议，我在这里添加了一些增量更改，这比标准库实现慢了大约 2.5 倍。但是，我这次仅对整数进行了测试，对于较大的数字，运行时间会增加。目前。我还不知道生成随机双数的技术，所以它还没有完全进行基准测试。为了使代码更健壮，我添加了对极端情况的更正。我所做的测试的平均误差约为。
double series_ln_better(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n == 0 ) return -1./0.; /* -inf */
if ( n < 0 ) return 0./0.; /* NaN*/
if ( n < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
/* the cutoff iteration is 650, as over e**650, term multiplication would
overflow. For larger numbers, the loop dominates the arctanh approximation
loop (with having 13-15 iterations on average for tested numbers so far */
for ( term = 1; term < n && lgVal < 650 ; term *= ee, lgVal++ );
if ( lgVal == 650 ){
n /= term;
for ( term = 1 ; term < n ; term *= ee, lgVal++ );
}
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
/* limiting the iteration for worst case scenario, maximum 24 iteration */
for ( i = 3 ; i < 50 ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}