# to the n-th power (long)

[Andre Dolenc]

[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 at sauna.hut.fi         SF-02150  Finland
----------------------------------------------------------------------------