I had written a code a while ago which attempted to calculate $log(x)$ without using library functions. Yesterday, I was reviewing the old code, and I tried to make it as fast as possible, (and correct). Here's my attempt so far:
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;
}
Here I am trying to find $a$ so that $e^{a}$ is just over n, and then I add the logarithm value of $\frac{n}{e^{a}}$, which is less than 1. At this point, the Taylor expansion of $log(1\ - \ x)$ can be used without worrying.
I have recently grown an interest in numerical analysis, and that's why I can't help asking the question, how much faster this code segment can be run in practice, while being correct enough? Do I need to switch to some other methods , for example, using continued fraction, like this?
The $log(x)$ function provided with C standard library is almost 5.1 times faster than this implementation.
UPDATE 1: Using the hyperbolic arctan series mentioned in Wikipedia, the computation seems to be almost 2.2 times slower than C standard library log function. Though, I have not extensively checked the performance, and for larger numbers, my current implementation seems to be REALLY slow. I want to check both of my implementation for error bound and average time for a wide range of numbers if I can manage. Here is my second effort.
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;
}
Any suggestion or criticism is appreciated.
UPDATE 2: Based on suggestions made below, I have added some incremental changes here, which is around 2.5 times slower than standard library implementation. However, I have tested it only for integers $\leq1e^{8}$ this time, for larger numbers the runtime would increase. For now. I do not yet know of techniques to generate random double numbers $\leq1e^{308}$, so it's not yet fully benchmarked. To make the code more robust, I have added corrections for corner cases. The average error for the tests I made is around $4e^{-15}$.
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;
}