ack/lang/cem/libcc.ansi/math/pow.c

/*
 * (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);
}