OrdinaryMagician/vkdoom_m
OrdinaryMagician/vkdoom_m
1
Fork 0
Code Activity

- moved most utility code to a subdirectory.

Browse source 
Again no code changes.
This commit is contained in:
Christoph Oelckers 2019-01-31 20:33:52 +01:00
commit 7241072b16
102 changed files with 46 additions and 40 deletions
Download patch file Download diff file Expand all files Collapse all files

315
src/utility/math/asin.c Normal file
View file

@ -0,0 +1,315 @@
/* asin.c
*
* Inverse circular sine
*
*
*
* SYNOPSIS:
*
* double x, y, asin();
*
* y = asin( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between -pi/2 and +pi/2 whose sine is x.
*
* A rational function of the form x + x**3 P(x**2)/Q(x**2)
* is used for |x| in the interval [0, 0.5]. If |x| > 0.5 it is
* transformed by the identity
*
* asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -1, 1 40000 2.6e-17 7.1e-18
* IEEE -1, 1 10^6 1.9e-16 5.4e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* asin domain |x| > 1 NAN
*
*/
/* acos()
*
* Inverse circular cosine
*
*
*
* SYNOPSIS:
*
* double x, y, acos();
*
* y = acos( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between 0 and pi whose cosine
* is x.
*
* Analytically, acos(x) = pi/2 - asin(x). However if |x| is
* near 1, there is cancellation error in subtracting asin(x)
* from pi/2. Hence if x < -0.5,
*
* acos(x) = pi - 2.0 * asin( sqrt((1+x)/2) );
*
* or if x > +0.5,
*
* acos(x) = 2.0 * asin( sqrt((1-x)/2) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -1, 1 50000 3.3e-17 8.2e-18
* IEEE -1, 1 10^6 2.2e-16 6.5e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* asin domain |x| > 1 NAN
*/
/* asin.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
/* arcsin(x) = x + x^3 P(x^2)/Q(x^2)
0 <= x <= 0.625
Peak relative error = 1.2e-18 */
#if UNK
static double P[6] = {
4.253011369004428248960E-3,
-6.019598008014123785661E-1,
5.444622390564711410273E0,
-1.626247967210700244449E1,
1.956261983317594739197E1,
-8.198089802484824371615E0,
};
static double Q[5] = {
/* 1.000000000000000000000E0, */
-1.474091372988853791896E1,
7.049610280856842141659E1,
-1.471791292232726029859E2,
1.395105614657485689735E2,
-4.918853881490881290097E1,
};
#endif
#if DEC
static short P[24] = {
0036213,0056330,0057244,0053234,
0140032,0015011,0114762,0160255,
0040656,0035130,0136121,0067313,
0141202,0014616,0170474,0101731,
0041234,0100076,0151674,0111310,
0141003,0025540,0033165,0077246,
};
static short Q[20] = {
/* 0040200,0000000,0000000,0000000, */
0141153,0155310,0055360,0072530,
0041614,0177001,0027764,0101237,
0142023,0026733,0064653,0133266,
0042013,0101264,0023775,0176351,
0141504,0140420,0050660,0036543,
};
#endif
#if IBMPC
static short P[24] = {
0x8ad3,0x0bd4,0x6b9b,0x3f71,
0x5c16,0x333e,0x4341,0xbfe3,
0x2dd9,0x178a,0xc74b,0x4015,
0x907b,0xde27,0x4331,0xc030,
0x9259,0xda77,0x9007,0x4033,
0xafd5,0x06ce,0x656c,0xc020,
};
static short Q[20] = {
/* 0x0000,0x0000,0x0000,0x3ff0, */
0x0eab,0x0b5e,0x7b59,0xc02d,
0x9054,0x25fe,0x9fc0,0x4051,
0x76d7,0x6d35,0x65bb,0xc062,
0xbf9d,0x84ff,0x7056,0x4061,
0x07ac,0x0a36,0x9822,0xc048,
};
#endif
#if MIEEE
static short P[24] = {
0x3f71,0x6b9b,0x0bd4,0x8ad3,
0xbfe3,0x4341,0x333e,0x5c16,
0x4015,0xc74b,0x178a,0x2dd9,
0xc030,0x4331,0xde27,0x907b,
0x4033,0x9007,0xda77,0x9259,
0xc020,0x656c,0x06ce,0xafd5,
};
static short Q[20] = {
/* 0x3ff0,0x0000,0x0000,0x0000, */
0xc02d,0x7b59,0x0b5e,0x0eab,
0x4051,0x9fc0,0x25fe,0x9054,
0xc062,0x65bb,0x6d35,0x76d7,
0x4061,0x7056,0x84ff,0xbf9d,
0xc048,0x9822,0x0a36,0x07ac,
};
#endif
/* arcsin(1-x) = pi/2 - sqrt(2x)(1+R(x))
0 <= x <= 0.5
Peak relative error = 4.2e-18 */
#if UNK
static double R[5] = {
2.967721961301243206100E-3,
-5.634242780008963776856E-1,
6.968710824104713396794E0,
-2.556901049652824852289E1,
2.853665548261061424989E1,
};
static double S[4] = {
/* 1.000000000000000000000E0, */
-2.194779531642920639778E1,
1.470656354026814941758E2,
-3.838770957603691357202E2,
3.424398657913078477438E2,
};
#endif
#if DEC
static short R[20] = {
0036102,0077034,0142164,0174103,
0140020,0036222,0147711,0044173,
0040736,0177655,0153631,0171523,
0141314,0106525,0060015,0055474,
0041344,0045422,0003630,0040344,
};
static short S[16] = {
/* 0040200,0000000,0000000,0000000, */
0141257,0112425,0132772,0166136,
0042023,0010315,0075523,0175020,
0142277,0170104,0126203,0017563,
0042253,0034115,0102662,0022757,
};
#endif
#if IBMPC
static short R[20] = {
0x9f08,0x988e,0x4fc3,0x3f68,
0x290f,0x59f9,0x0792,0xbfe2,
0x3e6a,0xbaf3,0xdff5,0x401b,
0xab68,0xac01,0x91aa,0xc039,
0x081d,0x40f3,0x8962,0x403c,
};
static short S[16] = {
/* 0x0000,0x0000,0x0000,0x3ff0, */
0x5d8c,0xb6bf,0xf2a2,0xc035,
0x7f42,0xaf6a,0x6219,0x4062,
0x63ee,0x9590,0xfe08,0xc077,
0x44be,0xb0b6,0x6709,0x4075,
};
#endif
#if MIEEE
static short R[20] = {
0x3f68,0x4fc3,0x988e,0x9f08,
0xbfe2,0x0792,0x59f9,0x290f,
0x401b,0xdff5,0xbaf3,0x3e6a,
0xc039,0x91aa,0xac01,0xab68,
0x403c,0x8962,0x40f3,0x081d,
};
static short S[16] = {
/* 0x3ff0,0x0000,0x0000,0x0000, */
0xc035,0xf2a2,0xb6bf,0x5d8c,
0x4062,0x6219,0xaf6a,0x7f42,
0xc077,0xfe08,0x9590,0x63ee,
0x4075,0x6709,0xb0b6,0x44be,
};
#endif
/* pi/2 = PIO2 + MOREBITS. */
#ifdef DEC
#define MOREBITS 5.721188726109831840122E-18
#else
#define MOREBITS 6.123233995736765886130E-17
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double c_sqrt ( double );
double c_asin ( double );
#else
double c_sqrt(), polevl(), p1evl();
double c_asin();
#endif
extern double PIO2, PIO4, NAN;
double c_asin(x)
double x;
{
double a, p, z, zz;
short sign;
if( x > 0 )
{
sign = 1;
a = x;
}
else
{
sign = -1;
a = -x;
}
if( a > 1.0 )
{
mtherr( "asin", DOMAIN );
return( NAN );
}
if( a > 0.625 )
{
/* arcsin(1-x) = pi/2 - sqrt(2x)(1+R(x)) */
zz = 1.0 - a;
p = zz * polevl( zz, R, 4)/p1evl( zz, S, 4);
zz = c_sqrt(zz+zz);
z = PIO4 - zz;
zz = zz * p - MOREBITS;
z = z - zz;
z = z + PIO4;
}
else
{
if( a < 1.0e-8 )
{
return(x);
}
zz = a * a;
z = zz * polevl( zz, P, 5)/p1evl( zz, Q, 5);
z = a * z + a;
}
if( sign < 0 )
z = -z;
return(z);
}
double c_acos(x)
double x;
{
if( (x < -1.0) || (x > 1.0) )
{
mtherr( "acos", DOMAIN );
return( NAN );
}
return PIO2 - c_asin(x) + MOREBITS;
}

393
src/utility/math/atan.c Normal file
View file

