mingw-w64-crt/math/powil.c

/**
 * This file has no copyright assigned and is placed in the Public Domain.
 * This file is part of the w64 mingw-runtime package.
 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
 */
#include "cephes_mconf.h"

#ifndef _SET_ERRNO
#define _SET_ERRNO(x)
#endif
long double __powil (long double x,int nn);

long double __powil (long double x,int nn)
{
	long double w, y;
	long double s;
	int n, e, sign, asign, lx;

	if (x == 0.0L)
	{
		if (nn == 0)
			return (1.0L);
		else if (nn < 0)
			return (INFINITYL);
		else
			return (0.0L);
	}

	if (nn == 0)
		return (1.0L);

	if (x < 0.0L)
	{
		asign = -1;
		x = -x;
	}
	else
		asign = 0;

	if (nn < 0)
	{
		sign = -1;
		n = -nn;
	}
	else
	{
		sign = 1;
		n = nn;
	}

	/* Overflow detection */

	/* Calculate approximate logarithm of answer */
	s = x;
	s = frexpl(s, &lx);
	e = (lx - 1)*n;
	if ((e == 0) || (e > 64) || (e < -64))
	{
		s = (s - 7.0710678118654752e-1L) / (s +  7.0710678118654752e-1L);
		s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;
	}
	else
	{
		s = LOGE2L * e;
	}

	if (s > MAXLOGL)
	{
		mtherr("__powil", OVERFLOW);
		_SET_ERRNO(ERANGE);
		y = INFINITYL;
		goto done;
	}

	if (s < MINLOGL)
	{
		mtherr("__powil", UNDERFLOW);
		_SET_ERRNO(ERANGE);
		return (0.0L);
	}
	/* Handle tiny denormal answer, but with less accuracy
	 * since roundoff error in 1.0/x will be amplified.
	 * The precise demarcation should be the gradual underflow threshold.
	 */
	if (s < (-MAXLOGL+2.0L))
	{
		x = 1.0L/x;
		sign = -sign;
	}

	/* First bit of the power */
	if (n & 1)
		y = x;
	else
	{
		y = 1.0L;
		asign = 0;
	}

	w = x;
	n >>= 1;
	while (n)
	{
		w = w * w;	/* arg to the 2-to-the-kth power */
		if (n & 1)	/* if that bit is set, then include in product */
			y *= w;
		n >>= 1;
	}

done:

	if (asign)
		y = -y; /* odd power of negative number */
	if (sign < 0)
		y = 1.0L/y;
	return (y);
}