mirror of
git://sourceware.org/git/newlib-cygwin.git
synced 2025-02-20 16:01:10 +08:00
* mingwex/math/powil.c: Rename powil to __powil.
* mingwex/math/powl.c: Adjust declaration and call to __powil. Remove comment on powil. * mingwex/math/powi.c: New file. * mingwex/math/powif.c: New file. * mingwex/math/pow.c: New file. * mingwex/math/cephes_mconf.h. Add double and float versions of constants. (polevl): Add double precision function. (p1evl): Likewise. * mingwex/Makefile.in (MATH_DISTFILES): Add pow.c, powi.c, powif.c. (MATH_OBJS): Add pow.o, powi.o, powif.o.
This commit is contained in:
Danny Smith
parent
2bacbfb1d1
commit
66451d9590
@ -1,3 +1,19 @@
|
||||
2002-10-07 Danny Smith <dannysmith@users.sourceforge.net>
|
||||
|
||||
* mingwex/math/powil.c: Rename powil to __powil.
|
||||
* mingwex/math/powl.c: Adjust declaration and call
|
||||
to __powil. Remove comment on powil.
|
||||
* mingwex/math/powi.c: New file.
|
||||
* mingwex/math/powif.c: New file.
|
||||
* mingwex/math/pow.c: New file.
|
||||
* mingwex/math/cephes_mconf.h. Add double and float
|
||||
versions of constants.
|
||||
(polevl): Add double precision function.
|
||||
(p1evl): Likewise.
|
||||
* mingwex/Makefile.in (MATH_DISTFILES): Add pow.c,
|
||||
powi.c, powif.c.
|
||||
(MATH_OBJS): Add pow.o, powi.o, powif.o.
|
||||
|
||||
2002-10-03 Danny Smith <dannysmith@users.sourceforge.net>
|
||||
|
||||
* include/cytpe.h (_imp____mbcur_max): Add missing ';'.
|
||||
|
||||
@ -50,8 +50,9 @@ MATH_DISTFILES = \
|
||||
log10f.S log10l.S log1p.S log1pf.S log1pl.S log2.S log2f.S \
|
||||
log2l.S logb.c logbf.c logbl.c logf.S logl.S lrint.c lrintf.c \
|
||||
lrintl.c lround.c lroundf.c lroundl.c modff.c modfl.c \
|
||||
nearbyint.S nearbyintf.S nearbyintl.S nextafterf.c powf.c powil.c \
|
||||
powl.c remainder.S remainderf.S remainderl.S remquo.S \
|
||||
nearbyint.S nearbyintf.S nearbyintl.S nextafterf.c \
|
||||
pow.c powf.c powi.c powif.c powil.c powl.c \
|
||||
remainder.S remainderf.S remainderl.S remquo.S \
|
||||
remquof.S remquol.S rint.c rintf.c rintl.c round.c roundf.c \
|
||||
roundl.c scalbn.S scalbnf.S scalbnl.S signbit.c signbitf.c \
|
||||
signbitl.c sinf.S sinhf.c sinhl.c sinl.S sqrtf.c sqrtl.c \
|
||||
@ -115,8 +116,9 @@ MATH_OBJS = \
|
||||
log10f.o log10l.o log1p.o log1pf.o log1pl.o log2.o log2f.o \
|
||||
log2l.o logb.o logbf.o logbl.o logf.o logl.o lrint.o lrintf.o \
|
||||
lrintl.o lround.o lroundf.o lroundl.o modff.o modfl.o \
|
||||
nearbyint.o nearbyintf.o nearbyintl.o nextafterf.o powf.o powil.o \
|
||||
powl.o remainder.o remainderf.o remainderl.o remquo.o \
|
||||
nearbyint.o nearbyintf.o nearbyintl.o nextafterf.o \
|
||||
pow.o powf.o powi.o powif.o powil.o powl.o \
|
||||
remainder.o remainderf.o remainderl.o remquo.o \
|
||||
remquof.o remquol.o rint.o rintf.o rintl.o round.o roundf.o \
|
||||
roundl.o scalbn.o scalbnf.o scalbnl.o signbit.o signbitf.o \
|
||||
signbitl.o sinf.o sinhf.o sinhl.o sinl.o sqrtf.o sqrtl.o \
|
||||
|
||||
@ -1,8 +1,45 @@
|
||||
#include <math.h>
|
||||
#include <errno.h>
|
||||
|
||||
|
||||
#define IBMPC 1
|
||||
#define ANSIPROT 1
|
||||
#define MINUSZERO 1
|
||||
#define INFINITIES 1
|
||||
#define NANS 1
|
||||
#define DENORMAL 1
|
||||
#define VOLATILE
|
||||
#define mtherr(fname, code)
|
||||
#define XPD 0,
|
||||
|
||||
#ifdef _CEPHES_USE_ERRNO
|
||||
#define _SET_ERRNO(x) errno = (x)
|
||||
#else
|
||||
#define _SET_ERRNO(x)
|
||||
#endif
|
||||
|
||||
/* constants used by cephes functions */
|
||||
|
||||
/* double */
|
||||
#define MAXNUM 1.7976931348623158E308
|
||||
#define MAXLOG 7.09782712893383996843E2
|
||||
#define MINLOG -7.08396418532264106224E2
|
||||
#define LOGE2 6.93147180559945309417E-1
|
||||
#define LOG2E 1.44269504088896340736
|
||||
#define PI 3.14159265358979323846
|
||||
#define PIO2 1.57079632679489661923
|
||||
#define PIO4 7.85398163397448309616E-1
|
||||
|
||||
#define NEGZERO (-0.0)
|
||||
extern double __INF;
|
||||
#undef INFINITY
|
||||
#define INFINITY (__INF)
|
||||
extern double __QNAN;
|
||||
#undef NAN
|
||||
#define NAN (__QNAN)
|
||||
|
||||
|
||||
/*long double*/
|
||||
#define MAXNUML 1.189731495357231765021263853E4932L
|
||||
#define MAXLOGL 1.1356523406294143949492E4L
|
||||
#define MINLOGL -1.13994985314888605586758E4L
|
||||
@ -17,27 +54,133 @@
|
||||
#define isnanl isnan
|
||||
#define signbitl signbit
|
||||
|
||||
|
||||
#define IBMPC 1
|
||||
#define ANSIPROT 1
|
||||
#define MINUSZERO 1
|
||||
#define INFINITIES 1
|
||||
#define NANS 1
|
||||
#define DENORMAL 1
|
||||
#define NEGZEROL (-0.0L)
|
||||
extern long double __INFL;
|
||||
#define INFINITYL (__INFL)
|
||||
extern long double __QNANL;
|
||||
#define NANL (__QNANL)
|
||||
#define VOLATILE
|
||||
#define mtherr(fname, code)
|
||||
#define XPD 0,
|
||||
|
||||
#ifdef _CEPHES_USE_ERRNO
|
||||
#define _SET_ERRNO(x) errno = (x)
|
||||
#else
|
||||
#define _SET_ERRNO(x)
|
||||
#endif
|
||||
/* float */
|
||||
|
||||
#define MAXNUMF 3.4028234663852885981170418348451692544e38F
|
||||
#define MAXLOGF 88.72283905206835F
|
||||
#define MINLOGF -103.278929903431851103F /* log(2^-149) */
|
||||
#define LOG2EF 1.44269504088896341F
|
||||
#define LOGE2F 0.693147180559945309F
|
||||
#define PIF 3.141592653589793238F
|
||||
#define PIO2F 1.5707963267948966192F
|
||||
#define PIO4F 0.7853981633974483096F
|
||||
|
||||
#define isfinitef isfinite
|
||||
#define isinff isinf
|
||||
#define isnanf isnan
|
||||
#define signbitf signbit
|
||||
|
||||
#define NEGZEROF (-0.0F)
|
||||
extern float __INFF;
|
||||
#define INFINITYF (__INFF)
|
||||
extern float __QNANF;
|
||||
#define NANF (__QNANF)
|
||||
|
||||
|
||||
/* double */
|
||||
|
||||
/*
|
||||
Cephes Math Library Release 2.2: July, 1992
|
||||
Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
|
||||
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
|
||||
*/
|
||||
|
||||
|
||||
/* polevl.c
|
||||
* p1evl.c
|
||||
*
|
||||
* Evaluate polynomial
|
||||
*
|
||||
*
|
||||
*
|
||||
* SYNOPSIS:
|
||||
*
|
||||
* int N;
|
||||
* double x, y, coef[N+1], polevl[];
|
||||
*
|
||||
* y = polevl( x, coef, N );
|
||||
*
|
||||
*
|
||||
*
|
||||
* DESCRIPTION:
|
||||
*
|
||||
* Evaluates polynomial of degree N:
|
||||
*
|
||||
* 2 N
|
||||
* y = C + C x + C x +...+ C x
|
||||
* 0 1 2 N
|
||||
*
|
||||
* Coefficients are stored in reverse order:
|
||||
*
|
||||
* coef[0] = C , ..., coef[N] = C .
|
||||
* N 0
|
||||
*
|
||||
* The function p1evl() assumes that coef[N] = 1.0 and is
|
||||
* omitted from the array. Its calling arguments are
|
||||
* otherwise the same as polevl().
|
||||
*
|
||||
*
|
||||
* SPEED:
|
||||
*
|
||||
* In the interest of speed, there are no checks for out
|
||||
* of bounds arithmetic. This routine is used by most of
|
||||
* the functions in the library. Depending on available
|
||||
* equipment features, the user may wish to rewrite the
|
||||
* program in microcode or assembly language.
|
||||
*
|
||||
*/
|
||||
|
||||
/* Polynomial evaluator:
|
||||
* P[0] x^n + P[1] x^(n-1) + ... + P[n]
|
||||
*/
|
||||
static __inline__ double polevl( x, p, n )
|
||||
double x;
|
||||
const void *p;
|
||||
int n;
|
||||
{
|
||||
register double y;
|
||||
register double *P = (double *)p;
|
||||
|
||||
y = *P++;
|
||||
do
|
||||
{
|
||||
y = y * x + *P++;
|
||||
}
|
||||
while( --n );
|
||||
return(y);
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* Polynomial evaluator:
|
||||
* x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n]
|
||||
*/
|
||||
static __inline__ double p1evl( x, p, n )
|
||||
double x;
|
||||
const void *p;
|
||||
int n;
|
||||
{
|
||||
register double y;
|
||||
register double *P = (double *)p;
|
||||
|
||||
n -= 1;
|
||||
y = x + *P++;
|
||||
do
|
||||
{
|
||||
y = y * x + *P++;
|
||||
}
|
||||
while( --n );
|
||||
return( y );
|
||||
}
|
||||
|
||||
|
||||
/* long double */
|
||||
/*
|
||||
Cephes Math Library Release 2.2: July, 1992
|
||||
Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
|
||||
|
||||
781
winsup/mingw/mingwex/math/pow.c
Normal file
781
winsup/mingw/mingwex/math/pow.c
Normal file
@ -0,0 +1,781 @@
|
||||
/* pow.c
|
||||
*
|
||||
* Power function
|
||||
*
|
||||
*
|
||||
*
|
||||
* SYNOPSIS:
|
||||
*
|
||||
* double x, y, z, pow();
|
||||
*
|
||||
* z = pow( x, y );
|
||||
*
|
||||
*
|
||||
*
|
||||
* DESCRIPTION:
|
||||
*
|
||||
* Computes x raised to the yth power. Analytically,
|
||||
*
|
||||
* x**y = exp( y log(x) ).
|
||||
*
|
||||
* Following Cody and Waite, this program uses a lookup table
|
||||
* of 2**-i/16 and pseudo extended precision arithmetic to
|
||||
* obtain an extra three bits of accuracy in both the logarithm
|
||||
* and the exponential.
|
||||
*
|
||||
*
|
||||
*
|
||||
* ACCURACY:
|
||||
*
|
||||
* Relative error:
|
||||
* arithmetic domain # trials peak rms
|
||||
* IEEE -26,26 30000 4.2e-16 7.7e-17
|
||||
* DEC -26,26 60000 4.8e-17 9.1e-18
|
||||
* 1/26 < x < 26, with log(x) uniformly distributed.
|
||||
* -26 < y < 26, y uniformly distributed.
|
||||
* IEEE 0,8700 30000 1.5e-14 2.1e-15
|
||||
* 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
|
||||
*
|
||||
*
|
||||
* ERROR MESSAGES:
|
||||
*
|
||||
* message condition value returned
|
||||
* pow overflow x**y > MAXNUM INFINITY
|
||||
* pow underflow x**y < 1/MAXNUM 0.0
|
||||
* pow domain x<0 and y noninteger 0.0
|
||||
*
|
||||
*/
|
||||
|
||||
/*
|
||||
Cephes Math Library Release 2.8: June, 2000
|
||||
Copyright 1984, 1995, 2000 by Stephen L. Moshier
|
||||
*/
|
||||
|
||||
/*
|
||||
Modified for mingw
|
||||
2002-09-27 Danny Smith <dannysmith@users.sourceforge.net>
|
||||
*/
|
||||
|
||||
#ifdef __MINGW32__
|
||||
#include "cephes_mconf.h"
|
||||
#else
|
||||
#include "mconf.h"
|
||||
static char fname[] = {"pow"};
|
||||
#endif
|
||||
|
||||
#ifndef _SET_ERRNO
|
||||
#define _SET_ERRNO(x)
|
||||
#endif
|
||||
|
||||
#define SQRTH 0.70710678118654752440
|
||||
|
||||
#ifdef UNK
|
||||
static double P[] = {
|
||||
4.97778295871696322025E-1,
|
||||
3.73336776063286838734E0,
|
||||
7.69994162726912503298E0,
|
||||
4.66651806774358464979E0
|
||||
};
|
||||
static double Q[] = {
|
||||
/* 1.00000000000000000000E0, */
|
||||
9.33340916416696166113E0,
|
||||
2.79999886606328401649E1,
|
||||
3.35994905342304405431E1,
|
||||
1.39995542032307539578E1
|
||||
};
|
||||
/* 2^(-i/16), IEEE precision */
|
||||
static double A[] = {
|
||||
1.00000000000000000000E0,
|
||||
9.57603280698573700036E-1,
|
||||
9.17004043204671215328E-1,
|
||||
8.78126080186649726755E-1,
|
||||
8.40896415253714502036E-1,
|
||||
8.05245165974627141736E-1,
|
||||
7.71105412703970372057E-1,
|
||||
7.38413072969749673113E-1,
|
||||
7.07106781186547572737E-1,
|
||||
6.77127773468446325644E-1,
|
||||
6.48419777325504820276E-1,
|
||||
6.20928906036742001007E-1,
|
||||
5.94603557501360513449E-1,
|
||||
5.69394317378345782288E-1,
|
||||
5.45253866332628844837E-1,
|
||||
5.22136891213706877402E-1,
|
||||
5.00000000000000000000E-1
|
||||
};
|
||||
static double B[] = {
|
||||
0.00000000000000000000E0,
|
||||
1.64155361212281360176E-17,
|
||||
4.09950501029074826006E-17,
|
||||
3.97491740484881042808E-17,
|
||||
-4.83364665672645672553E-17,
|
||||
1.26912513974441574796E-17,
|
||||
1.99100761573282305549E-17,
|
||||
-1.52339103990623557348E-17,
|
||||
0.00000000000000000000E0
|
||||
};
|
||||
static double R[] = {
|
||||
1.49664108433729301083E-5,
|
||||
1.54010762792771901396E-4,
|
||||
1.33335476964097721140E-3,
|
||||
9.61812908476554225149E-3,
|
||||
5.55041086645832347466E-2,
|
||||
2.40226506959099779976E-1,
|
||||
6.93147180559945308821E-1
|
||||
};
|
||||
|
||||
#define douba(k) A[k]
|
||||
#define doubb(k) B[k]
|
||||
#define MEXP 16383.0
|
||||
#ifdef DENORMAL
|
||||
#define MNEXP -17183.0
|
||||
#else
|
||||
#define MNEXP -16383.0
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#ifdef DEC
|
||||
static unsigned short P[] = {
|
||||
0037776,0156313,0175332,0163602,
|
||||
0040556,0167577,0052366,0174245,
|
||||
0040766,0062753,0175707,0055564,
|
||||
0040625,0052035,0131344,0155636,
|
||||
};
|
||||
static unsigned short Q[] = {
|
||||
/*0040200,0000000,0000000,0000000,*/
|
||||
0041025,0052644,0154404,0105155,
|
||||
0041337,0177772,0007016,0047646,
|
||||
0041406,0062740,0154273,0020020,
|
||||
0041137,0177054,0106127,0044555,
|
||||
};
|
||||
static unsigned short A[] = {
|
||||
0040200,0000000,0000000,0000000,
|
||||
0040165,0022575,0012444,0103314,
|
||||
0040152,0140306,0163735,0022071,
|
||||
0040140,0146336,0166052,0112341,
|
||||
0040127,0042374,0145326,0116553,
|
||||
0040116,0022214,0012437,0102201,
|
||||
0040105,0063452,0010525,0003333,
|
||||
0040075,0004243,0117530,0006067,
|
||||
0040065,0002363,0031771,0157145,
|
||||
0040055,0054076,0165102,0120513,
|
||||
0040045,0177326,0124661,0050471,
|
||||
0040036,0172462,0060221,0120422,
|
||||
0040030,0033760,0050615,0134251,
|
||||
0040021,0141723,0071653,0010703,
|
||||
0040013,0112701,0161752,0105727,
|
||||
0040005,0125303,0063714,0044173,
|
||||
0040000,0000000,0000000,0000000
|
||||
};
|
||||
static unsigned short B[] = {
|
||||
0000000,0000000,0000000,0000000,
|
||||
0021473,0040265,0153315,0140671,
|
||||
0121074,0062627,0042146,0176454,
|
||||
0121413,0003524,0136332,0066212,
|
||||
0121767,0046404,0166231,0012553,
|
||||
0121257,0015024,0002357,0043574,
|
||||
0021736,0106532,0043060,0056206,
|
||||
0121310,0020334,0165705,0035326,
|
||||
0000000,0000000,0000000,0000000
|
||||
};
|
||||
|
||||
static unsigned short R[] = {
|
||||
0034173,0014076,0137624,0115771,
|
||||
0035041,0076763,0003744,0111311,
|
||||
0035656,0141766,0041127,0074351,
|
||||
0036435,0112533,0073611,0116664,
|
||||
0037143,0054106,0134040,0152223,
|
||||
0037565,0176757,0176026,0025551,
|
||||
0040061,0071027,0173721,0147572
|
||||
};
|
||||
|
||||
/*
|
||||
static double R[] = {
|
||||
0.14928852680595608186e-4,
|
||||
0.15400290440989764601e-3,
|
||||
0.13333541313585784703e-2,
|
||||
0.96181290595172416964e-2,
|
||||
0.55504108664085595326e-1,
|
||||
0.24022650695909537056e0,
|
||||
0.69314718055994529629e0
|
||||
};
|
||||
*/
|
||||
#define douba(k) (*(double *)&A[(k)<<2])
|
||||
#define doubb(k) (*(double *)&B[(k)<<2])
|
||||
#define MEXP 2031.0
|
||||
#define MNEXP -2031.0
|
||||
#endif
|
||||
|
||||
#ifdef IBMPC
|
||||
static const unsigned short P[] = {
|
||||
0x5cf0,0x7f5b,0xdb99,0x3fdf,
|
||||
0xdf15,0xea9e,0xddef,0x400d,
|
||||
0xeb6f,0x7f78,0xccbd,0x401e,
|
||||
0x9b74,0xb65c,0xaa83,0x4012,
|
||||
};
|
||||
static const unsigned short Q[] = {
|
||||
/*0x0000,0x0000,0x0000,0x3ff0,*/
|
||||
0x914e,0x9b20,0xaab4,0x4022,
|
||||
0xc9f5,0x41c1,0xffff,0x403b,
|
||||
0x6402,0x1b17,0xccbc,0x4040,
|
||||
0xe92e,0x918a,0xffc5,0x402b,
|
||||
};
|
||||
static const unsigned short A[] = {
|
||||
0x0000,0x0000,0x0000,0x3ff0,
|
||||
0x90da,0xa2a4,0xa4af,0x3fee,
|
||||
0xa487,0xdcfb,0x5818,0x3fed,
|
||||
0x529c,0xdd85,0x199b,0x3fec,
|
||||
0xd3ad,0x995a,0xe89f,0x3fea,
|
||||
0xf090,0x82a3,0xc491,0x3fe9,
|
||||
0xa0db,0x422a,0xace5,0x3fe8,
|
||||
0x0187,0x73eb,0xa114,0x3fe7,
|
||||
0x3bcd,0x667f,0xa09e,0x3fe6,
|
||||
0x5429,0xdd48,0xab07,0x3fe5,
|
||||
0x2a27,0xd536,0xbfda,0x3fe4,
|
||||
0x3422,0x4c12,0xdea6,0x3fe3,
|
||||
0xb715,0x0a31,0x06fe,0x3fe3,
|
||||
0x6238,0x6e75,0x387a,0x3fe2,
|
||||
0x517b,0x3c7d,0x72b8,0x3fe1,
|
||||
0x890f,0x6cf9,0xb558,0x3fe0,
|
||||
0x0000,0x0000,0x0000,0x3fe0
|
||||
};
|
||||
static const unsigned short B[] = {
|
||||
0x0000,0x0000,0x0000,0x0000,
|
||||
0x3707,0xd75b,0xed02,0x3c72,
|
||||
0xcc81,0x345d,0xa1cd,0x3c87,
|
||||
0x4b27,0x5686,0xe9f1,0x3c86,
|
||||
0x6456,0x13b2,0xdd34,0xbc8b,
|
||||
0x42e2,0xafec,0x4397,0x3c6d,
|
||||
0x82e4,0xd231,0xf46a,0x3c76,
|
||||
0x8a76,0xb9d7,0x9041,0xbc71,
|
||||
0x0000,0x0000,0x0000,0x0000
|
||||
};
|
||||
static const unsigned short R[] = {
|
||||
0x937f,0xd7f2,0x6307,0x3eef,
|
||||
0x9259,0x60fc,0x2fbe,0x3f24,
|
||||
0xef1d,0xc84a,0xd87e,0x3f55,
|
||||
0x33b7,0x6ef1,0xb2ab,0x3f83,
|
||||
0x1a92,0xd704,0x6b08,0x3fac,
|
||||
0xc56d,0xff82,0xbfbd,0x3fce,
|
||||
0x39ef,0xfefa,0x2e42,0x3fe6
|
||||
};
|
||||
|
||||
#define douba(k) (*(double *)&A[(k)<<2])
|
||||
#define doubb(k) (*(double *)&B[(k)<<2])
|
||||
#define MEXP 16383.0
|
||||
#ifdef DENORMAL
|
||||
#define MNEXP -17183.0
|
||||
#else
|
||||
#define MNEXP -16383.0
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#ifdef MIEEE
|
||||
static unsigned short P[] = {
|
||||
0x3fdf,0xdb99,0x7f5b,0x5cf0,
|
||||
0x400d,0xddef,0xea9e,0xdf15,
|
||||
0x401e,0xccbd,0x7f78,0xeb6f,
|
||||
0x4012,0xaa83,0xb65c,0x9b74
|
||||
};
|
||||
static unsigned short Q[] = {
|
||||
0x4022,0xaab4,0x9b20,0x914e,
|
||||
0x403b,0xffff,0x41c1,0xc9f5,
|
||||
0x4040,0xccbc,0x1b17,0x6402,
|
||||
0x402b,0xffc5,0x918a,0xe92e
|
||||
};
|
||||
static unsigned short A[] = {
|
||||
0x3ff0,0x0000,0x0000,0x0000,
|
||||
0x3fee,0xa4af,0xa2a4,0x90da,
|
||||
0x3fed,0x5818,0xdcfb,0xa487,
|
||||
0x3fec,0x199b,0xdd85,0x529c,
|
||||
0x3fea,0xe89f,0x995a,0xd3ad,
|
||||
0x3fe9,0xc491,0x82a3,0xf090,
|
||||
0x3fe8,0xace5,0x422a,0xa0db,
|
||||
0x3fe7,0xa114,0x73eb,0x0187,
|
||||
0x3fe6,0xa09e,0x667f,0x3bcd,
|
||||
0x3fe5,0xab07,0xdd48,0x5429,
|
||||
0x3fe4,0xbfda,0xd536,0x2a27,
|
||||
0x3fe3,0xdea6,0x4c12,0x3422,
|
||||
0x3fe3,0x06fe,0x0a31,0xb715,
|
||||
0x3fe2,0x387a,0x6e75,0x6238,
|
||||
0x3fe1,0x72b8,0x3c7d,0x517b,
|
||||
0x3fe0,0xb558,0x6cf9,0x890f,
|
||||
0x3fe0,0x0000,0x0000,0x0000
|
||||
};
|
||||
static unsigned short B[] = {
|
||||
0x0000,0x0000,0x0000,0x0000,
|
||||
0x3c72,0xed02,0xd75b,0x3707,
|
||||
0x3c87,0xa1cd,0x345d,0xcc81,
|
||||
0x3c86,0xe9f1,0x5686,0x4b27,
|
||||
0xbc8b,0xdd34,0x13b2,0x6456,
|
||||
0x3c6d,0x4397,0xafec,0x42e2,
|
||||
0x3c76,0xf46a,0xd231,0x82e4,
|
||||
0xbc71,0x9041,0xb9d7,0x8a76,
|
||||
0x0000,0x0000,0x0000,0x0000
|
||||
};
|
||||
static unsigned short R[] = {
|
||||
0x3eef,0x6307,0xd7f2,0x937f,
|
||||
0x3f24,0x2fbe,0x60fc,0x9259,
|
||||
0x3f55,0xd87e,0xc84a,0xef1d,
|
||||
0x3f83,0xb2ab,0x6ef1,0x33b7,
|
||||
0x3fac,0x6b08,0xd704,0x1a92,
|
||||
0x3fce,0xbfbd,0xff82,0xc56d,
|
||||
0x3fe6,0x2e42,0xfefa,0x39ef
|
||||
};
|
||||
|
||||
#define douba(k) (*(double *)&A[(k)<<2])
|
||||
#define doubb(k) (*(double *)&B[(k)<<2])
|
||||
#define MEXP 16383.0
|
||||
#ifdef DENORMAL
|
||||
#define MNEXP -17183.0
|
||||
#else
|
||||
#define MNEXP -16383.0
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* log2(e) - 1 */
|
||||
#define LOG2EA 0.44269504088896340736
|
||||
|
||||
#define F W
|
||||
#define Fa Wa
|
||||
#define Fb Wb
|
||||
#define G W
|
||||
#define Ga Wa
|
||||
#define Gb u
|
||||
#define H W
|
||||
#define Ha Wb
|
||||
#define Hb Wb
|
||||
|
||||
#ifdef __MINGW32__
|
||||
static __inline__ double reduc( double );
|
||||
extern double __powi ( double, int );
|
||||
extern double pow ( double x, double y);
|
||||
|
||||
#else /* __MINGW32__ */
|
||||
|
||||
#ifdef ANSIPROT
|
||||
extern double floor ( double );
|
||||
extern double fabs ( double );
|
||||
extern double frexp ( double, int * );
|
||||
extern double ldexp ( double, int );
|
||||
extern double polevl ( double, void *, int );
|
||||
extern double p1evl ( double, void *, int );
|
||||
extern double __powi ( double, int );
|
||||
extern int signbit ( double );
|
||||
extern int isnan ( double );
|
||||
extern int isfinite ( double );
|
||||
static double reduc ( double );
|
||||
#else
|
||||
double floor(), fabs(), frexp(), ldexp();
|
||||
double polevl(), p1evl(), __powi();
|
||||
int signbit(), isnan(), isfinite();
|
||||
static double reduc();
|
||||
#endif
|
||||
extern double MAXNUM;
|
||||
#ifdef INFINITIES
|
||||
extern double INFINITY;
|
||||
#endif
|
||||
#ifdef NANS
|
||||
extern double NAN;
|
||||
#endif
|
||||
#ifdef MINUSZERO
|
||||
extern double NEGZERO;
|
||||
#endif
|
||||
|
||||
#endif /* __MINGW32__ */
|
||||
|
||||
double pow( x, y )
|
||||
double x, y;
|
||||
{
|
||||
double w, z, W, Wa, Wb, ya, yb, u;
|
||||
/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
|
||||
double aw, ay, wy;
|
||||
int e, i, nflg, iyflg, yoddint;
|
||||
|
||||
if( y == 0.0 )
|
||||
return( 1.0 );
|
||||
#ifdef NANS
|
||||
if( isnan(x) || isnan(y) )
|
||||
{
|
||||
_SET_ERRNO (EDOM);
|
||||
return( x + y );
|
||||
}
|
||||
#endif
|
||||
if( y == 1.0 )
|
||||
return( x );
|
||||
|
||||
|
||||
#ifdef INFINITIES
|
||||
if( !isfinite(y) && (x == 1.0 || x == -1.0) )
|
||||
{
|
||||
mtherr( "pow", DOMAIN );
|
||||
#ifdef NANS
|
||||
return( NAN );
|
||||
#else
|
||||
return( INFINITY );
|
||||
#endif
|
||||
}
|
||||
#endif
|
||||
|
||||
if( x == 1.0 )
|
||||
return( 1.0 );
|
||||
|
||||
if( y >= MAXNUM )
|
||||
{
|
||||
_SET_ERRNO (ERANGE);
|
||||
#ifdef INFINITIES
|
||||
if( x > 1.0 )
|
||||
return( INFINITY );
|
||||
#else
|
||||
if( x > 1.0 )
|
||||
return( MAXNUM );
|
||||
#endif
|
||||
if( x > 0.0 && x < 1.0 )
|
||||
return( 0.0);
|
||||
if( x < -1.0 )
|
||||
{
|
||||
#ifdef INFINITIES
|
||||
return( INFINITY );
|
||||
#else
|
||||
return( MAXNUM );
|
||||
#endif
|
||||
}
|
||||
if( x > -1.0 && x < 0.0 )
|
||||
return( 0.0 );
|
||||
}
|
||||
if( y <= -MAXNUM )
|
||||
{
|
||||
_SET_ERRNO (ERANGE);
|
||||
if( x > 1.0 )
|
||||
return( 0.0 );
|
||||
#ifdef INFINITIES
|
||||
if( x > 0.0 && x < 1.0 )
|
||||
return( INFINITY );
|
||||
#else
|
||||
if( x > 0.0 && x < 1.0 )
|
||||
return( MAXNUM );
|
||||
#endif
|
||||
if( x < -1.0 )
|
||||
return( 0.0 );
|
||||
#ifdef INFINITIES
|
||||
if( x > -1.0 && x < 0.0 )
|
||||
return( INFINITY );
|
||||
#else
|
||||
if( x > -1.0 && x < 0.0 )
|
||||
return( MAXNUM );
|
||||
#endif
|
||||
}
|
||||
if( x >= MAXNUM )
|
||||
{
|
||||
#if INFINITIES
|
||||
if( y > 0.0 )
|
||||
return( INFINITY );
|
||||
#else
|
||||
if( y > 0.0 )
|
||||
return( MAXNUM );
|
||||
#endif
|
||||
return(0.0);
|
||||
}
|
||||
/* Set iyflg to 1 if y is an integer. */
|
||||
iyflg = 0;
|
||||
w = floor(y);
|
||||
if( w == y )
|
||||
iyflg = 1;
|
||||
|
||||
/* Test for odd integer y. */
|
||||
yoddint = 0;
|
||||
if( iyflg )
|
||||
{
|
||||
ya = fabs(y);
|
||||
ya = floor(0.5 * ya);
|
||||
yb = 0.5 * fabs(w);
|
||||
if( ya != yb )
|
||||
yoddint = 1;
|
||||
}
|
||||
|
||||
if( x <= -MAXNUM )
|
||||
{
|
||||
if( y > 0.0 )
|
||||
{
|
||||
#ifdef INFINITIES
|
||||
if( yoddint )
|
||||
return( -INFINITY );
|
||||
return( INFINITY );
|
||||
#else
|
||||
if( yoddint )
|
||||
return( -MAXNUM );
|
||||
return( MAXNUM );
|
||||
#endif
|
||||
}
|
||||
if( y < 0.0 )
|
||||
{
|
||||
#ifdef MINUSZERO
|
||||
if( yoddint )
|
||||
return( NEGZERO );
|
||||
#endif
|
||||
return( 0.0 );
|
||||
}
|
||||
}
|
||||
|
||||
nflg = 0; /* flag = 1 if x<0 raised to integer power */
|
||||
if( x <= 0.0 )
|
||||
{
|
||||
if( x == 0.0 )
|
||||
{
|
||||
if( y < 0.0 )
|
||||
{
|
||||
#ifdef MINUSZERO
|
||||
if( signbit(x) && yoddint )
|
||||
return( -INFINITY );
|
||||
#endif
|
||||
#ifdef INFINITIES
|
||||
return( INFINITY );
|
||||
#else
|
||||
return( MAXNUM );
|
||||
#endif
|
||||
}
|
||||
if( y > 0.0 )
|
||||
{
|
||||
#ifdef MINUSZERO
|
||||
if( signbit(x) && yoddint )
|
||||
return( NEGZERO );
|
||||
#endif
|
||||
return( 0.0 );
|
||||
}
|
||||
return( 1.0 );
|
||||
}
|
||||
else
|
||||
{
|
||||
if( iyflg == 0 )
|
||||
{ /* noninteger power of negative number */
|
||||
mtherr( fname, DOMAIN );
|
||||
_SET_ERRNO (EDOM);
|
||||
#ifdef NANS
|
||||
return(NAN);
|
||||
#else
|
||||
return(0.0L);
|
||||
#endif
|
||||
}
|
||||
nflg = 1;
|
||||
}
|
||||
}
|
||||
|
||||
/* Integer power of an integer. */
|
||||
|
||||
if( iyflg )
|
||||
{
|
||||
i = w;
|
||||
w = floor(x);
|
||||
if( (w == x) && (fabs(y) < 32768.0) )
|
||||
{
|
||||
w = __powi( x, (int) y );
|
||||
return( w );
|
||||
}
|
||||
}
|
||||
|
||||
if( nflg )
|
||||
x = fabs(x);
|
||||
|
||||
/* For results close to 1, use a series expansion. */
|
||||
w = x - 1.0;
|
||||
aw = fabs(w);
|
||||
ay = fabs(y);
|
||||
wy = w * y;
|
||||
ya = fabs(wy);
|
||||
if((aw <= 1.0e-3 && ay <= 1.0)
|
||||
|| (ya <= 1.0e-3 && ay >= 1.0))
|
||||
{
|
||||
z = (((((w*(y-5.)/720. + 1./120.)*w*(y-4.) + 1./24.)*w*(y-3.)
|
||||
+ 1./6.)*w*(y-2.) + 0.5)*w*(y-1.) )*wy + wy + 1.;
|
||||
goto done;
|
||||
}
|
||||
/* These are probably too much trouble. */
|
||||
#if 0
|
||||
w = y * log(x);
|
||||
if (aw > 1.0e-3 && fabs(w) < 1.0e-3)
|
||||
{
|
||||
z = ((((((
|
||||
w/7. + 1.)*w/6. + 1.)*w/5. + 1.)*w/4. + 1.)*w/3. + 1.)*w/2. + 1.)*w + 1.;
|
||||
goto done;
|
||||
}
|
||||
|
||||
if(ya <= 1.0e-3 && aw <= 1.0e-4)
|
||||
{
|
||||
z = (((((
|
||||
wy*1./720.
|
||||
+ (-w*1./48. + 1./120.) )*wy
|
||||
+ ((w*17./144. - 1./12.)*w + 1./24.) )*wy
|
||||
+ (((-w*5./16. + 7./24.)*w - 1./4.)*w + 1./6.) )*wy
|
||||
+ ((((w*137./360. - 5./12.)*w + 11./24.)*w - 1./2.)*w + 1./2.) )*wy
|
||||
+ (((((-w*1./6. + 1./5.)*w - 1./4)*w + 1./3.)*w -1./2.)*w ) )*wy
|
||||
+ wy + 1.0;
|
||||
goto done;
|
||||
}
|
||||
#endif
|
||||
|
||||
/* separate significand from exponent */
|
||||
x = frexp( x, &e );
|
||||
|
||||
#if 0
|
||||
/* For debugging, check for gross overflow. */
|
||||
if( (e * y) > (MEXP + 1024) )
|
||||
goto overflow;
|
||||
#endif
|
||||
|
||||
/* Find significand of x in antilog table A[]. */
|
||||
i = 1;
|
||||
if( x <= douba(9) )
|
||||
i = 9;
|
||||
if( x <= douba(i+4) )
|
||||
i += 4;
|
||||
if( x <= douba(i+2) )
|
||||
i += 2;
|
||||
if( x >= douba(1) )
|
||||
i = -1;
|
||||
i += 1;
|
||||
|
||||
|
||||
/* Find (x - A[i])/A[i]
|
||||
* in order to compute log(x/A[i]):
|
||||
*
|
||||
* log(x) = log( a x/a ) = log(a) + log(x/a)
|
||||
*
|
||||
* log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a
|
||||
*/
|
||||
x -= douba(i);
|
||||
x -= doubb(i/2);
|
||||
x /= douba(i);
|
||||
|
||||
|
||||
/* rational approximation for log(1+v):
|
||||
*
|
||||
* log(1+v) = v - v**2/2 + v**3 P(v) / Q(v)
|
||||
*/
|
||||
z = x*x;
|
||||
w = x * ( z * polevl( x, P, 3 ) / p1evl( x, Q, 4 ) );
|
||||
w = w - ldexp( z, -1 ); /* w - 0.5 * z */
|
||||
|
||||
/* Convert to base 2 logarithm:
|
||||
* multiply by log2(e)
|
||||
*/
|
||||
w = w + LOG2EA * w;
|
||||
/* Note x was not yet added in
|
||||
* to above rational approximation,
|
||||
* so do it now, while multiplying
|
||||
* by log2(e).
|
||||
*/
|
||||
z = w + LOG2EA * x;
|
||||
z = z + x;
|
||||
|
||||
/* Compute exponent term of the base 2 logarithm. */
|
||||
w = -i;
|
||||
w = ldexp( w, -4 ); /* divide by 16 */
|
||||
w += e;
|
||||
/* Now base 2 log of x is w + z. */
|
||||
|
||||
/* Multiply base 2 log by y, in extended precision. */
|
||||
|
||||
/* separate y into large part ya
|
||||
* and small part yb less than 1/16
|
||||
*/
|
||||
ya = reduc(y);
|
||||
yb = y - ya;
|
||||
|
||||
|
||||
F = z * y + w * yb;
|
||||
Fa = reduc(F);
|
||||
Fb = F - Fa;
|
||||
|
||||
G = Fa + w * ya;
|
||||
Ga = reduc(G);
|
||||
Gb = G - Ga;
|
||||
|
||||
H = Fb + Gb;
|
||||
Ha = reduc(H);
|
||||
w = ldexp( Ga+Ha, 4 );
|
||||
|
||||
/* Test the power of 2 for overflow */
|
||||
if( w > MEXP )
|
||||
{
|
||||
#ifndef INFINITIES
|
||||
mtherr( fname, OVERFLOW );
|
||||
#endif
|
||||
#ifdef INFINITIES
|
||||
if( nflg && yoddint )
|
||||
return( -INFINITY );
|
||||
return( INFINITY );
|
||||
#else
|
||||
if( nflg && yoddint )
|
||||
return( -MAXNUM );
|
||||
return( MAXNUM );
|
||||
#endif
|
||||
}
|
||||
|
||||
if( w < (MNEXP - 1) )
|
||||
{
|
||||
#ifndef DENORMAL
|
||||
mtherr( fname, UNDERFLOW );
|
||||
#endif
|
||||
#ifdef MINUSZERO
|
||||
if( nflg && yoddint )
|
||||
return( NEGZERO );
|
||||
#endif
|
||||
return( 0.0 );
|
||||
}
|
||||
|
||||
e = w;
|
||||
Hb = H - Ha;
|
||||
|
||||
if( Hb > 0.0 )
|
||||
{
|
||||
e += 1;
|
||||
Hb -= 0.0625;
|
||||
}
|
||||
|
||||
/* Now the product y * log2(x) = Hb + e/16.0.
|
||||
*
|
||||
* Compute base 2 exponential of Hb,
|
||||
* where -0.0625 <= Hb <= 0.
|
||||
*/
|
||||
z = Hb * polevl( Hb, R, 6 ); /* z = 2**Hb - 1 */
|
||||
|
||||
/* Express e/16 as an integer plus a negative number of 16ths.
|
||||
* Find lookup table entry for the fractional power of 2.
|
||||
*/
|
||||
if( e < 0 )
|
||||
i = 0;
|
||||
else
|
||||
i = 1;
|
||||
i = e/16 + i;
|
||||
e = 16*i - e;
|
||||
w = douba( e );
|
||||
z = w + w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
|
||||
z = ldexp( z, i ); /* multiply by integer power of 2 */
|
||||
|
||||
done:
|
||||
|
||||
/* Negate if odd integer power of negative number */
|
||||
if( nflg && yoddint )
|
||||
{
|
||||
#ifdef MINUSZERO
|
||||
if( z == 0.0 )
|
||||
z = NEGZERO;
|
||||
else
|
||||
#endif
|
||||
z = -z;
|
||||
}
|
||||
return( z );
|
||||
}
|
||||
|
||||
|
||||
/* Find a multiple of 1/16 that is within 1/16 of x. */
|
||||
static __inline__ double reduc(x)
|
||||
double x;
|
||||
{
|
||||
double t;
|
||||
|
||||
t = ldexp( x, 4 );
|
||||
t = floor( t );
|
||||
t = ldexp( t, -4 );
|
||||
return(t);
|
||||
}
|
||||
200
winsup/mingw/mingwex/math/powi.c
Normal file
200
winsup/mingw/mingwex/math/powi.c
Normal file
@ -0,0 +1,200 @@
|
||||
/* 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);
|
||||
}
|
||||
198
winsup/mingw/mingwex/math/powif.c
Normal file
198
winsup/mingw/mingwex/math/powif.c
Normal file
@ -0,0 +1,198 @@
|
||||
/* powi.c
|
||||
*
|
||||
* Real raised to integer power
|
||||
*
|
||||
*
|
||||
*
|
||||
* SYNOPSIS:
|
||||
*
|
||||
* float x, y, __powif();
|
||||
* 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 float from powi.c and adapted to mingw
|
||||
2002-10-01 Danny Smith <dannysmith@users.sourceforge.net>
|
||||
*/
|
||||
|
||||
#ifdef __MINGW32__
|
||||
#include "cephes_mconf.h"
|
||||
#else
|
||||
#include "mconf.h"
|
||||
#ifdef ANSIPROT
|
||||
extern float logf ( float );
|
||||
extern float frexpf ( float, int * );
|
||||
extern int signbitf ( float );
|
||||
#else
|
||||
float logf(), frexpf();
|
||||
int signbitf();
|
||||
#endif
|
||||
extern float NEGZEROF, INFINITYF, MAXNUMF, MAXLOGF, MINLOGF, LOGE2F;
|
||||
#endif /* __MINGW32__ */
|
||||
|
||||
#ifndef _SET_ERRNO
|
||||
#define _SET_ERRNO(x)
|
||||
#endif
|
||||
|
||||
float __powif( float x, int nn )
|
||||
{
|
||||
int n, e, sign, asign, lx;
|
||||
float w, y, s;
|
||||
|
||||
/* See pow.c for these tests. */
|
||||
if( x == 0.0F )
|
||||
{
|
||||
if( nn == 0 )
|
||||
return( 1.0F );
|
||||
else if( nn < 0 )
|
||||
return( INFINITYF );
|
||||
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 = frexpf( 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 * LOGE2F;
|
||||
}
|
||||
else
|
||||
{
|
||||
s = LOGE2F * e;
|
||||
}
|
||||
|
||||
if( s > MAXLOGF )
|
||||
{
|
||||
mtherr( "__powif", OVERFLOW );
|
||||
_SET_ERRNO(ERANGE);
|
||||
y = INFINITYF;
|
||||
goto done;
|
||||
}
|
||||
|
||||
#if DENORMAL
|
||||
if( s < MINLOGF )
|
||||
{
|
||||
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 < (-MAXLOGF+2.0)) && (sign < 0) )
|
||||
{
|
||||
x = 1.0/x;
|
||||
sign = -sign;
|
||||
}
|
||||
#else
|
||||
/* do not produce denormal answer */
|
||||
if( s < -MAXLOGF )
|
||||
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 = NEGZEROF;
|
||||
else
|
||||
y = -y;
|
||||
}
|
||||
return(y);
|
||||
}
|
||||
@ -1,4 +1,4 @@
|
||||
/* powil.c
|
||||
/* __powil.c
|
||||
*
|
||||
* Real raised to integer power, long double precision
|
||||
*
|
||||
@ -6,10 +6,10 @@
|
||||
*
|
||||
* SYNOPSIS:
|
||||
*
|
||||
* long double x, y, powil();
|
||||
* long double x, y, __powil();
|
||||
* int n;
|
||||
*
|
||||
* y = powil( x, n );
|
||||
* y = __powil( x, n );
|
||||
*
|
||||
*
|
||||
*
|
||||
@ -36,7 +36,7 @@
|
||||
*
|
||||
*/
|
||||
|
||||
/* powil.c */
|
||||
/* __powil.c */
|
||||
|
||||
/*
|
||||
Cephes Math Library Release 2.2: December, 1990
|
||||
@ -66,7 +66,7 @@ long double frexpl();
|
||||
#define _SET_ERRNO(x)
|
||||
#endif
|
||||
|
||||
long double powil( x, nn )
|
||||
long double __powil( x, nn )
|
||||
long double x;
|
||||
int nn;
|
||||
{
|
||||
@ -126,7 +126,7 @@ else
|
||||
|
||||
if( s > MAXLOGL )
|
||||
{
|
||||
mtherr( "powil", OVERFLOW );
|
||||
mtherr( "__powil", OVERFLOW );
|
||||
_SET_ERRNO(ERANGE);
|
||||
y = INFINITYL;
|
||||
goto done;
|
||||
@ -134,7 +134,7 @@ if( s > MAXLOGL )
|
||||
|
||||
if( s < MINLOGL )
|
||||
{
|
||||
mtherr( "powil", UNDERFLOW );
|
||||
mtherr( "__powil", UNDERFLOW );
|
||||
_SET_ERRNO(ERANGE);
|
||||
return(0.0L);
|
||||
}
|
||||
|
||||
@ -382,7 +382,7 @@ static long double w, W, Wa, Wb, ya, yb, u;
|
||||
|
||||
#ifdef __MINGW32__
|
||||
static __inline__ long double reducl( long double );
|
||||
extern long double powil ( long double, int );
|
||||
extern long double __powil ( long double, int );
|
||||
extern long double powl ( long double x, long double y);
|
||||
#else
|
||||
#ifdef ANSIPROT
|
||||
@ -392,14 +392,14 @@ extern long double frexpl ( long double, int * );
|
||||
extern long double ldexpl ( long double, int );
|
||||
extern long double polevll ( long double, void *, int );
|
||||
extern long double p1evll ( long double, void *, int );
|
||||
extern long double powil ( long double, int );
|
||||
extern long double __powil ( long double, int );
|
||||
extern int isnanl ( long double );
|
||||
extern int isfinitel ( long double );
|
||||
static long double reducl( long double );
|
||||
extern int signbitl ( long double );
|
||||
#else
|
||||
long double floorl(), fabsl(), frexpl(), ldexpl();
|
||||
long double polevll(), p1evll(), powil();
|
||||
long double polevll(), p1evll(), __powil();
|
||||
static long double reducl();
|
||||
int isnanl(), isfinitel(), signbitl();
|
||||
#endif
|
||||
@ -603,7 +603,7 @@ if( iyflg )
|
||||
w = floorl(x);
|
||||
if( (w == x) && (fabsl(y) < 32768.0) )
|
||||
{
|
||||
w = powil( x, (int) y );
|
||||
w = __powil( x, (int) y );
|
||||
return( w );
|
||||
}
|
||||
}
|
||||
@ -764,40 +764,3 @@ t = ldexpl( t, -LNXT );
|
||||
return(t);
|
||||
}
|
||||
|
||||
/* powil.c
|
||||
*
|
||||
* Real raised to integer power, long double precision
|
||||
*
|
||||
*
|
||||
*
|
||||
* SYNOPSIS:
|
||||
*
|
||||
* long double x, y, powil();
|
||||
* int n;
|
||||
*
|
||||
* y = powil( 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
|
||||
* IEEE .001,1000 -1022,1023 50000 4.3e-17 7.8e-18
|
||||
* IEEE 1,2 -1022,1023 20000 3.9e-17 7.6e-18
|
||||
* IEEE .99,1.01 0,8700 10000 3.6e-16 7.2e-17
|
||||
*
|
||||
* Returns INFINITY on overflow, zero on underflow.
|
||||
*
|
||||
*/
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user