201 lines
3.1 KiB
C
201 lines
3.1 KiB
C
|
/* powi.c
|
|||
|
|
*
|
|||
|
|
* Real raised to integer power
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* SYNOPSIS:
|
|||
|
|
*
|
|||
|
|
* double x, y, __powi();
|
|||
|
|
* int n;
|
|||
|
|
*
|
|||
|
|
* y = __powi( x, n );
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* DESCRIPTION:
|
|||
|
|
*
|
|||
|
|
* Returns argument x raised to the nth power.
|
|||
|
|
* The routine efficiently decomposes n as a sum of powers of
|
|||
|
|
* two. The desired power is a product of two-to-the-kth
|
|||
|
|
* powers of x. Thus to compute the 32767 power of x requires
|
|||
|
|
* 28 multiplications instead of 32767 multiplications.
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* ACCURACY:
|
|||
|
|
*
|
|||
|
|
*
|
|||
|
|
* Relative error:
|
|||
|
|
* arithmetic x domain n domain # trials peak rms
|
|||
|
|
* DEC .04,26 -26,26 100000 2.7e-16 4.3e-17
|
|||
|
|
* IEEE .04,26 -26,26 50000 2.0e-15 3.8e-16
|
|||
|
|
* IEEE 1,2 -1022,1023 50000 8.6e-14 1.6e-14
|
|||
|
|
*
|
|||
|
|
* Returns MAXNUM on overflow, zero on underflow.
|
|||
|
|
*
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
/* powi.c */
|
|||
|
|
|
|||
|
|
/*
|
|||
|
|
Cephes Math Library Release 2.8: June, 2000
|
|||
|
|
Copyright 1984, 1995, 2000 by Stephen L. Moshier
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
/*
|
|||
|
|
Modified for mingw
|
|||
|
|
2002-07-22 Danny Smith <dannysmith@users.sourceforge.net>
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
#ifdef __MINGW32__
|
|||
|
|
#include "cephes_mconf.h"
|
|||
|
|
#else
|
|||
|
|
#include "mconf.h"
|
|||
|
|
#ifdef ANSIPROT
|
|||
|
|
extern double log ( double );
|
|||
|
|
extern double frexp ( double, int * );
|
|||
|
|
extern int signbit ( double );
|
|||
|
|
#else
|
|||
|
|
double log(), frexp();
|
|||
|
|
int signbit();
|
|||
|
|
#endif
|
|||
|
|
extern double NEGZERO, INFINITY, MAXNUM, MAXLOG, MINLOG, LOGE2;
|
|||
|
|
#endif /* __MINGW32__ */
|
|||
|
|
|
|||
|
|
#ifndef _SET_ERRNO
|
|||
|
|
#define _SET_ERRNO(x)
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
double __powi( x, nn )
|
|||
|
|
double x;
|
|||
|
|
int nn;
|
|||
|
|
{
|
|||
|
|
int n, e, sign, asign, lx;
|
|||
|
|
double w, y, s;
|
|||
|
|
|
|||
|
|
/* See pow.c for these tests. */
|
|||
|
|
if( x == 0.0 )
|
|||
|
|
{
|
|||
|
|
if( nn == 0 )
|
|||
|
|
return( 1.0 );
|
|||
|
|
else if( nn < 0 )
|
|||
|
|
return( INFINITY );
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
if( nn & 1 )
|
|||
|
|
return( x );
|
|||
|
|
else
|
|||
|
|
return( 0.0 );
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if( nn == 0 )
|
|||
|
|
return( 1.0 );
|
|||
|
|
|
|||
|
|
if( nn == -1 )
|
|||
|
|
return( 1.0/x );
|
|||
|
|
|
|||
|
|
if( x < 0.0 )
|
|||
|
|
{
|
|||
|
|
asign = -1;
|
|||
|
|
x = -x;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
asign = 0;
|
|||
|
|
|
|||
|
|
|
|||
|
|
if( nn < 0 )
|
|||
|
|
{
|
|||
|
|
sign = -1;
|
|||
|
|
n = -nn;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
sign = 1;
|
|||
|
|
n = nn;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Even power will be positive. */
|
|||
|
|
if( (n & 1) == 0 )
|
|||
|
|
asign = 0;
|
|||
|
|
|
|||
|
|
/* Overflow detection */
|
|||
|
|
|
|||
|
|
/* Calculate approximate logarithm of answer */
|
|||
|
|
s = frexp( x, &lx );
|
|||
|
|
e = (lx - 1)*n;
|
|||
|
|
if( (e == 0) || (e > 64) || (e < -64) )
|
|||
|
|
{
|
|||
|
|
s = (s - 7.0710678118654752e-1) / (s + 7.0710678118654752e-1);
|
|||
|
|
s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
s = LOGE2 * e;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if( s > MAXLOG )
|
|||
|
|
{
|
|||
|
|
mtherr( "powi", OVERFLOW );
|
|||
|
|
_SET_ERRNO(ERANGE);
|
|||
|
|
y = INFINITY;
|
|||
|
|
goto done;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
#if DENORMAL
|
|||
|
|
if( s < MINLOG )
|
|||
|
|
{
|
|||
|
|
y = 0.0;
|
|||
|
|
goto done;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* 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 < (-MAXLOG+2.0)) && (sign < 0) )
|
|||
|
|
{
|
|||
|
|
x = 1.0/x;
|
|||
|
|
sign = -sign;
|
|||
|
|
}
|
|||
|
|
#else
|
|||
|
|
/* do not produce denormal answer */
|
|||
|
|
if( s < -MAXLOG )
|
|||
|
|
return(0.0);
|
|||
|
|
#endif
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* First bit of the power */
|
|||
|
|
if( n & 1 )
|
|||
|
|
y = x;
|
|||
|
|
|
|||
|
|
else
|
|||
|
|
y = 1.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;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if( sign < 0 )
|
|||
|
|
y = 1.0/y;
|
|||
|
|
|
|||
|
|
done:
|
|||
|
|
|
|||
|
|
if( asign )
|
|||
|
|
{
|
|||
|
|
/* odd power of negative number */
|
|||
|
|
if( y == 0.0 )
|
|||
|
|
y = NEGZERO;
|
|||
|
|
else
|
|||
|
|
y = -y;
|
|||
|
|
}
|
|||
|
|
return(y);
|
|||
|
|
}
|