@ -0,0 +1,393 @@
/* atan.c
*
* Inverse circular tangent
* (arctangent)
*
*
*
* SYNOPSIS:
*
* double x, y, atan();
*
* y = atan( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between -pi/2 and +pi/2 whose tangent
* is x.
*
* Range reduction is from three intervals into the interval
* from zero to 0.66. The approximant uses a rational
* function of degree 4/5 of the form x + x**3 P(x)/Q(x).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -10, 10 50000 2.4e-17 8.3e-18
* IEEE -10, 10 10^6 1.8e-16 5.0e-17
*
*/
/* atan2()
*
* Quadrant correct inverse circular tangent
*
*
*
* SYNOPSIS:
*
* double x, y, z, atan2();
*
* z = atan2( y, x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle whose tangent is y/x.
* Define compile time symbol ANSIC = 1 for ANSI standard,
* range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range
* 0 to 2PI, args (x,y).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10, 10 10^6 2.5e-16 6.9e-17
* See atan.c.
*
*/
/* atan.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
/* arctan(x) = x + x^3 P(x^2)/Q(x^2)
0 <= x <= 0.66
Peak relative error = 2.6e-18 */
#ifdef UNK
static double P[5] = {
-8.750608600031904122785E-1,
-1.615753718733365076637E1,
-7.500855792314704667340E1,
-1.228866684490136173410E2,
-6.485021904942025371773E1,
};
static double Q[5] = {
/* 1.000000000000000000000E0, */
2.485846490142306297962E1,
1.650270098316988542046E2,
4.328810604912902668951E2,
4.853903996359136964868E2,
1.945506571482613964425E2,
};
/* tan( 3*pi/8 ) */
static double T3P8 = 2.41421356237309504880;
#endif
#ifdef DEC
static short P[20] = {
0140140,0001775,0007671,0026242,
0141201,0041242,0155534,0001715,
0141626,0002141,0132100,0011625,
0141765,0142771,0064055,0150453,
0141601,0131517,0164507,0062164,
};
static short Q[20] = {
/* 0040200,0000000,0000000,0000000, */
0041306,0157042,0154243,0000742,
0042045,0003352,0016707,0150452,
0042330,0070306,0113425,0170730,
0042362,0130770,0116602,0047520,
0042102,0106367,0156753,0013541,
};
/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {040432,0101171,0114774,0167462,};
#define T3P8 *(double *)T3P8A
#endif
#ifdef IBMPC
static short P[20] = {
0x2594,0xa1f7,0x007f,0xbfec,
0x807a,0x5b6b,0x2854,0xc030,
0x0273,0x3688,0xc08c,0xc052,
0xba25,0x2d05,0xb8bf,0xc05e,
0xec8e,0xfd28,0x3669,0xc050,
};
static short Q[20] = {
/* 0x0000,0x0000,0x0000,0x3ff0, */
0x603c,0x5b14,0xdbc4,0x4038,
0xfa25,0x43b8,0xa0dd,0x4064,
0xbe3b,0xd2e2,0x0e18,0x407b,
0x49ea,0x13b0,0x563f,0x407e,
0x62ec,0xfbbd,0x519e,0x4068,
};
/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {0x9de6,0x333f,0x504f,0x4003};
#define T3P8 *(double *)T3P8A
#endif
#ifdef MIEEE
static short P[20] = {
0xbfec,0x007f,0xa1f7,0x2594,
0xc030,0x2854,0x5b6b,0x807a,
0xc052,0xc08c,0x3688,0x0273,
0xc05e,0xb8bf,0x2d05,0xba25,
0xc050,0x3669,0xfd28,0xec8e,
};
static short Q[20] = {
/* 0x3ff0,0x0000,0x0000,0x0000, */
0x4038,0xdbc4,0x5b14,0x603c,
0x4064,0xa0dd,0x43b8,0xfa25,
0x407b,0x0e18,0xd2e2,0xbe3b,
0x407e,0x563f,0x13b0,0x49ea,
0x4068,0x519e,0xfbbd,0x62ec,
};
/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {
0x4003,0x504f,0x333f,0x9de6
};
#define T3P8 *(double *)T3P8A
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double atan ( double );
extern double fabs ( double );
extern int signbit ( double );
extern int isnan ( double );
#else
double polevl(), p1evl(), atan(), fabs();
int signbit(), isnan();
#endif
extern double PI, PIO2, PIO4, INFINITY, NEGZERO, MAXNUM;
/* pi/2 = PIO2 + MOREBITS. */
#ifdef DEC
#define MOREBITS 5.721188726109831840122E-18
#else
#define MOREBITS 6.123233995736765886130E-17
#endif
double c_atan(x)
double x;
{
double y, z;
short sign, flag;
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
#ifdef INFINITIES
if(x == INFINITY)
return(PIO2);
if(x == -INFINITY)
return(-PIO2);
#endif
/* make argument positive and save the sign */
sign = 1;
if( x < 0.0 )
{
sign = -1;
x = -x;
}
/* range reduction */
flag = 0;
if( x > T3P8 )
{
y = PIO2;
flag = 1;
x = -( 1.0/x );
}
else if( x <= 0.66 )
{
y = 0.0;
}
else
{
y = PIO4;
flag = 2;
x = (x-1.0)/(x+1.0);
}
z = x * x;
z = z * polevl( z, P, 4 ) / p1evl( z, Q, 5 );
z = x * z + x;
if( flag == 2 )
z += 0.5 * MOREBITS;
else if( flag == 1 )
z += MOREBITS;
y = y + z;
if( sign < 0 )
y = -y;
return(y);
}
/* atan2 */
#ifdef ANSIC
double c_atan2( y, x )
#else
double c_atan2( x, y )
#endif
double x, y;
{
double z, w;
short code;
code = 0;
#ifdef NANS
if( isnan(x) )
return(x);
if( isnan(y) )
return(y);
#endif
#ifdef MINUSZERO
if( y == 0.0 )
{
if( signbit(y) )
{
if( x > 0.0 )
z = y;
else if( x < 0.0 )
z = -PI;
else
{
if( signbit(x) )
z = -PI;
else
z = y;
}
}
else /* y is +0 */
{
if( x == 0.0 )
{
if( signbit(x) )
z = PI;
else
z = 0.0;
}
else if( x > 0.0 )
z = 0.0;
else
z = PI;
}
return z;
}
if( x == 0.0 )
{
if( y > 0.0 )
z = PIO2;
else
z = -PIO2;
return z;
}
#endif /* MINUSZERO */
#ifdef INFINITIES
if( x == INFINITY )
{
if( y == INFINITY )
z = 0.25 * PI;
else if( y == -INFINITY )
z = -0.25 * PI;
else if( y < 0.0 )
z = NEGZERO;
else
z = 0.0;
return z;
}
if( x == -INFINITY )
{
if( y == INFINITY )
z = 0.75 * PI;
else if( y <= -INFINITY )
z = -0.75 * PI;
else if( y >= 0.0 )
z = PI;
else
z = -PI;
return z;
}
if( y == INFINITY )
return( PIO2 );
if( y == -INFINITY )
return( -PIO2 );
#endif
if( x < 0.0 )
code = 2;
if( y < 0.0 )
code |= 1;
#ifdef INFINITIES
if( x == 0.0 )
#else
if( fabs(x) <= (fabs(y) / MAXNUM) )
#endif
{
if( code & 1 )
{
#if ANSIC
return( -PIO2 );
#else
return( 3.0*PIO2 );
#endif
}
if( y == 0.0 )
return( 0.0 );
return( PIO2 );
}
if( y == 0.0 )
{
if( code & 2 )
return( PI );
return( 0.0 );
}
switch( code )
{
#if ANSIC
default:
case 0:
case 1: w = 0.0; break;
case 2: w = PI; break;
case 3: w = -PI; break;
#else
default:
case 0: w = 0.0; break;
case 1: w = 2.0 * PI; break;
case 2:
case 3: w = PI; break;
#endif
}
z = w + c_atan( y/x );
#ifdef MINUSZERO
if( z == 0.0 && y < 0 )
z = NEGZERO;
#endif
return( z );
}

149
src/utility/math/cmath.h Normal file
View file

@ -0,0 +1,149 @@
#ifndef __CMATH_H
#define __CMATH_H
#include "xs_Float.h"
#define USE_CUSTOM_MATH // we want repreducably reliable results, even at the cost of performance
#define USE_FAST_MATH // use faster table-based sin and cos variants with limited precision (sufficient for Doom gameplay)
extern"C"
{
double c_asin(double);
double c_acos(double);
double c_atan(double);
double c_atan2(double, double);
double c_sin(double);
double c_cos(double);
double c_tan(double);
double c_cot(double);
double c_sqrt(double);
double c_sinh(double);
double c_cosh(double);
double c_tanh(double);
double c_exp(double);
double c_log(double);
double c_log10(double);
double c_pow(double, double);
}
// This uses a sine table with linear interpolation
// For in-game calculations this is precise enough
// and this code is more than 10x faster than the
// Cephes sin and cos function.
struct FFastTrig
{
static const int TBLPERIOD = 8192;
static const int BITSHIFT = 19;
static const int REMAINDER = (1 << BITSHIFT) - 1;
float sinetable[2049];
double sinq1(unsigned);
public:
FFastTrig();
double sin(unsigned);
double cos(unsigned);
};
extern FFastTrig fasttrig;
// This must use xs_Float to guarantee proper integer wraparound.
#define DEG2BAM(f) ((unsigned)xs_CRoundToInt((f) * (0x40000000/90.)))
#define RAD2BAM(f) ((unsigned)xs_CRoundToInt((f) * (0x80000000/3.14159265358979323846)))
inline double fastcosdeg(double v)
{
return fasttrig.cos(DEG2BAM(v));
}
inline double fastsindeg(double v)
{
return fasttrig.sin(DEG2BAM(v));
}
inline double fastcos(double v)
{
return fasttrig.cos(RAD2BAM(v));
}
inline double fastsin(double v)
{
return fasttrig.sin(RAD2BAM(v));
}
// these are supposed to be local to this file.
#undef DEG2BAM
#undef RAD2BAM
inline double sindeg(double v)
{
#ifdef USE_CUSTOM_MATH
return c_sin(v * (3.14159265358979323846 / 180.));
#else
return sin(v * (3.14159265358979323846 / 180.));
#endif
}
inline double cosdeg(double v)
{
#ifdef USE_CUSTOM_MATH
return c_cos(v * (3.14159265358979323846 / 180.));
#else
return cos(v * (3.14159265358979323846 / 180.));
#endif
}
#ifndef USE_CUSTOM_MATH
#define g_asin asin
#define g_acos acos
#define g_atan atan
#define g_atan2 atan2
#define g_sin sin
#define g_cos cos
#define g_sindeg sindeg
#define g_cosdeg cosdeg
#define g_tan tan
#define g_cot cot
#define g_sqrt sqrt
#define g_sinh sinh
#define g_cosh cosh
#define g_tanh tanh
#define g_exp exp
#define g_log log
#define g_log10 log10
#define g_pow pow
#else
#define g_asin c_asin
#define g_acos c_acos
#define g_atan c_atan
#define g_atan2 c_atan2
#ifndef USE_FAST_MATH
#define g_sindeg sindeg
#define g_cosdeg cosdeg
#define g_sin c_sin
#define g_cos c_cos
#else
#define g_sindeg fastsindeg
#define g_cosdeg fastcosdeg
#define g_sin fastsin
#define g_cos fastcos
#endif
#define g_tan c_tan
#define g_cot c_cot
#define g_sqrt c_sqrt
#define g_sinh c_sinh
#define g_cosh c_cosh
#define g_tanh c_tanh
#define g_exp c_exp
#define g_log c_log
#define g_log10 c_log10
#define g_pow c_pow
#endif
#endif

252
src/utility/math/const.c Normal file
View file

