ado@sauna.hut.fi (Andre Dolenc) (11/05/90)
[sci.math.num-analysis readers: the function
double pow(double x,double n) = x**n is not always present in C
run-time libraries and we are trying to define one when "n" is an
integer.]
(There are two parallel discussions: the behavior of "pow" and what is
a correct implementation when "n" is an integer.)
------------------------
double
eval_x_n (value, n)
double value;
int n;
/* This algorithm was taken straight from D.E.Knuth "The Art of
* Computer Programming", vol.2, Seminumerical Algorithms,
* section 4.6.3 (pg 398), Addison-Wesley, 1969.
* It is the "binary method", similar to one already posted */
{
double acc, pow_x; /* add "register"s for efficiency */
unsigned int odd, N;
/* a couple, but not all, boundary conditions... */
if ((value == 0.0) && (n <0)) {SignalE (FPE_DIV_ZERO);}
if ((value == 0.0) && (n == 0)) {SignalE (FPE_NaN);}
switch (n)
{
case 0: return 1.0;
case 1: return value;
case 2: return SQUARE(value);
case 3: return SQUARE(value)*value;
case 4: return acc = SQUARE(value), SQUARE(acc);
/* include here other common cases in your application. */
default: break;
}
acc = 1.0, pow_x = value, N = ABS(n); /* wrong if n == -MININT */
LoopForever /* use your favorite macro; "for" or "while" */
{
odd = N&1, N >>= 1;
if (odd) { acc *= pow_x; if (N == 0) return n < 0 ? 1.0/acc : acc;}
pow_x *= pow_x;
}
/*NOTREACHED*/
}
------------------------
Comments:
(1) Before the subject of efficiency can be discussed we must first
determine the correctness. The above algorithm cannot, in my
opinion, be qualified as correct. If not, let's see:
(2) The original problem was to provide a function pow(x,n), with
"n" a double. How can one reliably determine that "n" is an
integer?
I suggest the following:
if (ABS(n - (int) n) < eps)
then /* n is an int */
return eval_x_n (x, (n < 0: (int) n-eps : (int) n+eps));
else
if (x >=0 )
then return (exp(log(x),n)) /* let log(x) handle the 0... */
else /* error */
(3) eval_x_n does not handle floating-point exceptions (FPE).
We must add a signal handler and exit from there using
a longjmp.
(4) As mentioned in the code, if n == MIN_INT, it fails.
(5) The definition of "pow" at the boundaries is not standard.
Some implementations define pow(x,0) = 1, regardless of "x".
(I can post the reasons given for this if there is interest.)
Others return NaN if x is NaN. What now?
(6) Note: under Ultrix the "pow" function is for multiple *integer*
precision arithmetic. It computes x**n mod m, x,n,m are large
integers. Beware!
Regards to you netters,
----------------------------------------------------------------------------
Andre' Dolenc Helsinki University of Technology
Research Associate TKO, Otakaari 1A, Espoo
Email (Internet): ado@sauna.hut.fi SF-02150 Finland
----------------------------------------------------------------------------rctthdb@dutrun.UUCP (H. de Bruijn) (11/23/90)
The function x**n (x and n integer) is also non-present in standard FIG-FORTH. Here's a possible implementation (just for fun): : ^ ( n x -> x**n ) 1 SWAP ROT BEGIN DUP 0 <> WHILE 2 /MOD >R IF SWAP OVER * SWAP THEN DUP * R> REPEAT DROP DROP ; -- * Han de Bruijn; applications engineer | "A little bit of Physics * * TUD Computing Centre; P.O. Box 354 | would be NO idleness in * * 2600 AJ Delft; The Netherlands. | Mathematics" (HdB). * * E-mail: rctthdb@dutrun.tudelft.nl ---| Fax: +31 15 78 37 87 -----*