Implementation of pow(3m) function
George V. Reilly
gvr at cs.brown.edu
Sat Aug 4 13:11:37 AEST 1990
In article <13474 at smoke.BRL.MIL> gwyn at smoke.BRL.MIL (Doug Gwyn) writes:
>In article <1990Aug2.120706.25713 at bnrgate.bnr.ca> mleech at bcarh342.bnr.ca (Marcus Leech) writes:
>>I'm wondering about the efficiency of "standard" implementations of this
>> (and other) math library functions for C. Does pow(3m) conventionally
>> use a logarithm table, or is it (ugh!) iterative, or some B.O.T?
>
>All common implementations of pow(x,y) basically perform exp(log(x)*y),
>but they add some twists designed to avoid unnecessary overflow, or in
>the case of 4.3BSD, to implement the IEEE FP funny numbers. Thus there
>is often a considerable amount of additional bookkeeping computation
>beyond the log() and exp() evaluations that one would expect.
There are many common special cases where calling pow is unwise:
when y = 0.5, use sqrt instead; when y is a small positive or negative
integer, just do the multiplications by hand; when y is of the form
1/k (i.e., you're trying to find the k-th root of x), use the
Newton-Raphson method instead.
The following routine can be used to calculate x^n, where n is a non-negative
integer. Instead of requiring the n-1 multiplications that the naive
algorithm would use, it requires lg n multiplications.
double
power(double x, unsigned n)
{
double t = 1.0;
while (n > 0)
if (n & 1) {
t *= x; n--;
} else {
x *= x; n >>= 1;
}
return t;
}
Caveat: I would expect any decent implementation of pow to do this, though
this power function should be slightly faster.
If n is known in advance, just code the multiplications by hand,
using the above algorithm; e.g., to calculate x^6:
temp = x*x; temp *= x; temp *= temp;
The Newton-Raphson method can be used to find the roots of an equation
f(x) = 0 in an iterative fashion. Given x_i, the approximation to
the root at the i-th iteration, we compute
x_{i+1} = x_i - f(x_i) / f'(x_i)
It's important to make a good guess for x_0, or this may fail to converge,
or take a large number of iterations to converge---see any good calculus
book or numerical analysis text for further details.
To calculate n^(1/k), the k-th root of n, we use f(x) = x^k - n, which
yields f'(x) = k * x ^ (k-1), so
x_{i+1} = x_i - f(x_i) / f'(x_i)
= x_i - (x_i ^ k - n) / (k * x_i ^ (k-1))
= ((k - 1)/k) * x_i + (n/k) * (1/(x_i ^ (k-1)))
For example, when k = 5, this becomes
x_{i+1} = 0.8 * x_i + n' /(x_i ^ 4)
where n' = 0.2 * n and x_i ^ 4 should be calculated as (x_i ^ 2) ^ 2.
A reasonably good value for x_0, the initial guess, is n/k, when |n| > 1.0,
and min(n*k, 1.0), otherwise.
Caveat: you should find that the sqrt library function is faster
than a general-purpose Newton-Raphson root-finding routine.
________________
George V. Reilly gvr at cs.brown.edu
uunet!brunix!gvr gvr at browncs.bitnet Box 1910, Brown U, Prov, RI 02912
More information about the Comp.lang.c
mailing list