Path: utzoo!utgpu!news-server.csri.toronto.edu!mailrus!cs.utexas.edu!uunet!brunix!gvr
From: gvr@cs.brown.edu (George V. Reilly)
Newsgroups: comp.lang.c
Subject: Re: Implementation of pow(3m) function
Keywords: math-library exponentiation C
Message-ID: <46584@brunix.UUCP>
Date: 4 Aug 90 03:11:37 GMT
References: <1990Aug2.120706.25713@bnrgate.bnr.ca> <13474@smoke.BRL.MIL>
Sender: news@brunix.UUCP
Reply-To: gvr@cs.brown.edu (George V. Reilly)
Organization: Brown University Department of Computer Science
Lines: 70

In article <13474@smoke.BRL.MIL> gwyn@smoke.BRL.MIL (Doug Gwyn) writes:
>In article <1990Aug2.120706.25713@bnrgate.bnr.ca> mleech@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@cs.brown.edu
uunet!brunix!gvr   gvr@browncs.bitnet	Box 1910, Brown U, Prov, RI 02912