@ -0,0 +1,252 @@
/* const.c
*
* Globally declared constants
*
*
*
* SYNOPSIS:
*
* extern double nameofconstant;
*
*
*
*
* DESCRIPTION:
*
* This file contains a number of mathematical constants and
* also some needed size parameters of the computer arithmetic.
* The values are supplied as arrays of hexadecimal integers
* for IEEE arithmetic; arrays of octal constants for DEC
* arithmetic; and in a normal decimal scientific notation for
* other machines. The particular notation used is determined
* by a symbol (DEC, IBMPC, or UNK) defined in the include file
* mconf.h.
*
* The default size parameters are as follows.
*
* For DEC and UNK modes:
* MACHEP = 1.38777878078144567553E-17 2**-56
* MAXLOG = 8.8029691931113054295988E1 log(2**127)
* MINLOG = -8.872283911167299960540E1 log(2**-128)
* MAXNUM = 1.701411834604692317316873e38 2**127
*
* For IEEE arithmetic (IBMPC):
* MACHEP = 1.11022302462515654042E-16 2**-53
* MAXLOG = 7.09782712893383996843E2 log(2**1024)
* MINLOG = -7.08396418532264106224E2 log(2**-1022)
* MAXNUM = 1.7976931348623158E308 2**1024
*
* The global symbols for mathematical constants are
* PI = 3.14159265358979323846 pi
* PIO2 = 1.57079632679489661923 pi/2
* PIO4 = 7.85398163397448309616E-1 pi/4
* SQRT2 = 1.41421356237309504880 sqrt(2)
* SQRTH = 7.07106781186547524401E-1 sqrt(2)/2
* LOG2E = 1.4426950408889634073599 1/log(2)
* SQ2OPI = 7.9788456080286535587989E-1 sqrt( 2/pi )
* LOGE2 = 6.93147180559945309417E-1 log(2)
* LOGSQ2 = 3.46573590279972654709E-1 log(2)/2
* THPIO4 = 2.35619449019234492885 3*pi/4
* TWOOPI = 6.36619772367581343075535E-1 2/pi
*
* These lists are subject to change.
*/
/* const.c */
/*
Cephes Math Library Release 2.3: March, 1995
Copyright 1984, 1995 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
#if 1
double MACHEP = 1.11022302462515654042E-16; /* 2**-53 */
#else
double MACHEP = 1.38777878078144567553E-17; /* 2**-56 */
#endif
double UFLOWTHRESH = 2.22507385850720138309E-308; /* 2**-1022 */
#ifdef DENORMAL
double MAXLOG = 7.09782712893383996732E2; /* log(MAXNUM) */
/* double MINLOG = -7.44440071921381262314E2; */ /* log(2**-1074) */
double MINLOG = -7.451332191019412076235E2; /* log(2**-1075) */
#else
double MAXLOG = 7.08396418532264106224E2; /* log 2**1022 */
double MINLOG = -7.08396418532264106224E2; /* log 2**-1022 */
#endif
double MAXNUM = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
double PI = 3.14159265358979323846; /* pi */
double PIO2 = 1.57079632679489661923; /* pi/2 */
double PIO4 = 7.85398163397448309616E-1; /* pi/4 */
double SQRT2 = 1.41421356237309504880; /* sqrt(2) */
double SQRTH = 7.07106781186547524401E-1; /* sqrt(2)/2 */
double LOG2E = 1.4426950408889634073599; /* 1/log(2) */
double SQ2OPI = 7.9788456080286535587989E-1; /* sqrt( 2/pi ) */
double LOGE2 = 6.93147180559945309417E-1; /* log(2) */
double LOGSQ2 = 3.46573590279972654709E-1; /* log(2)/2 */
double THPIO4 = 2.35619449019234492885; /* 3*pi/4 */
double TWOOPI = 6.36619772367581343075535E-1; /* 2/pi */
#ifdef INFINITIES
double INFINITY = 1.0/0.0; /* 99e999; */
#else
double INFINITY = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
#endif
#ifdef NANS
double NAN = 1.0/0.0 - 1.0/0.0;
#else
double NAN = 0.0;
#endif
#ifdef MINUSZERO
double NEGZERO = -0.0;
#else
double NEGZERO = 0.0;
#endif
#endif
#ifdef IBMPC
/* 2**-53 = 1.11022302462515654042E-16 */
unsigned short MACHEP[4] = {0x0000,0x0000,0x0000,0x3ca0};
unsigned short UFLOWTHRESH[4] = {0x0000,0x0000,0x0000,0x0010};
#ifdef DENORMAL
/* log(MAXNUM) = 7.09782712893383996732224E2 */
unsigned short MAXLOG[4] = {0x39ef,0xfefa,0x2e42,0x4086};
/* log(2**-1074) = - -7.44440071921381262314E2 */
/*unsigned short MINLOG[4] = {0x71c3,0x446d,0x4385,0xc087};*/
unsigned short MINLOG[4] = {0x3052,0xd52d,0x4910,0xc087};
#else
/* log(2**1022) = 7.08396418532264106224E2 */
unsigned short MAXLOG[4] = {0xbcd2,0xdd7a,0x232b,0x4086};
/* log(2**-1022) = - 7.08396418532264106224E2 */
unsigned short MINLOG[4] = {0xbcd2,0xdd7a,0x232b,0xc086};
#endif
/* 2**1024*(1-MACHEP) = 1.7976931348623158E308 */
unsigned short MAXNUM[4] = {0xffff,0xffff,0xffff,0x7fef};
unsigned short PI[4] = {0x2d18,0x5444,0x21fb,0x4009};
unsigned short PIO2[4] = {0x2d18,0x5444,0x21fb,0x3ff9};
unsigned short PIO4[4] = {0x2d18,0x5444,0x21fb,0x3fe9};
unsigned short SQRT2[4] = {0x3bcd,0x667f,0xa09e,0x3ff6};
unsigned short SQRTH[4] = {0x3bcd,0x667f,0xa09e,0x3fe6};
unsigned short LOG2E[4] = {0x82fe,0x652b,0x1547,0x3ff7};
unsigned short SQ2OPI[4] = {0x3651,0x33d4,0x8845,0x3fe9};
unsigned short LOGE2[4] = {0x39ef,0xfefa,0x2e42,0x3fe6};
unsigned short LOGSQ2[4] = {0x39ef,0xfefa,0x2e42,0x3fd6};
unsigned short THPIO4[4] = {0x21d2,0x7f33,0xd97c,0x4002};
unsigned short TWOOPI[4] = {0xc883,0x6dc9,0x5f30,0x3fe4};
#ifdef INFINITIES
unsigned short INFINITY[4] = {0x0000,0x0000,0x0000,0x7ff0};
#else
unsigned short INFINITY[4] = {0xffff,0xffff,0xffff,0x7fef};
#endif
#ifdef NANS
unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x7ffc};
#else
unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#ifdef MINUSZERO
unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x8000};
#else
unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#endif
#ifdef MIEEE
/* 2**-53 = 1.11022302462515654042E-16 */
unsigned short MACHEP[4] = {0x3ca0,0x0000,0x0000,0x0000};
unsigned short UFLOWTHRESH[4] = {0x0010,0x0000,0x0000,0x0000};
#ifdef DENORMAL
/* log(2**1024) = 7.09782712893383996843E2 */
unsigned short MAXLOG[4] = {0x4086,0x2e42,0xfefa,0x39ef};
/* log(2**-1074) = - -7.44440071921381262314E2 */
/* unsigned short MINLOG[4] = {0xc087,0x4385,0x446d,0x71c3}; */
unsigned short MINLOG[4] = {0xc087,0x4910,0xd52d,0x3052};
#else
/* log(2**1022) = 7.08396418532264106224E2 */
unsigned short MAXLOG[4] = {0x4086,0x232b,0xdd7a,0xbcd2};
/* log(2**-1022) = - 7.08396418532264106224E2 */
unsigned short MINLOG[4] = {0xc086,0x232b,0xdd7a,0xbcd2};
#endif
/* 2**1024*(1-MACHEP) = 1.7976931348623158E308 */
unsigned short MAXNUM[4] = {0x7fef,0xffff,0xffff,0xffff};
unsigned short PI[4] = {0x4009,0x21fb,0x5444,0x2d18};
unsigned short PIO2[4] = {0x3ff9,0x21fb,0x5444,0x2d18};
unsigned short PIO4[4] = {0x3fe9,0x21fb,0x5444,0x2d18};
unsigned short SQRT2[4] = {0x3ff6,0xa09e,0x667f,0x3bcd};
unsigned short SQRTH[4] = {0x3fe6,0xa09e,0x667f,0x3bcd};
unsigned short LOG2E[4] = {0x3ff7,0x1547,0x652b,0x82fe};
unsigned short SQ2OPI[4] = {0x3fe9,0x8845,0x33d4,0x3651};
unsigned short LOGE2[4] = {0x3fe6,0x2e42,0xfefa,0x39ef};
unsigned short LOGSQ2[4] = {0x3fd6,0x2e42,0xfefa,0x39ef};
unsigned short THPIO4[4] = {0x4002,0xd97c,0x7f33,0x21d2};
unsigned short TWOOPI[4] = {0x3fe4,0x5f30,0x6dc9,0xc883};
#ifdef INFINITIES
unsigned short INFINITY[4] = {0x7ff0,0x0000,0x0000,0x0000};
#else
unsigned short INFINITY[4] = {0x7fef,0xffff,0xffff,0xffff};
#endif
#ifdef NANS
unsigned short NAN[4] = {0x7ff8,0x0000,0x0000,0x0000};
#else
unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#ifdef MINUSZERO
unsigned short NEGZERO[4] = {0x8000,0x0000,0x0000,0x0000};
#else
unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#endif
#ifdef DEC
/* 2**-56 = 1.38777878078144567553E-17 */
unsigned short MACHEP[4] = {0022200,0000000,0000000,0000000};
unsigned short UFLOWTHRESH[4] = {0x0080,0x0000,0x0000,0x0000};
/* log 2**127 = 88.029691931113054295988 */
unsigned short MAXLOG[4] = {041660,007463,0143742,025733,};
/* log 2**-128 = -88.72283911167299960540 */
unsigned short MINLOG[4] = {0141661,071027,0173721,0147572,};
/* 2**127 = 1.701411834604692317316873e38 */
unsigned short MAXNUM[4] = {077777,0177777,0177777,0177777,};
unsigned short PI[4] = {040511,007732,0121041,064302,};
unsigned short PIO2[4] = {040311,007732,0121041,064302,};
unsigned short PIO4[4] = {040111,007732,0121041,064302,};
unsigned short SQRT2[4] = {040265,002363,031771,0157145,};
unsigned short SQRTH[4] = {040065,002363,031771,0157144,};
unsigned short LOG2E[4] = {040270,0125073,024534,013761,};
unsigned short SQ2OPI[4] = {040114,041051,0117241,0131204,};
unsigned short LOGE2[4] = {040061,071027,0173721,0147572,};
unsigned short LOGSQ2[4] = {037661,071027,0173721,0147572,};
unsigned short THPIO4[4] = {040426,0145743,0174631,007222,};
unsigned short TWOOPI[4] = {040042,0174603,067116,042025,};
/* Approximate infinity by MAXNUM. */
unsigned short INFINITY[4] = {077777,0177777,0177777,0177777,};
unsigned short NAN[4] = {0000000,0000000,0000000,0000000};
#ifdef MINUSZERO
unsigned short NEGZERO[4] = {0000000,0000000,0000000,0100000};
#else
unsigned short NEGZERO[4] = {0000000,0000000,0000000,0000000};
#endif
#endif
#ifndef UNK
extern unsigned short MACHEP[];
extern unsigned short UFLOWTHRESH[];
extern unsigned short MAXLOG[];
extern unsigned short UNDLOG[];
extern unsigned short MINLOG[];
extern unsigned short MAXNUM[];
extern unsigned short PI[];
extern unsigned short PIO2[];
extern unsigned short PIO4[];
extern unsigned short SQRT2[];
extern unsigned short SQRTH[];
extern unsigned short LOG2E[];
extern unsigned short SQ2OPI[];
extern unsigned short LOGE2[];
extern unsigned short LOGSQ2[];
extern unsigned short THPIO4[];
extern unsigned short TWOOPI[];
extern unsigned short INFINITY[];
extern unsigned short NAN[];
extern unsigned short NEGZERO[];
#endif

83
src/utility/math/cosh.c Normal file
View file

@ -0,0 +1,83 @@
/* cosh.c
*
* Hyperbolic cosine
*
*
*
* SYNOPSIS:
*
* double x, y, cosh();
*
* y = cosh( x );
*
*
*
* DESCRIPTION:
*
* Returns hyperbolic cosine of argument in the range MINLOG to
* MAXLOG.
*
* cosh(x) = ( exp(x) + exp(-x) )/2.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC +- 88 50000 4.0e-17 7.7e-18
* IEEE +-MAXLOG 30000 2.6e-16 5.7e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* cosh overflow |x| > MAXLOG MAXNUM
*
*
*/
/* cosh.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1985, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double c_exp ( double );
extern int isnan ( double );
extern int isfinite ( double );
#else
double c_exp();
int isnan(), isfinite();
#endif
extern double MAXLOG, INFINITY, LOGE2;
double c_cosh(x)
double x;
{
double y;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
if( x < 0 )
x = -x;
if( x > (MAXLOG + LOGE2) )
{
mtherr( "cosh", OVERFLOW );
return( INFINITY );
}
if( x >= (MAXLOG - LOGE2) )
{
y = c_exp(0.5 * x);
y = (0.5 * y) * y;
return(y);
}
y = c_exp(x);
y = 0.5 * (y + 1.0 / y);
return( y );
}

182
src/utility/math/exp.c Normal file
View file

@ -0,0 +1,182 @@
/* exp.c
*
* Exponential function
*
*
*
* SYNOPSIS:
*
* double x, y, exp();
*
* y = exp( x );
*
*
*
* DESCRIPTION:
*
* Returns e (2.71828...) raised to the x power.
*
* Range reduction is accomplished by separating the argument
* into an integer k and fraction f such that
*
* x k f
* e = 2 e.
*
* A Pade' form 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
* of degree 2/3 is used to approximate exp(f) in the basic
* interval [-0.5, 0.5].
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0, MAXLOG 38000 3.0e-17 6.2e-18
* IEEE +- 708 40000 2.0e-16 5.6e-17
*
*
* Error amplification in the exponential function can be
* a serious matter. The error propagation involves
* exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
* which shows that a 1 lsb error in representing X produces
* a relative error of X times 1 lsb in the function.
* While the routine gives an accurate result for arguments
* that are exactly represented by a double precision
* computer number, the result contains amplified roundoff
* error for large arguments not exactly represented.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* exp underflow x < MINLOG 0.0
* exp overflow x > MAXLOG MAXNUM
*
*/
/*
Cephes Math Library Release 2.2: January, 1991
Copyright 1984, 1991 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* Exponential function */
#include "mconf.h"
static char fname[] = {"exp"};
#ifdef UNK
static double P[] = {
1.26177193074810590878E-4,
3.02994407707441961300E-2,
9.99999999999999999910E-1,
};
static double Q[] = {
3.00198505138664455042E-6,
2.52448340349684104192E-3,
2.27265548208155028766E-1,
2.00000000000000000009E0,
};
static double C1 = 6.93145751953125E-1;
static double C2 = 1.42860682030941723212E-6;
#endif
#ifdef DEC
static short P[] = {
0035004,0047156,0127442,0057502,
0036770,0033210,0063121,0061764,
0040200,0000000,0000000,0000000,
};
static short Q[] = {
0033511,0072665,0160662,0176377,
0036045,0070715,0124105,0132777,
0037550,0134114,0142077,0001637,
0040400,0000000,0000000,0000000,
};
static short sc1[] = {0040061,0071000,0000000,0000000};
#define C1 (*(double *)sc1)
static short sc2[] = {0033277,0137216,0075715,0057117};
#define C2 (*(double *)sc2)
#endif
#ifdef IBMPC
static short P[] = {
0x4be8,0xd5e4,0x89cd,0x3f20,
0x2c7e,0x0cca,0x06d1,0x3f9f,
0x0000,0x0000,0x0000,0x3ff0,
};
static short Q[] = {
0x5fa0,0xbc36,0x2eb6,0x3ec9,
0xb6c0,0xb508,0xae39,0x3f64,
0xe074,0x9887,0x1709,0x3fcd,
0x0000,0x0000,0x0000,0x4000,
};
static short sc1[] = {0x0000,0x0000,0x2e40,0x3fe6};
#define C1 (*(double *)sc1)
static short sc2[] = {0xabca,0xcf79,0xf7d1,0x3eb7};
#define C2 (*(double *)sc2)
#endif
#ifdef MIEEE
static short P[] = {
0x3f20,0x89cd,0xd5e4,0x4be8,
0x3f9f,0x06d1,0x0cca,0x2c7e,
0x3ff0,0x0000,0x0000,0x0000,
};
static short Q[] = {
0x3ec9,0x2eb6,0xbc36,0x5fa0,
0x3f64,0xae39,0xb508,0xb6c0,
0x3fcd,0x1709,0x9887,0xe074,
0x4000,0x0000,0x0000,0x0000,
};
static short sc1[] = {0x3fe6,0x2e40,0x0000,0x0000};
#define C1 (*(double *)sc1)
static short sc2[] = {0x3eb7,0xf7d1,0xcf79,0xabca};
#define C2 (*(double *)sc2)
#endif
extern double LOGE2, LOG2E, MAXLOG, MINLOG, MAXNUM;
double c_exp(x)
double x;
{
double px, xx;
int n;
double polevl(), floor(), ldexp();
if( x > MAXLOG)
{
mtherr( fname, OVERFLOW );
return( MAXNUM );
}
if( x < MINLOG )
{
mtherr( fname, UNDERFLOW );
return(0.0);
}
/* Express e**x = e**g 2**n
* = e**g e**( n loge(2) )
* = e**( g + n loge(2) )
*/
px = floor( LOG2E * x + 0.5 ); /* floor() truncates toward -infinity. */
n = (int)px;
x -= px * C1;
x -= px * C2;
/* rational approximation for exponential
* of the fractional part:
* e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
*/
xx = x * x;
px = x * polevl( xx, P, 2 );
x = px/( polevl( xx, Q, 3 ) - px );
x = 1.0 + ldexp( x, 1 );
/* multiply by power of 2 */
x = ldexp( x, n );
return(x);
}

