123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100 |
- /*
- * (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 <math.h>
- #include <float.h>
- #include <errno.h>
- #include <limits.h>
- 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);
- }
|