lib/math/pow.c

   1 /*
   2  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3  * See the copyright notice in the ACK home directory, in the file "Copyright".
   4  *
   5  * Author: Ceriel J.H. Jacobs
   6  */
   7 /* $Header$ */
   8
   9 #include        <math.h>
  10 #include        <float.h>
  11 #include        <errno.h>
  12 #include        "localmath.h"
  13
  14 double
  15 pow(double x, double y)
  16 {
  17         /*      Simple version for now. The Cody and Waite book has
  18                 a very complicated, much more precise version, but
  19                 this version has machine-dependent arrays A1 and A2,
  20                 and I don't know yet how to solve this ???
  21         */
  22         double dummy;
  23         int     result_neg = 0;
  24
  25         if ((x == 0 && y == 0) ||
  26             (x < 0 && modf(y, &dummy) != 0)) {
  27                 errno = EDOM;
  28                 return 0;
  29         }
  30
  31         if (x == 0) return x;
  32
  33         if (x < 0) {
  34                 if (modf(y/2.0, &dummy) != 0) {
  35                         /* y was odd */
  36                         result_neg = 1;
  37                 }
  38                 x = -x;
  39         }
  40         x = log(x);
  41
  42         if (x < 0) {
  43                 x = -x;
  44                 y = -y;
  45         }
  46         /* Beware of overflow in the multiplication */
  47         if (x > 1.0 && y > DBL_MAX/x) {
  48                 errno = ERANGE;
  49                 return result_neg ? -HUGE_VAL : HUGE_VAL;
  50         }
  51
  52         x = exp(x * y);
  53         return result_neg ? -x : x;
  54 }