106
src/utility/math/fastsin.cpp Normal file
View file

@ -0,0 +1,106 @@
/*
** fastsin.cpp
** a table/linear interpolation-based sine function that is both
** precise and fast enough for most purposes.
**
**---------------------------------------------------------------------------
** Copyright 2015 Christoph Oelckers
** All rights reserved.
**
** Redistribution and use in source and binary forms, with or without
** modification, are permitted provided that the following conditions
** are met:
**
** 1. Redistributions of source code must retain the above copyright
** notice, this list of conditions and the following disclaimer.
** 2. Redistributions in binary form must reproduce the above copyright
** notice, this list of conditions and the following disclaimer in the
** documentation and/or other materials provided with the distribution.
** 3. The name of the author may not be used to endorse or promote products
** derived from this software without specific prior written permission.
**
** THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
** IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
** OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
** IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
** INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
** NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
** DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
** THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
** (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
** THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
**---------------------------------------------------------------------------
**
*/
#include <math.h>
#include "cmath.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
FFastTrig fasttrig;
FFastTrig::FFastTrig()
{
const double pimul = M_PI * 2 / TBLPERIOD;
for (int i = 0; i < 2049; i++)
{
sinetable[i] = (float)c_sin(i*pimul);
}
}
__forceinline double FFastTrig::sinq1(unsigned bangle)
{
unsigned int index = bangle >> BITSHIFT;
if ((bangle &= (REMAINDER)) == 0) // This is to avoid precision problems at 180°
{
return double(sinetable[index]);
}
else
{
return (double(sinetable[index]) * (REMAINDER - bangle) + double(sinetable[index + 1]) * bangle) * (1. / REMAINDER);
}
}
double FFastTrig::sin(unsigned bangle)
{
switch (bangle & 0xc0000000)
{
default:
return sinq1(bangle);
case 0x40000000:
return sinq1(0x80000000 - bangle);
case 0x80000000:
return -sinq1(bangle - 0x80000000);
case 0xc0000000:
return -sinq1(0 - bangle);
}
}
double FFastTrig::cos(unsigned bangle)
{
switch (bangle & 0xc0000000)
{
default:
return sinq1(0x40000000 - bangle);
case 0x40000000:
return -sinq1(bangle - 0x40000000);
case 0x80000000:
return -sinq1(0xc0000000 - bangle);
case 0xc0000000:
return sinq1(bangle - 0xc0000000);
}
}

237
src/utility/math/isnan.c Normal file
View file

@ -0,0 +1,237 @@
/* isnan()
* signbit()
* isfinite()
*
* Floating point numeric utilities
*
*
*
* SYNOPSIS:
*
* double ceil(), floor(), frexp(), ldexp();
* int signbit(), isnan(), isfinite();
* double x, y;
* int expnt, n;
*
* y = floor(x);
* y = ceil(x);
* y = frexp( x, &expnt );
* y = ldexp( x, n );
* n = signbit(x);
* n = isnan(x);
* n = isfinite(x);
*
*
*
* DESCRIPTION:
*
* All four routines return a double precision floating point
* result.
*
* floor() returns the largest integer less than or equal to x.
* It truncates toward minus infinity.
*
* ceil() returns the smallest integer greater than or equal
* to x. It truncates toward plus infinity.
*
* frexp() extracts the exponent from x. It returns an integer
* power of two to expnt and the significand between 0.5 and 1
* to y. Thus x = y * 2**expn.
*
* ldexp() multiplies x by 2**n.
*
* signbit(x) returns 1 if the sign bit of x is 1, else 0.
*
* These functions are part of the standard C run time library
* for many but not all C compilers. The ones supplied are
* written in C for either DEC or IEEE arithmetic. They should
* be used only if your compiler library does not already have
* them.
*
* The IEEE versions assume that denormal numbers are implemented
* in the arithmetic. Some modifications will be required if
* the arithmetic has abrupt rather than gradual underflow.
*/
/*
Cephes Math Library Release 2.3: March, 1995
Copyright 1984, 1995 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
/* ceil(), floor(), frexp(), ldexp() may need to be rewritten. */
#undef UNK
#if BIGENDIAN
#define MIEEE 1
#else
#define IBMPC 1
#endif
#endif
/* Return 1 if the sign bit of x is 1, else 0. */
int signbit(x)
double x;
{
union
{
double d;
short s[4];
int i[2];
} u;
u.d = x;
if( sizeof(int) == 4 )
{
#ifdef IBMPC
return( u.i[1] < 0 );
#endif
#ifdef DEC
return( u.s[3] < 0 );
#endif
#ifdef MIEEE
return( u.i[0] < 0 );
#endif
}
else
{
#ifdef IBMPC
return( u.s[3] < 0 );
#endif
#ifdef DEC
return( u.s[3] < 0 );
#endif
#ifdef MIEEE
return( u.s[0] < 0 );
#endif
}
}
/* Return 1 if x is a number that is Not a Number, else return 0. */
int isnan(x)
double x;
{
#ifdef NANS
union
{
double d;
unsigned short s[4];
unsigned int i[2];
} u;
u.d = x;
if( sizeof(int) == 4 )
{
#ifdef IBMPC
if( ((u.i[1] & 0x7ff00000) == 0x7ff00000)
&& (((u.i[1] & 0x000fffff) != 0) || (u.i[0] != 0)))
return 1;
#endif
#ifdef DEC
if( (u.s[1] & 0x7fff) == 0)
{
if( (u.s[2] | u.s[1] | u.s[0]) != 0 )
return(1);
}
#endif
#ifdef MIEEE
if( ((u.i[0] & 0x7ff00000) == 0x7ff00000)
&& (((u.i[0] & 0x000fffff) != 0) || (u.i[1] != 0)))
return 1;
#endif
return(0);
}
else
{ /* size int not 4 */
#ifdef IBMPC
if( (u.s[3] & 0x7ff0) == 0x7ff0)
{
if( ((u.s[3] & 0x000f) | u.s[2] | u.s[1] | u.s[0]) != 0 )
return(1);
}
#endif
#ifdef DEC
if( (u.s[3] & 0x7fff) == 0)
{
if( (u.s[2] | u.s[1] | u.s[0]) != 0 )
return(1);
}
#endif
#ifdef MIEEE
if( (u.s[0] & 0x7ff0) == 0x7ff0)
{
if( ((u.s[0] & 0x000f) | u.s[1] | u.s[2] | u.s[3]) != 0 )
return(1);
}
#endif
return(0);
} /* size int not 4 */
#else
/* No NANS. */
return(0);
#endif
}
/* Return 1 if x is not infinite and is not a NaN. */
int isfinite(x)
double x;
{
#ifdef INFINITIES
union
{
double d;
unsigned short s[4];
unsigned int i[2];
} u;
u.d = x;
if( sizeof(int) == 4 )
{
#ifdef IBMPC
if( (u.i[1] & 0x7ff00000) != 0x7ff00000)
return 1;
#endif
#ifdef DEC
if( (u.s[3] & 0x7fff) != 0)
return 1;
#endif
#ifdef MIEEE
if( (u.i[0] & 0x7ff00000) != 0x7ff00000)
return 1;
#endif
return(0);
}
else
{
#ifdef IBMPC
if( (u.s[3] & 0x7ff0) != 0x7ff0)
return 1;
#endif
#ifdef DEC
if( (u.s[3] & 0x7fff) != 0)
return 1;
#endif
#ifdef MIEEE
if( (u.s[0] & 0x7ff0) != 0x7ff0)
return 1;
#endif
return(0);
}
#else
/* No INFINITY. */
return(1);
#endif
}

341
src/utility/math/log.c Normal file
View file

