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
----------------------------------------------------------------------------ado@sauna.hut.fi (Andre Dolenc) (11/05/90)
(1) Since my last posting I have learned that the source code to the
math functions (libm) of 4.3 BSD are publicly available and may
be obtained by anon. ftp from several sites, eg.
funic.funet.fi in ~ft/pub/unix/bsd-sources/src/usr.lib/libm
which is a copy of an anon. ftp directory in uunet.uu.net.
Read the copyright notice before using them.
(2) Looking into "pow" one learns that to figure out if "y"
in "pow(x,y)" is an integer is reduced to:
t=drem(x,two);
if ((t == zero) || (copysign(t,one) == one)) then "y is int".
Now, the code which implements "drem" is not trivial. Therefore,
I withdraw my question "How do we decide if y is an integer?".
(3) My previous posting contains (at least) one error. The multiple
precision function "pow" is defined under SunOs in "misc. lib
functions", MP(3X), *not* under Ultrix.
The problem is that the name obviously conflicts with the math
function "double pow(double x, double y)". Isn't dum to have
2 (system) library functions with the same name?
In the hope of not having wasted too much bandwidth,
----------------------------------------------------------------------------
Andre' Dolenc Helsinki University of Technology
Research Associate TKO, Otakaari 1A, Espoo
Email (Internet): ado@sauna.hut.fi SF-02150 Finland
----------------------------------------------------------------------------