/* * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. * See the copyright notice in the ACK home directory, in the file "Copyright". * * Author: Ceriel J.H. Jacobs */ /* $Id$ */ #include #include #include #include double pow(double x, double y) { double y_intpart, y_fractpart, fp; int negexp, negx; int ex, newexp; unsigned long yi; if (x == 1.0) return x; if (x == 0 && y <= 0) { errno = EDOM; return 0; } if (y == 0) return 1.0; if (y < 0) { y = -y; negexp = 1; } else negexp = 0; y_fractpart = modf(y, &y_intpart); if (y_fractpart != 0) { if (x < 0) { errno = EDOM; return 0; } } negx = 0; if (x < 0) { x = -x; negx = 1; } if (y_intpart > ULONG_MAX) { if (negx && modf(y_intpart/2.0, &y_fractpart) == 0) { negx = 0; } x = log(x); /* Beware of overflow in the multiplication */ if (x > 1.0 && y > DBL_MAX/x) { errno = ERANGE; return HUGE_VAL; } if (negexp) y = -y; if (negx) return -exp(x*y); return exp(x * y); } if (y_fractpart != 0) { fp = exp(y_fractpart * log(x)); } else fp = 1.0; yi = y_intpart; if (! (yi & 1)) negx = 0; x = frexp(x, &ex); newexp = 0; for (;;) { if (yi & 1) { fp *= x; newexp += ex; } yi >>= 1; if (yi == 0) break; x *= x; ex <<= 1; if (x < 0.5) { x += x; ex -= 1; } } if (negexp) { fp = 1.0/fp; newexp = -newexp; } if (negx) { return -ldexp(fp, newexp); } return ldexp(fp, newexp); }