@ -0,0 +1,341 @@
/* log.c
*
* Natural logarithm
*
*
*
* SYNOPSIS:
*
* double x, y, log();
*
* y = log( x );
*
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z**3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 150000 1.44e-16 5.06e-17
* IEEE +-MAXNUM 30000 1.20e-16 4.78e-17
* DEC 0, 10 170000 1.8e-17 6.3e-18
*
* In the tests over the interval [+-MAXNUM], the logarithms
* of the random arguments were uniformly distributed over
* [0, MAXLOG].
*
* ERROR MESSAGES:
*
* log singularity: x = 0; returns -INFINITY
* log domain: x < 0; returns NAN
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
static char fname[] = {"log"};
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
*/
#ifdef UNK
static double P[] = {
1.01875663804580931796E-4,
4.97494994976747001425E-1,
4.70579119878881725854E0,
1.44989225341610930846E1,
1.79368678507819816313E1,
7.70838733755885391666E0,
};
static double Q[] = {
/* 1.00000000000000000000E0, */
1.12873587189167450590E1,
4.52279145837532221105E1,
8.29875266912776603211E1,
7.11544750618563894466E1,
2.31251620126765340583E1,
};
#endif
#ifdef DEC
static unsigned short P[] = {
0037777,0127270,0162547,0057274,
0041001,0054665,0164317,0005341,
0041451,0034104,0031640,0105773,
0041677,0011276,0123617,0160135,
0041701,0126603,0053215,0117250,
0041420,0115777,0135206,0030232,
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041220,0144332,0045272,0174241,
0041742,0164566,0035720,0130431,
0042246,0126327,0166065,0116357,
0042372,0033420,0157525,0124560,
0042271,0167002,0066537,0172303,
0041730,0164777,0113711,0044407,
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x1bb0,0x93c3,0xb4c2,0x3f1a,
0x52f2,0x3f56,0xd6f5,0x3fdf,
0x6911,0xed92,0xd2ba,0x4012,
0xeb2e,0xc63e,0xff72,0x402c,
0xc84d,0x924b,0xefd6,0x4031,
0xdcf8,0x7d7e,0xd563,0x401e,
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0xef8e,0xae97,0x9320,0x4026,
0xc033,0x4e19,0x9d2c,0x4046,
0xbdbd,0xa326,0xbf33,0x4054,
0xae21,0xeb5e,0xc9e2,0x4051,
0x25b2,0x9e1f,0x200a,0x4037,
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x3f1a,0xb4c2,0x93c3,0x1bb0,
0x3fdf,0xd6f5,0x3f56,0x52f2,
0x4012,0xd2ba,0xed92,0x6911,
0x402c,0xff72,0xc63e,0xeb2e,
0x4031,0xefd6,0x924b,0xc84d,
0x401e,0xd563,0x7d7e,0xdcf8,
};
static unsigned short Q[] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0x4026,0x9320,0xae97,0xef8e,
0x4046,0x9d2c,0x4e19,0xc033,
0x4054,0xbf33,0xa326,0xbdbd,
0x4051,0xc9e2,0xeb5e,0xae21,
0x4037,0x200a,0x9e1f,0x25b2,
};
#endif
/* Coefficients for log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
*/
#ifdef UNK
static double R[3] = {
-7.89580278884799154124E-1,
1.63866645699558079767E1,
-6.41409952958715622951E1,
};
static double S[3] = {
/* 1.00000000000000000000E0,*/
-3.56722798256324312549E1,
3.12093766372244180303E2,
-7.69691943550460008604E2,
};
#endif
#ifdef DEC
static unsigned short R[12] = {
0140112,0020756,0161540,0072035,
0041203,0013743,0114023,0155527,
0141600,0044060,0104421,0050400,
};
static unsigned short S[12] = {
/*0040200,0000000,0000000,0000000,*/
0141416,0130152,0017543,0064122,
0042234,0006000,0104527,0020155,
0142500,0066110,0146631,0174731,
};
#endif
#ifdef IBMPC
static unsigned short R[12] = {
0x0e84,0xdc6c,0x443d,0xbfe9,
0x7b6b,0x7302,0x62fc,0x4030,
0x2a20,0x1122,0x0906,0xc050,
};
static unsigned short S[12] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x6d0a,0x43ec,0xd60d,0xc041,
0xe40e,0x112a,0x8180,0x4073,
0x3f3b,0x19b3,0x0d89,0xc088,
};
#endif
#ifdef MIEEE
static unsigned short R[12] = {
0xbfe9,0x443d,0xdc6c,0x0e84,
0x4030,0x62fc,0x7302,0x7b6b,
0xc050,0x0906,0x1122,0x2a20,
};
static unsigned short S[12] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0xc041,0xd60d,0x43ec,0x6d0a,
0x4073,0x8180,0x112a,0xe40e,
0xc088,0x0d89,0x19b3,0x3f3b,
};
#endif
#ifdef ANSIPROT
extern double frexp ( double, int * );
extern double ldexp ( double, int );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double frexp(), ldexp(), polevl(), p1evl();
int isnan(), isfinite();
#endif
#define SQRTH 0.70710678118654752440
extern double INFINITY, NAN;
double c_log(x)
double x;
{
int e;
#ifdef DEC
short *q;
#endif
double y, z;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
#ifdef INFINITIES
if( x == INFINITY )
return(x);
#endif
/* Test for domain */
if( x <= 0.0 )
{
if( x == 0.0 )
{
mtherr( fname, SING );
return( -INFINITY );
}
else
{
mtherr( fname, DOMAIN );
return( NAN );
}
}
/* separate mantissa from exponent */
#ifdef DEC
q = (short *)&x;
e = *q; /* short containing exponent */
e = ((e >> 7) & 0377) - 0200; /* the exponent */
*q &= 0177; /* strip exponent from x */
*q |= 040000; /* x now between 0.5 and 1 */
#endif
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
#ifdef IBMPC
x = frexp( x, &e );
/*
q = (short *)&x;
q += 3;
e = *q;
e = ((e >> 4) & 0x0fff) - 0x3fe;
*q &= 0x0f;
*q |= 0x3fe0;
*/
#endif
/* Equivalent C language standard library function: */
#ifdef UNK
x = frexp( x, &e );
#endif
#ifdef MIEEE
x = frexp( x, &e );
#endif
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if( (e > 2) || (e < -2) )
{
if( x < SQRTH )
{ /* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5;
y = 0.5 * z + 0.5;
}
else
{ /* 2 (x-1)/(x+1) */
z = x - 0.5;
z -= 0.5;
y = 0.5 * x + 0.5;
}
x = z / y;
/* rational form */
z = x*x;
z = x * ( z * polevl( z, R, 2 ) / p1evl( z, S, 3 ) );
y = e;
z = z - y * 2.121944400546905827679e-4;
z = z + x;
z = z + e * 0.693359375;
goto ldone;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if( x < SQRTH )
{
e -= 1;
x = ldexp( x, 1 ) - 1.0; /* 2x - 1 */
}
else
{
x = x - 1.0;
}
/* rational form */
z = x*x;
#if DEC
y = x * ( z * polevl( x, P, 5 ) / p1evl( x, Q, 6 ) );
#else
y = x * ( z * polevl( x, P, 5 ) / p1evl( x, Q, 5 ) );
#endif
if( e )
y = y - e * 2.121944400546905827679e-4;
y = y - ldexp( z, -1 ); /* y - 0.5 * z */
z = x + y;
if( e )
z = z + e * 0.693359375;
ldone:
return( z );
}

250
src/utility/math/log10.c Normal file
View file

@ -0,0 +1,250 @@
/* log10.c
*
* Common logarithm
*
*
*
* SYNOPSIS:
*
* double x, y, log10();
*
* y = log10( x );
*
*
*
* DESCRIPTION:
*
* Returns logarithm to the base 10 of x.
*
* The argument is separated into its exponent and fractional
* parts. The logarithm of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 30000 1.5e-16 5.0e-17
* IEEE 0, MAXNUM 30000 1.4e-16 4.8e-17
* DEC 1, MAXNUM 50000 2.5e-17 6.0e-18
*
* In the tests over the interval [1, MAXNUM], the logarithms
* of the random arguments were uniformly distributed over
* [0, MAXLOG].
*
* ERROR MESSAGES:
*
* log10 singularity: x = 0; returns -INFINITY
* log10 domain: x < 0; returns NAN
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
static char fname[] = {"log10"};
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
*/
#ifdef UNK
static double P[] = {
4.58482948458143443514E-5,
4.98531067254050724270E-1,
6.56312093769992875930E0,
2.97877425097986925891E1,
6.06127134467767258030E1,
5.67349287391754285487E1,
1.98892446572874072159E1
};
static double Q[] = {
/* 1.00000000000000000000E0, */
1.50314182634250003249E1,
8.27410449222435217021E1,
2.20664384982121929218E2,
3.07254189979530058263E2,
2.14955586696422947765E2,
5.96677339718622216300E1
};
#endif
#ifdef DEC
static unsigned short P[] = {
0034500,0046473,0051374,0135174,
0037777,0037566,0145712,0150321,
0040722,0002426,0031543,0123107,
0041356,0046513,0170752,0004346,
0041562,0071553,0023536,0163343,
0041542,0170221,0024316,0114216,
0041237,0016454,0046611,0104602
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041160,0100260,0067736,0102424,
0041645,0075552,0036563,0147072,
0042134,0125025,0021132,0025320,
0042231,0120211,0046030,0103271,
0042126,0172241,0052151,0120426,
0041556,0125702,0072116,0047103
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x974f,0x6a5f,0x09a7,0x3f08,
0x5a1a,0xd979,0xe7ee,0x3fdf,
0x74c9,0xc66c,0x40a2,0x401a,
0x411d,0x7e3d,0xc9a9,0x403d,
0xdcdc,0x64eb,0x4e6d,0x404e,
0xd312,0x2519,0x5e12,0x404c,
0x3130,0x89b1,0xe3a5,0x4033
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0xd0a2,0x0dfb,0x1016,0x402e,
0x79c7,0x47ae,0xaf6d,0x4054,
0x455a,0xa44b,0x9542,0x406b,
0x10d7,0x2983,0x3411,0x4073,
0x3423,0x2a8d,0xde94,0x406a,
0xc9c8,0x4e89,0xd578,0x404d
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x3f08,0x09a7,0x6a5f,0x974f,
0x3fdf,0xe7ee,0xd979,0x5a1a,
0x401a,0x40a2,0xc66c,0x74c9,
0x403d,0xc9a9,0x7e3d,0x411d,
0x404e,0x4e6d,0x64eb,0xdcdc,
0x404c,0x5e12,0x2519,0xd312,
0x4033,0xe3a5,0x89b1,0x3130
};
static unsigned short Q[] = {
0x402e,0x1016,0x0dfb,0xd0a2,
0x4054,0xaf6d,0x47ae,0x79c7,
0x406b,0x9542,0xa44b,0x455a,
0x4073,0x3411,0x2983,0x10d7,
0x406a,0xde94,0x2a8d,0x3423,
0x404d,0xd578,0x4e89,0xc9c8
};
#endif
#define SQRTH 0.70710678118654752440
#define L102A 3.0078125E-1
#define L102B 2.48745663981195213739E-4
#define L10EA 4.3359375E-1
#define L10EB 7.00731903251827651129E-4
#ifdef ANSIPROT
extern double frexp ( double, int * );
extern double ldexp ( double, int );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double frexp(), ldexp(), polevl(), p1evl();
int isnan(), isfinite();
#endif
extern double LOGE2, SQRT2, INFINITY, NAN;
double c_log10(x)
double x;
{
VOLATILE double z;
double y;
#ifdef DEC
short *q;
#endif
int e;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
#ifdef INFINITIES
if( x == INFINITY )
return(x);
#endif
/* Test for domain */
if( x <= 0.0 )
{
if( x == 0.0 )
{
mtherr( fname, SING );
return( -INFINITY );
}
else
{
mtherr( fname, DOMAIN );
return( NAN );
}
}
/* separate mantissa from exponent */
#ifdef DEC
q = (short *)&x;
e = *q; /* short containing exponent */
e = ((e >> 7) & 0377) - 0200; /* the exponent */
*q &= 0177; /* strip exponent from x */
*q |= 040000; /* x now between 0.5 and 1 */
#endif
#ifdef IBMPC
x = frexp( x, &e );
/*
q = (short *)&x;
q += 3;
e = *q;
e = ((e >> 4) & 0x0fff) - 0x3fe;
*q &= 0x0f;
*q |= 0x3fe0;
*/
#endif
/* Equivalent C language standard library function: */
#ifdef UNK
x = frexp( x, &e );
#endif
#ifdef MIEEE
x = frexp( x, &e );
#endif
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if( x < SQRTH )
{
e -= 1;
x = ldexp( x, 1 ) - 1.0; /* 2x - 1 */
}
else
{
x = x - 1.0;
}
/* rational form */
z = x*x;
y = x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) );
y = y - ldexp( z, -1 ); /* y - 0.5 * x**2 */
/* multiply log of fraction by log10(e)
* and base 2 exponent by log10(2)
*/
z = (x + y) * L10EB; /* accumulate terms in order of size */
z += y * L10EA;
z += x * L10EA;
z += e * L102B;
z += e * L102A;
return( z );
}

199
src/utility/math/mconf.h Normal file
View file

@ -0,0 +1,199 @@
/* mconf.h
*
* Common include file for math routines
*
*
*
* SYNOPSIS:
*
* #include "mconf.h"
*
*
*
* DESCRIPTION:
*
* This file contains definitions for error codes that are
* passed to the common error handling routine mtherr()
* (which see).
*
* The file also includes a conditional assembly definition
* for the type of computer arithmetic (IEEE, DEC, Motorola
* IEEE, or UNKnown).
*
* For Digital Equipment PDP-11 and VAX computers, certain
* IBM systems, and others that use numbers with a 56-bit
* significand, the symbol DEC should be defined. In this
* mode, most floating point constants are given as arrays
* of octal integers to eliminate decimal to binary conversion
* errors that might be introduced by the compiler.
*
* For little-endian computers, such as IBM PC, that follow the
* IEEE Standard for Binary Floating Point Arithmetic (ANSI/IEEE
* Std 754-1985), the symbol IBMPC should be defined. These
* numbers have 53-bit significands. In this mode, constants
* are provided as arrays of hexadecimal 16 bit integers.
*
* Big-endian IEEE format is denoted MIEEE. On some RISC
* systems such as Sun SPARC, double precision constants
* must be stored on 8-byte address boundaries. Since integer
* arrays may be aligned differently, the MIEEE configuration
* may fail on such machines.
*
* To accommodate other types of computer arithmetic, all
* constants are also provided in a normal decimal radix
* which one can hope are correctly converted to a suitable
* format by the available C language compiler. To invoke
* this mode, define the symbol UNK.
*
* An important difference among these modes is a predefined
* set of machine arithmetic constants for each. The numbers
* MACHEP (the machine roundoff error), MAXNUM (largest number
* represented), and several other parameters are preset by
* the configuration symbol. Check the file const.c to
* ensure that these values are correct for your computer.
*
* Configurations NANS, INFINITIES, MINUSZERO, and DENORMAL
* may fail on many systems. Verify that they are supposed
* to work on your computer.
*/
/*
Cephes Math Library Release 2.3: June, 1995
Copyright 1984, 1987, 1989, 1995 by Stephen L. Moshier
*/
/* Define if the `long double' type works. */
//#define HAVE_LONG_DOUBLE 0
/* Define as the return type of signal handlers (int or void). */
#define RETSIGTYPE void
/* Define if you have the ANSI C header files. */
#define STDC_HEADERS 1
/* Define if your processor stores words with the most significant
byte first (like Motorola and SPARC, unlike Intel and VAX). */
/* #undef WORDS_BIGENDIAN */
/* Define if floating point words are bigendian. */
/* #undef FLOAT_WORDS_BIGENDIAN */
/* The number of bytes in a int. */
#define SIZEOF_INT 4
/* Define if you have the <string.h> header file. */
#define HAVE_STRING_H 1
/* Name of package */
#define PACKAGE "cephes"
/* Version number of package */
#define VERSION "2.7"
/* Constant definitions for math error conditions
*/
#define DOMAIN 1 /* argument domain error */
#define SING 2 /* argument singularity */
#define OVERFLOW 3 /* overflow range error */
#define UNDERFLOW 4 /* underflow range error */
#define TLOSS 5 /* total loss of precision */
#define PLOSS 6 /* partial loss of precision */
#define EDOM 33
#define ERANGE 34
/* Complex numeral. */
typedef struct
{
double r;
double i;
} cmplx;
#ifdef HAVE_LONG_DOUBLE
/* Long double complex numeral. */
typedef struct
{
long double r;
long double i;
} cmplxl;
#endif
/* Type of computer arithmetic */
/* PDP-11, Pro350, VAX:
*/
/* #define DEC 1 */
/* Intel IEEE, low order words come first:
*/
//#define IBMPC 1
/* Motorola IEEE, high order words come first
* (Sun 680x0 workstation):
*/
/* #define MIEEE 1 */
/* UNKnown arithmetic, invokes coefficients given in
* normal decimal format. Beware of range boundary
* problems (MACHEP, MAXLOG, etc. in const.c) and
* roundoff problems in pow.c:
* (Sun SPARCstation)
*/
#define UNK 1
/* If you define UNK, then be sure to set BIGENDIAN properly. */
#ifdef FLOAT_WORDS_BIGENDIAN
#define BIGENDIAN 1
#else
#define BIGENDIAN 0
#endif
/* Define this `volatile' if your compiler thinks
* that floating point arithmetic obeys the associative
* and distributive laws. It will defeat some optimizations
* (but probably not enough of them).
*
* #define VOLATILE volatile
*/
#define VOLATILE
/* For 12-byte long doubles on an i386, pad a 16-bit short 0
* to the end of real constants initialized by integer arrays.
*
* #define XPD 0,
*
* Otherwise, the type is 10 bytes long and XPD should be
* defined blank (e.g., Microsoft C).
*
* #define XPD
*/
#define XPD 0,
/* Define to support tiny denormal numbers, else undefine. */
#define DENORMAL 1
/* Define to ask for infinity support, else undefine. */
#define INFINITIES 1
/* Define to ask for support of numbers that are Not-a-Number,
else undefine. This may automatically define INFINITIES in some files. */
#define NANS 1
/* Define to distinguish between -0.0 and +0.0. */
#define MINUSZERO 1
/* Define 1 for ANSI C atan2() function
See atan.c and clog.c. */
#define ANSIC 1
/* Get ANSI function prototypes, if you want them. */
#if 1
/* #ifdef __STDC__ */
#define ANSIPROT 1
int mtherr ( char *, int );
#else
int mtherr();
#endif
/* Variable for error reporting. See mtherr.c. */
extern int merror;

102
src/utility/math/mtherr.c Normal file
View file

@ -0,0 +1,102 @@
/* mtherr.c
*
* Library common error handling routine
*
*
*
* SYNOPSIS:
*
* char *fctnam;
* int code;
* int mtherr();
*
* mtherr( fctnam, code );
*
*
*
* DESCRIPTION:
*
* This routine may be called to report one of the following
* error conditions (in the include file mconf.h).
*
* Mnemonic Value Significance
*
* DOMAIN 1 argument domain error
* SING 2 function singularity
* OVERFLOW 3 overflow range error
* UNDERFLOW 4 underflow range error
* TLOSS 5 total loss of precision
* PLOSS 6 partial loss of precision
* EDOM 33 Unix domain error code
* ERANGE 34 Unix range error code
*
* The default version of the file prints the function name,
* passed to it by the pointer fctnam, followed by the
* error condition. The display is directed to the standard
* output device. The routine then returns to the calling
* program. Users may wish to modify the program to abort by
* calling exit() under severe error conditions such as domain
* errors.
*
* Since all error conditions pass control to this function,
* the display may be easily changed, eliminated, or directed
* to an error logging device.
*
* SEE ALSO:
*
* mconf.h
*
*/
/*
Cephes Math Library Release 2.0: April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include <stdio.h>
#include "mconf.h"
int merror = 0;
/* Notice: the order of appearance of the following
* messages is bound to the error codes defined
* in mconf.h.
*/
static char *ermsg[7] = {
"unknown", /* error code 0 */
"domain", /* error code 1 */
"singularity", /* et seq. */
"overflow",
"underflow",
"total loss of precision",
"partial loss of precision"
};
int mtherr( name, code )
char *name;
int code;
{
/* Display string passed by calling program,
* which is supposed to be the name of the
* function in which the error occurred:
*/
printf( "\n%s ", name );
/* Set global error message word */
merror = code;
/* Display error message defined
* by the code argument.
*/
if( (code <= 0) || (code >= 7) )
code = 0;
printf( "%s error\n", ermsg[code] );
/* Return to calling
* program
*/
return( 0 );
}

97
src/utility/math/polevl.c Normal file
View file

@ -0,0 +1,97 @@
/* 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.
*
*/
/*
Cephes Math Library Release 2.1: December, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
double polevl( x, coef, N )
double x;
double coef[];
int N;
{
double ans;
int i;
double *p;
p = coef;
ans = *p++;
i = N;
do
ans = ans * x + *p++;
while( --i );
return( ans );
}
/* p1evl() */
/* N
* Evaluate polynomial when coefficient of x is 1.0.
* Otherwise same as polevl.
*/
double p1evl( x, coef, N )
double x;
double coef[];
int N;
{
double ans;
double *p;
int i;
p = coef;
ans = x + *p++;
i = N-1;
do
ans = ans * x + *p++;
while( --i );
return( ans );
}

756
src/utility/math/pow.c Normal file
View file

@ -0,0 +1,756 @@
/* 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
*/
#include "mconf.h"
static char fname[] = {"pow"};
#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 unsigned short P[] = {
0x5cf0,0x7f5b,0xdb99,0x3fdf,
0xdf15,0xea9e,0xddef,0x400d,
0xeb6f,0x7f78,0xccbd,0x401e,
0x9b74,0xb65c,0xaa83,0x4012,
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x914e,0x9b20,0xaab4,0x4022,
0xc9f5,0x41c1,0xffff,0x403b,
0x6402,0x1b17,0xccbc,0x4040,
0xe92e,0x918a,0xffc5,0x402b,
};
static 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 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 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 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 c_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(), c_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
double c_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) )
return( x );
if( isnan(y) )
return( 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 )
{
#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 )
{
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 );
#ifdef NANS
return(NAN);
#else
return(0.0L);
#endif
}
nflg = 1;
}
}
/* Integer power of an integer. */
if( iyflg )
{
i = (int)w;
w = floor(x);
if( (w == x) && (fabs(y) < 32768.0) )
{
w = c_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 = (int)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 double reduc(x)
double x;
{
double t;
t = ldexp( x, 4 );
t = floor( t );
t = ldexp( t, -4 );
return(t);
}

186
src/utility/math/powi.c Normal file
View file

@ -0,0 +1,186 @@
/* 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
*/
#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;
double c_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 );
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);
}

15
src/utility/math/readme.txt Normal file
View file

@ -0,0 +1,15 @@
Some software in this archive may be from the book _Methods and
Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
International, 1989) or from the Cephes Mathematical Library, a
commercial product. In either event, it is copyrighted by the author.
What you see here may be used freely but it comes with no support or
guarantee.
The two known misprints in the book are repaired here in the
source listings for the gamma function and the incomplete beta
integral.
Stephen L. Moshier
moshier@na-net.ornl.gov

387
src/utility/math/sin.c Normal file
View file

@ -0,0 +1,387 @@
/* sin.c
*
* Circular sine
*
*
*
* SYNOPSIS:
*
* double x, y, sin();
*
* y = sin( x );
*
*
*
* DESCRIPTION:
*
* Range reduction is into intervals of pi/4. The reduction
* error is nearly eliminated by contriving an extended precision
* modular arithmetic.
*
* Two polynomial approximating functions are employed.
* Between 0 and pi/4 the sine is approximated by
* x + x**3 P(x**2).
* Between pi/4 and pi/2 the cosine is represented as
* 1 - x**2 Q(x**2).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0, 10 150000 3.0e-17 7.8e-18
* IEEE -1.07e9,+1.07e9 130000 2.1e-16 5.4e-17
*
* ERROR MESSAGES:
*
* message condition value returned
* sin total loss x > 1.073741824e9 0.0
*
* Partial loss of accuracy begins to occur at x = 2**30
* = 1.074e9. The loss is not gradual, but jumps suddenly to
* about 1 part in 10e7. Results may be meaningless for
* x > 2**49 = 5.6e14. The routine as implemented flags a
* TLOSS error for x > 2**30 and returns 0.0.
*/
/* cos.c
*
* Circular cosine
*
*
*
* SYNOPSIS:
*
* double x, y, cos();
*
* y = cos( x );
*
*
*
* DESCRIPTION:
*
* Range reduction is into intervals of pi/4. The reduction
* error is nearly eliminated by contriving an extended precision
* modular arithmetic.
*
* Two polynomial approximating functions are employed.
* Between 0 and pi/4 the cosine is approximated by
* 1 - x**2 Q(x**2).
* Between pi/4 and pi/2 the sine is represented as
* x + x**3 P(x**2).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1.07e9,+1.07e9 130000 2.1e-16 5.4e-17
* DEC 0,+1.07e9 17000 3.0e-17 7.2e-18
*/
/* sin.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1985, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
static double sincof[] = {
1.58962301576546568060E-10,
-2.50507477628578072866E-8,
2.75573136213857245213E-6,
-1.98412698295895385996E-4,
8.33333333332211858878E-3,
-1.66666666666666307295E-1,
};
static double coscof[6] = {
-1.13585365213876817300E-11,
2.08757008419747316778E-9,
-2.75573141792967388112E-7,
2.48015872888517045348E-5,
-1.38888888888730564116E-3,
4.16666666666665929218E-2,
};
static double DP1 = 7.85398125648498535156E-1;
static double DP2 = 3.77489470793079817668E-8;
static double DP3 = 2.69515142907905952645E-15;
/* static double lossth = 1.073741824e9; */
#endif
#ifdef DEC
static unsigned short sincof[] = {
0030056,0143750,0177214,0163153,
0131727,0027455,0044510,0175352,
0033470,0167432,0131752,0042414,
0135120,0006400,0146776,0174027,
0036410,0104210,0104207,0137202,
0137452,0125252,0125252,0125103,
};
static unsigned short coscof[24] = {
0127107,0151115,0002060,0152325,
0031017,0072353,0155161,0174053,
0132623,0171173,0172542,0057056,
0034320,0006400,0147102,0023652,
0135666,0005540,0133012,0076213,
0037052,0125252,0125252,0125126,
};
/* 7.853981629014015197753906250000E-1 */
static unsigned short P1[] = {0040111,0007732,0120000,0000000,};
/* 4.960467869796758577649598009884E-10 */
static unsigned short P2[] = {0030410,0055060,0100000,0000000,};
/* 2.860594363054915898381331279295E-18 */
static unsigned short P3[] = {0021523,0011431,0105056,0001560,};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
#endif
#ifdef IBMPC
static unsigned short sincof[] = {
0x9ccd,0x1fd1,0xd8fd,0x3de5,
0x1f5d,0xa929,0xe5e5,0xbe5a,
0x48a1,0x567d,0x1de3,0x3ec7,
0xdf03,0x19bf,0x01a0,0xbf2a,
0xf7d0,0x1110,0x1111,0x3f81,
0x5548,0x5555,0x5555,0xbfc5,
};
static unsigned short coscof[24] = {
0x1a9b,0xa086,0xfa49,0xbda8,
0x3f05,0x7b4e,0xee9d,0x3e21,
0x4bc6,0x7eac,0x7e4f,0xbe92,
0x44f5,0x19c8,0x01a0,0x3efa,
0x4f91,0x16c1,0xc16c,0xbf56,
0x554b,0x5555,0x5555,0x3fa5,
};
/*
7.85398125648498535156E-1,
3.77489470793079817668E-8,
2.69515142907905952645E-15,
*/
static unsigned short P1[] = {0x0000,0x4000,0x21fb,0x3fe9};
static unsigned short P2[] = {0x0000,0x0000,0x442d,0x3e64};
static unsigned short P3[] = {0x5170,0x98cc,0x4698,0x3ce8};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
#endif
#ifdef MIEEE
static unsigned short sincof[] = {
0x3de5,0xd8fd,0x1fd1,0x9ccd,
0xbe5a,0xe5e5,0xa929,0x1f5d,
0x3ec7,0x1de3,0x567d,0x48a1,
0xbf2a,0x01a0,0x19bf,0xdf03,
0x3f81,0x1111,0x1110,0xf7d0,
0xbfc5,0x5555,0x5555,0x5548,
};
static unsigned short coscof[24] = {
0xbda8,0xfa49,0xa086,0x1a9b,
0x3e21,0xee9d,0x7b4e,0x3f05,
0xbe92,0x7e4f,0x7eac,0x4bc6,
0x3efa,0x01a0,0x19c8,0x44f5,
0xbf56,0xc16c,0x16c1,0x4f91,
0x3fa5,0x5555,0x5555,0x554b,
};
static unsigned short P1[] = {0x3fe9,0x21fb,0x4000,0x0000};
static unsigned short P2[] = {0x3e64,0x442d,0x0000,0x0000};
static unsigned short P3[] = {0x3ce8,0x4698,0x98cc,0x5170};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double floor ( double );
extern double ldexp ( double, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double polevl(), floor(), ldexp();
int isnan(), isfinite();
#endif
extern double PIO4;
static double lossth = 1.073741824e9;
#ifdef NANS
extern double NAN;
#endif
#ifdef INFINITIES
extern double INFINITY;
#endif
double c_sin(x)
double x;
{
double y, z, zz;
int j, sign;
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
#ifdef NANS
if( isnan(x) )
return(x);
if( !isfinite(x) )
{
mtherr( "sin", DOMAIN );
return(NAN);
}
#endif
/* make argument positive but save the sign */
sign = 1;
if( x < 0 )
{
x = -x;
sign = -1;
}
if( x > lossth )
{
mtherr( "sin", TLOSS );
return(0.0);
}
y = floor( x/PIO4 ); /* integer part of x/PIO4 */
/* strip high bits of integer part to prevent integer overflow */
z = ldexp( y, -4 );
z = floor(z); /* integer part of y/8 */
z = y - ldexp( z, 4 ); /* y - 16 * (y/16) */
j = (int)z; /* convert to integer for tests on the phase angle */
/* map zeros to origin */
if( j & 1 )
{
j += 1;
y += 1.0;
}
j = j & 07; /* octant modulo 360 degrees */
/* reflect in x axis */
if( j > 3)
{
sign = -sign;
j -= 4;
}
/* Extended precision modular arithmetic */
z = ((x - y * DP1) - y * DP2) - y * DP3;
zz = z * z;
if( (j==1) || (j==2) )
{
y = 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
}
else
{
/* y = z + z * (zz * polevl( zz, sincof, 5 ));*/
y = z + z * z * z * polevl( zz, sincof, 5 );
}
if(sign < 0)
y = -y;
return(y);
}
double c_cos(x)
double x;
{
double y, z, zz;
int i;
int j, sign;
#ifdef NANS
if( isnan(x) )
return(x);
if( !isfinite(x) )
{
mtherr( "cos", DOMAIN );
return(NAN);
}
#endif
/* make argument positive */
sign = 1;
if( x < 0 )
x = -x;
if( x > lossth )
{
mtherr( "cos", TLOSS );
return(0.0);
}
y = floor( x/PIO4 );
z = ldexp( y, -4 );
z = floor(z); /* integer part of y/8 */
z = y - ldexp( z, 4 ); /* y - 16 * (y/16) */
/* integer and fractional part modulo one octant */
i = (int)z;
if( i & 1 ) /* map zeros to origin */
{
i += 1;
y += 1.0;
}
j = i & 07;
if( j > 3)
{
j -=4;
sign = -sign;
}
if( j > 1 )
sign = -sign;
/* Extended precision modular arithmetic */
z = ((x - y * DP1) - y * DP2) - y * DP3;
zz = z * z;
if( (j==1) || (j==2) )
{
/* y = z + z * (zz * polevl( zz, sincof, 5 ));*/
y = z + z * z * z * polevl( zz, sincof, 5 );
}
else
{
y = 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
}
if(sign < 0)
y = -y;
return(y);
}
/* Degrees, minutes, seconds to radians: */
/* 1 arc second, in radians = 4.8481368110953599358991410e-5 */
#ifdef DEC
static unsigned short P648[] = {034513,054170,0176773,0116043,};
#define P64800 *(double *)P648
#else
static double P64800 = 4.8481368110953599358991410e-5;
#endif
double radian(d,m,s)
double d,m,s;
{
return( ((d*60.0 + m)*60.0 + s)*P64800 );
}

148
src/utility/math/sinh.c Normal file
View file

@ -0,0 +1,148 @@
/* sinh.c
*
* Hyperbolic sine
*
*
*
* SYNOPSIS:
*
* double x, y, sinh();
*
* y = sinh( x );
*
*
*
* DESCRIPTION:
*
* Returns hyperbolic sine of argument in the range MINLOG to
* MAXLOG.
*
* The range is partitioned into two segments. If |x| <= 1, a
* rational function of the form x + x**3 P(x)/Q(x) is employed.
* Otherwise the calculation is sinh(x) = ( exp(x) - exp(-x) )/2.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC +- 88 50000 4.0e-17 7.7e-18
* IEEE +-MAXLOG 30000 2.6e-16 5.7e-17
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
static double P[] = {
-7.89474443963537015605E-1,
-1.63725857525983828727E2,
-1.15614435765005216044E4,
-3.51754964808151394800E5
};
static double Q[] = {
/* 1.00000000000000000000E0,*/
-2.77711081420602794433E2,
3.61578279834431989373E4,
-2.11052978884890840399E6
};
#endif
#ifdef DEC
static unsigned short P[] = {
0140112,0015377,0042731,0163255,
0142043,0134721,0146177,0123761,
0143464,0122706,0034353,0006017,
0144653,0140536,0157665,0054045
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0142212,0155404,0133513,0022040,
0044015,0036723,0173271,0011053,
0145400,0150407,0023710,0001034
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x3cd6,0xe8bb,0x435f,0xbfe9,
0xf4fe,0x398f,0x773a,0xc064,
0x6182,0xc71d,0x94b8,0xc0c6,
0xab05,0xdbf6,0x782b,0xc115
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x6484,0x96e9,0x5b60,0xc071,
0x2245,0x7ed7,0xa7ba,0x40e1,
0x0044,0xe4f9,0x1a20,0xc140
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0xbfe9,0x435f,0xe8bb,0x3cd6,
0xc064,0x773a,0x398f,0xf4fe,
0xc0c6,0x94b8,0xc71d,0x6182,
0xc115,0x782b,0xdbf6,0xab05
};
static unsigned short Q[] = {
0xc071,0x5b60,0x96e9,0x6484,
0x40e1,0xa7ba,0x7ed7,0x2245,
0xc140,0x1a20,0xe4f9,0x0044
};
#endif
#ifdef ANSIPROT
extern double fabs ( double );
extern double c_exp ( double );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
#else
double fabs(), c_exp(), polevl(), p1evl();
#endif
extern double INFINITY, MINLOG, MAXLOG, LOGE2;
double c_sinh(x)
double x;
{
double a;
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
a = fabs(x);
if( (x > (MAXLOG + LOGE2)) || (x > -(MINLOG-LOGE2) ) )
{
mtherr( "sinh", DOMAIN );
if( x > 0 )
return( INFINITY );
else
return( -INFINITY );
}
if( a > 1.0 )
{
if( a >= (MAXLOG - LOGE2) )
{
a = c_exp(0.5*a);
a = (0.5 * a) * a;
if( x < 0 )
a = -a;
return(a);
}
a = c_exp(a);
a = 0.5*a - (0.5/a);
if( x < 0 )
a = -a;
return(a);
}
a *= a;
return( x + x * a * (polevl(a,P,3)/p1evl(a,Q,3)) );
}

178
src/utility/math/sqrt.c Normal file
View file

@ -0,0 +1,178 @@
/* _sqrt.c
*
* Square root
*
*
*
* SYNOPSIS:
*
* double x, y, _sqrt();
*
* y = _sqrt( x );
*
*
*
* DESCRIPTION:
*
* Returns the square root of x.
*
* Range reduction involves isolating the power of two of the
* argument and using a polynomial approximation to obtain
* a rough value for the square root. Then Heron's iteration
* is used three times to converge to an accurate value.
*
*
*
* ACCURACY:
*
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0, 10 60000 2.1e-17 7.9e-18
* IEEE 0,1.7e308 30000 1.7e-16 6.3e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* _sqrt domain x < 0 0.0
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double frexp ( double, int * );
extern double ldexp ( double, int );
#else
double frexp(), ldexp();
#endif
extern double SQRT2; /* _sqrt2 = 1.41421356237309504880 */
double c_sqrt(x)
double x;
{
int e;
#ifndef UNK
short *q;
#endif
double z, w;
if( x <= 0.0 )
{
if( x < 0.0 )
mtherr( "_sqrt", DOMAIN );
return( 0.0 );
}
w = x;
/* separate exponent and significand */
#ifdef UNK
z = frexp( x, &e );
#endif
#ifdef DEC
q = (short *)&x;
e = ((*q >> 7) & 0377) - 0200;
*q &= 0177;
*q |= 040000;
z = x;
#endif
/* Note, frexp and ldexp are used in order to
* handle denormal numbers properly.
*/
#ifdef IBMPC
z = frexp( x, &e );
q = (short *)&x;
q += 3;
/*
e = ((*q >> 4) & 0x0fff) - 0x3fe;
*q &= 0x000f;
*q |= 0x3fe0;
z = x;
*/
#endif
#ifdef MIEEE
z = frexp( x, &e );
q = (short *)&x;
/*
e = ((*q >> 4) & 0x0fff) - 0x3fe;
*q &= 0x000f;
*q |= 0x3fe0;
z = x;
*/
#endif
/* approximate square root of number between 0.5 and 1
* relative error of approximation = 7.47e-3
*/
x = 4.173075996388649989089E-1 + 5.9016206709064458299663E-1 * z;
/* adjust for odd powers of 2 */
if( (e & 1) != 0 )
x *= SQRT2;
/* re-insert exponent */
#ifdef UNK
x = ldexp( x, (e >> 1) );
#endif
#ifdef DEC
*q += ((e >> 1) & 0377) << 7;
*q &= 077777;
#endif
#ifdef IBMPC
x = ldexp( x, (e >> 1) );
/*
*q += ((e >>1) & 0x7ff) << 4;
*q &= 077777;
*/
#endif
#ifdef MIEEE
x = ldexp( x, (e >> 1) );
/*
*q += ((e >>1) & 0x7ff) << 4;
*q &= 077777;
*/
#endif
/* Newton iterations: */
#ifdef UNK
x = 0.5*(x + w/x);
x = 0.5*(x + w/x);
x = 0.5*(x + w/x);
#endif
/* Note, assume the square root cannot be denormal,
* so it is safe to use integer exponent operations here.
*/
#ifdef DEC
x += w/x;
*q -= 0200;
x += w/x;
*q -= 0200;
x += w/x;
*q -= 0200;
#endif
#ifdef IBMPC
x += w/x;
*q -= 0x10;
x += w/x;
*q -= 0x10;
x += w/x;
*q -= 0x10;
#endif
#ifdef MIEEE
x += w/x;
*q -= 0x10;
x += w/x;
*q -= 0x10;
x += w/x;
*q -= 0x10;
#endif
return(x);
}

304
src/utility/math/tan.c Normal file
View file

@ -0,0 +1,304 @@
/* tan.c
*
* Circular tangent
*
*
*
* SYNOPSIS:
*
* double x, y, tan();
*
* y = tan( x );
*
*
*
* DESCRIPTION:
*
* Returns the circular tangent of the radian argument x.
*
* Range reduction is modulo pi/4. A rational function
* x + x**3 P(x**2)/Q(x**2)
* is employed in the basic interval [0, pi/4].
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC +-1.07e9 44000 4.1e-17 1.0e-17
* IEEE +-1.07e9 30000 2.9e-16 8.1e-17
*
* ERROR MESSAGES:
*
* message condition value returned
* tan total loss x > 1.073741824e9 0.0
*
*/
/* cot.c
*
* Circular cotangent
*
*
*
* SYNOPSIS:
*
* double x, y, cot();
*
* y = cot( x );
*
*
*
* DESCRIPTION:
*
* Returns the circular cotangent of the radian argument x.
*
* Range reduction is modulo pi/4. A rational function
* x + x**3 P(x**2)/Q(x**2)
* is employed in the basic interval [0, pi/4].
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE +-1.07e9 30000 2.9e-16 8.2e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* cot total loss x > 1.073741824e9 0.0
* cot singularity x = 0 INFINITY
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
yright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
static double P[] = {
-1.30936939181383777646E4,
1.15351664838587416140E6,
-1.79565251976484877988E7
};
static double Q[] = {
/* 1.00000000000000000000E0,*/
1.36812963470692954678E4,
-1.32089234440210967447E6,
2.50083801823357915839E7,
-5.38695755929454629881E7
};
static double DP1 = 7.853981554508209228515625E-1;
static double DP2 = 7.94662735614792836714E-9;
static double DP3 = 3.06161699786838294307E-17;
static double lossth = 1.073741824e9;
#endif
#ifdef DEC
static unsigned short P[] = {
0143514,0113306,0111171,0174674,
0045214,0147545,0027744,0167346,
0146210,0177526,0114514,0105660
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0043525,0142457,0072633,0025617,
0145241,0036742,0140525,0162256,
0046276,0146176,0013526,0143573,
0146515,0077401,0162762,0150607
};
/* 7.853981629014015197753906250000E-1 */
static unsigned short P1[] = {0040111,0007732,0120000,0000000,};
/* 4.960467869796758577649598009884E-10 */
static unsigned short P2[] = {0030410,0055060,0100000,0000000,};
/* 2.860594363054915898381331279295E-18 */
static unsigned short P3[] = {0021523,0011431,0105056,0001560,};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
static double lossth = 1.073741824e9;
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x3f38,0xd24f,0x92d8,0xc0c9,
0x9ddd,0xa5fc,0x99ec,0x4131,
0x9176,0xd329,0x1fea,0xc171
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x6572,0xeeb3,0xb8a5,0x40ca,
0xbc96,0x582a,0x27bc,0xc134,
0xd8ef,0xc2ea,0xd98f,0x4177,
0x5a31,0x3cbe,0xafe0,0xc189
};
/*
7.85398125648498535156E-1,
3.77489470793079817668E-8,
2.69515142907905952645E-15,
*/
static unsigned short P1[] = {0x0000,0x4000,0x21fb,0x3fe9};
static unsigned short P2[] = {0x0000,0x0000,0x442d,0x3e64};
static unsigned short P3[] = {0x5170,0x98cc,0x4698,0x3ce8};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
static double lossth = 1.073741824e9;
#endif
#ifdef MIEEE
static unsigned short P[] = {
0xc0c9,0x92d8,0xd24f,0x3f38,
0x4131,0x99ec,0xa5fc,0x9ddd,
0xc171,0x1fea,0xd329,0x9176
};
static unsigned short Q[] = {
0x40ca,0xb8a5,0xeeb3,0x6572,
0xc134,0x27bc,0x582a,0xbc96,
0x4177,0xd98f,0xc2ea,0xd8ef,
0xc189,0xafe0,0x3cbe,0x5a31
};
static unsigned short P1[] = {
0x3fe9,0x21fb,0x4000,0x0000
};
static unsigned short P2[] = {
0x3e64,0x442d,0x0000,0x0000
};
static unsigned short P3[] = {
0x3ce8,0x4698,0x98cc,0x5170,
};
#define DP1 *(double *)P1
#define DP2 *(double *)P2
#define DP3 *(double *)P3
static double lossth = 1.073741824e9;
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double floor ( double );
extern double ldexp ( double, int );
extern int isnan ( double );
extern int isfinite ( double );
static double tancot(double, int);
#else
double polevl(), p1evl(), floor(), ldexp();
static double tancot();
int isnan(), isfinite();
#endif
extern double PIO4;
extern double INFINITY;
extern double NAN;
double c_tan(x)
double x;
{
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
#ifdef NANS
if( isnan(x) )
return(x);
if( !isfinite(x) )
{
mtherr( "tan", DOMAIN );
return(NAN);
}
#endif
return( tancot(x,0) );
}
double c_cot(x)
double x;
{
if( x == 0.0 )
{
mtherr( "cot", SING );
return( INFINITY );
}
return( tancot(x,1) );
}
static double tancot( xx, cotflg )
double xx;
int cotflg;
{
double x, y, z, zz;
int j, sign;
/* make argument positive but save the sign */
if( xx < 0 )
{
x = -xx;
sign = -1;
}
else
{
x = xx;
sign = 1;
}
if( x > lossth )
{
if( cotflg )
mtherr( "cot", TLOSS );
else
mtherr( "tan", TLOSS );
return(0.0);
}
/* compute x mod PIO4 */
y = floor( x/PIO4 );
/* strip high bits of integer part */
z = ldexp( y, -3 );
z = floor(z); /* integer part of y/8 */
z = y - ldexp( z, 3 ); /* y - 16 * (y/16) */
/* integer and fractional part modulo one octant */
j = (int)z;
/* map zeros and singularities to origin */
if( j & 1 )
{
j += 1;
y += 1.0;
}
z = ((x - y * DP1) - y * DP2) - y * DP3;
zz = z * z;
if( zz > 1.0e-14 )
y = z + z * (zz * polevl( zz, P, 2 )/p1evl(zz, Q, 4));
else
y = z;
if( j & 2 )
{
if( cotflg )
y = -y;
else
y = -1.0/y;
}
else
{
if( cotflg )
y = 1.0/y;
}
if( sign < 0 )
y = -y;
return( y );
}

141
src/utility/math/tanh.c Normal file
View file

@ -0,0 +1,141 @@
/* tanh.c
*
* Hyperbolic tangent
*
*
*
* SYNOPSIS:
*
* double x, y, tanh();
*
* y = tanh( x );
*
*
*
* DESCRIPTION:
*
* Returns hyperbolic tangent of argument in the range MINLOG to
* MAXLOG.
*
* A rational function is used for |x| < 0.625. The form
* x + x**3 P(x)/Q(x) of Cody _& Waite is employed.
* Otherwise,
* tanh(x) = sinh(x)/cosh(x) = 1 - 2/(exp(2x) + 1).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -2,2 50000 3.3e-17 6.4e-18
* IEEE -2,2 30000 2.5e-16 5.8e-17
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
static double P[] = {
-9.64399179425052238628E-1,
-9.92877231001918586564E1,
-1.61468768441708447952E3
};
static double Q[] = {
/* 1.00000000000000000000E0,*/
1.12811678491632931402E2,
2.23548839060100448583E3,
4.84406305325125486048E3
};
#endif
#ifdef DEC
static unsigned short P[] = {
0140166,0161335,0053753,0075126,
0141706,0111520,0070463,0040552,
0142711,0153001,0101300,0025430
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041741,0117624,0051300,0156060,
0043013,0133720,0071251,0127717,
0043227,0060201,0021020,0020136
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x6f4b,0xaafd,0xdc5b,0xbfee,
0x682d,0x0e26,0xd26a,0xc058,
0x0563,0x3058,0x3ac0,0xc099
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x1b86,0x8a58,0x33f2,0x405c,
0x35fa,0x0e55,0x76fa,0x40a1,
0x040c,0x2442,0xec10,0x40b2
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0xbfee,0xdc5b,0xaafd,0x6f4b,
0xc058,0xd26a,0x0e26,0x682d,
0xc099,0x3ac0,0x3058,0x0563
};
static unsigned short Q[] = {
0x405c,0x33f2,0x8a58,0x1b86,
0x40a1,0x76fa,0x0e55,0x35fa,
0x40b2,0xec10,0x2442,0x040c
};
#endif
#ifdef ANSIPROT
extern double fabs ( double );
extern double c_exp ( double );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
#else
double fabs(), c_exp(), polevl(), p1evl();
#endif
extern double MAXLOG;
double c_tanh(x)
double x;
{
double s, z;
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
z = fabs(x);
if( z > 0.5 * MAXLOG )
{
if( x > 0 )
return( 1.0 );
else
return( -1.0 );
}
if( z >= 0.625 )
{
s = c_exp(2.0*z);
z = 1.0 - 2.0/(s + 1.0);
if( x < 0 )
z = -z;
}
else
{
if( x == 0.0 )
return(x);
s = x * x;
z = polevl( s, P, 2 )/p1evl(s, Q, 3);
z = x * s * z;
z = x + z;
}
return( z );
}