zdoom/qzdoom
0
0
Fork 0
mirror of https://github.com/ZDoom/qzdoom.git synced 2025-01-18 07:02:03 +00:00
Code Issues Projects Releases Packages Wiki Activity

- added custom math routines for reliability.

Browse source
This commit is contained in:
Christoph Oelckers 2016-03-11 15:45:47 +01:00
parent 7edd5e2dac
commit 161d03231a
50 changed files with 4025 additions and 120 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
CMakeLists.txt
View file

@ -164,7 +164,7 @@ if( MSVC )
string(REPLACE " /GR" " " CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS} )
else()
set( REL_LINKER_FLAGS "" )
set( ALL_C_FLAGS "" )
set( ALL_C_FLAGS "-ffp-contract=off" )
set( REL_C_FLAGS "" )
set( DEB_C_FLAGS "" )
endif()

16
src/CMakeLists.txt
View file

@ -1159,6 +1159,22 @@ add_executable( zdoom WIN32 MACOSX_BUNDLE
g_shared/sbar_mugshot.cpp
g_shared/shared_hud.cpp
g_shared/shared_sbar.cpp
math/asin.c
math/atan.c
math/const.c
math/cosh.c
math/exp.c
math/isnan.c
math/log.c
math/log10.c
math/mtherr.c
math/polevl.c
math/sin.c
math/sinh.c
math/sqrt.c
math/tan.c
math/tanh.c
math/fastsin.cpp
resourcefiles/ancientzip.cpp
resourcefiles/file_7z.cpp
resourcefiles/file_grp.cpp

3
src/b_move.cpp
View file

@ -21,6 +21,7 @@
#include "d_player.h"
#include "p_spec.h"
#include "p_checkposition.h"
#include "math/cmath.h"
static FRandom pr_botopendoor ("BotOpenDoor");
static FRandom pr_bottrywalk ("BotTryWalk");
@ -348,7 +349,7 @@ void DBot::Pitch (AActor *target)
double diff;
diff = target->Z() - player->mo->Z();
aim = atan(diff / (double)player->mo->AproxDistance(target));
aim = g_atan(diff / (double)player->mo->AproxDistance(target));
player->mo->pitch = -(int)(aim * ANGLE_180/M_PI);
}

23
src/fragglescript/t_func.cpp
View file

@ -69,6 +69,7 @@
#include "farchive.h"
#include "p_setup.h"
#include "p_spec.h"
#include "math/cmath.h"
static FRandom pr_script("FScript");
@ -1431,7 +1432,7 @@ void FParser::SF_PointToDist(void)
double y = floatvalue(t_argv[3]) - floatvalue(t_argv[1]);
t_return.type = svt_fixed;
t_return.value.f = FLOAT2FIXED(sqrt(x*x+y*y));
t_return.value.f = FLOAT2FIXED(g_sqrt(x*x+y*y));
}
}
@ -3830,7 +3831,7 @@ void FParser::SF_Sin()
if (CheckArgs(1))
{
t_return.type = svt_fixed;
t_return.value.f = FLOAT2FIXED(sin(floatvalue(t_argv[0])));
t_return.value.f = FLOAT2FIXED(g_sin(floatvalue(t_argv[0])));
}
}
@ -3840,7 +3841,7 @@ void FParser::SF_ASin()
if (CheckArgs(1))
{
t_return.type = svt_fixed;
t_return.value.f = FLOAT2FIXED(asin(floatvalue(t_argv[0])));
t_return.value.f = FLOAT2FIXED(g_asin(floatvalue(t_argv[0])));
}
}
@ -3850,7 +3851,7 @@ void FParser::SF_Cos()
if (CheckArgs(1))
{
t_return.type = svt_fixed;
t_return.value.f = FLOAT2FIXED(cos(floatvalue(t_argv[0])));
t_return.value.f = FLOAT2FIXED(g_cos(floatvalue(t_argv[0])));
}
}
@ -3860,7 +3861,7 @@ void FParser::SF_ACos()
if (CheckArgs(1))
{
t_return.type = svt_fixed;
t_return.value.f = FLOAT2FIXED(acos(floatvalue(t_argv[0])));
t_return.value.f = FLOAT2FIXED(g_acos(floatvalue(t_argv[0])));
}
}
@ -3870,7 +3871,7 @@ void FParser::SF_Tan()
if (CheckArgs(1))
{
t_return.type = svt_fixed;
t_return.value.f = FLOAT2FIXED(tan(floatvalue(t_argv[0])));
t_return.value.f = FLOAT2FIXED(g_tan(floatvalue(t_argv[0])));
}
}
@ -3880,7 +3881,7 @@ void FParser::SF_ATan()
if (CheckArgs(1))
{
t_return.type = svt_fixed;
t_return.value.f = FLOAT2FIXED(atan(floatvalue(t_argv[0])));
t_return.value.f = FLOAT2FIXED(g_atan(floatvalue(t_argv[0])));
}
}
@ -3890,7 +3891,7 @@ void FParser::SF_Exp()
if (CheckArgs(1))
{
t_return.type = svt_fixed;
t_return.value.f = FLOAT2FIXED(exp(floatvalue(t_argv[0])));
t_return.value.f = FLOAT2FIXED(g_exp(floatvalue(t_argv[0])));
}
}
@ -3899,7 +3900,7 @@ void FParser::SF_Log()
if (CheckArgs(1))
{
t_return.type = svt_fixed;
t_return.value.f = FLOAT2FIXED(log(floatvalue(t_argv[0])));
t_return.value.f = FLOAT2FIXED(g_log(floatvalue(t_argv[0])));
}
}
@ -3909,7 +3910,7 @@ void FParser::SF_Sqrt()
if (CheckArgs(1))
{
t_return.type = svt_fixed;
t_return.value.f = FLOAT2FIXED(sqrt(floatvalue(t_argv[0])));
t_return.value.f = FLOAT2FIXED(g_sqrt(floatvalue(t_argv[0])));
}
}
@ -4043,7 +4044,7 @@ void FParser::SF_SetCorona(void)
break;
case 6:
lspr[num].dynamic_radius = fval;
lspr[num].dynamic_sqrradius = sqrt(lspr[num].dynamic_radius);
lspr[num].dynamic_sqrradius = g_sqrt(lspr[num].dynamic_radius);
break;
default:
CONS_Printf("Error in setcorona\n");

2
src/g_heretic/a_hereticweaps.cpp
View file

@ -484,7 +484,7 @@ DEFINE_ACTION_FUNCTION(AActor, A_MacePL1Check)
self->velx = FixedMul(7*FRACUNIT, finecosine[angle]);
self->vely = FixedMul(7*FRACUNIT, finesine[angle]);
#else
double velscale = sqrt ((double)self->velx * (double)self->velx +
double velscale = g_sqrt ((double)self->velx * (double)self->velx +
(double)self->vely * (double)self->vely);
velscale = 458752 / velscale;
self->velx = (int)(self->velx * velscale);

2
src/g_hexen/a_flechette.cpp
View file

@ -117,7 +117,7 @@ bool AArtiPoisonBag3::Use (bool pickup)
angle_t orgpitch = angle_t(-Owner->pitch) >> ANGLETOFINESHIFT;
angle_t modpitch = angle_t(0xDC00000 - Owner->pitch) >> ANGLETOFINESHIFT;
angle_t angle = mo->angle >> ANGLETOFINESHIFT;
fixed_t speed = fixed_t(sqrt((double)mo->Speed*mo->Speed + (4.0*65536*4*65536)));
fixed_t speed = fixed_t(g_sqrt((double)mo->Speed*mo->Speed + (4.0*65536*4*65536)));
fixed_t xyscale = FixedMul(speed, finecosine[modpitch]);
mo->velz = FixedMul(speed, finesine[modpitch]);

3
src/g_shared/a_camera.cpp
View file

@ -37,6 +37,7 @@
#include "a_sharedglobal.h"
#include "p_local.h"
#include "farchive.h"
#include "math/cmath.h"
/*
== SecurityCamera
@ -180,7 +181,7 @@ void AAimingCamera::Tick ()
DVector2 vect(fv3.x, fv3.y);
double dz = Z() - tracer->Z() - tracer->height/2;
double dist = vect.Length();
double ang = dist != 0.f ? atan2 (dz, dist) : 0;
double ang = dist != 0.f ? g_atan2 (dz, dist) : 0;
int desiredpitch = (int)RAD2ANGLE(ang);
if (abs (desiredpitch - pitch) < MaxPitchChange)
{

2
src/g_shared/a_decals.cpp
View file

@ -414,7 +414,7 @@ static void GetWallStuff (side_t *wall, vertex_t *&v1, fixed_t &ldx, fixed_t &ld
static fixed_t Length (fixed_t dx, fixed_t dy)
{
return (fixed_t)sqrt ((double)dx*(double)dx+(double)dy*(double)dy);
return (fixed_t)g_sqrt ((double)dx*(double)dx+(double)dy*(double)dy);
}
static side_t *NextWall (const side_t *wall)

8
src/g_shared/a_movingcamera.cpp
View file

@ -433,8 +433,8 @@ bool APathFollower::Interpolate ()
double fdx = FIXED2DBL(dx);
double fdy = FIXED2DBL(dy);
double fdz = FIXED2DBL(-dz);
double dist = sqrt (fdx*fdx + fdy*fdy);
double ang = dist != 0.f ? atan2 (fdz, dist) : 0;
double dist = g_sqrt (fdx*fdx + fdy*fdy);
double ang = dist != 0.f ? g_atan2 (fdz, dist) : 0;
pitch = (fixed_t)RAD2ANGLE(ang);
}
}
@ -672,8 +672,8 @@ bool AMovingCamera::Interpolate ()
double dx = FIXED2DBL(X() - tracer->X());
double dy = FIXED2DBL(Y() - tracer->Y());
double dz = FIXED2DBL(Z() - tracer->Z() - tracer->height/2);
double dist = sqrt (dx*dx + dy*dy);
double ang = dist != 0.f ? atan2 (dz, dist) : 0;
double dist = g_sqrt (dx*dx + dy*dy);
double ang = dist != 0.f ? g_atan2 (dz, dist) : 0;
pitch = RAD2ANGLE(ang);
}

2
src/g_shared/a_quake.cpp
View file

@ -175,7 +175,7 @@ fixed_t DEarthquake::GetModWave(double waveMultiplier) const
//Named waveMultiplier because that's as the name implies: adds more waves per second.
double time = m_Countdown - FIXED2DBL(r_TicFrac);
return FLOAT2FIXED(sin(waveMultiplier * time * (M_PI * 2 / TICRATE)));
return FLOAT2FIXED(g_sin(waveMultiplier * time * (M_PI * 2 / TICRATE)));
}
//==========================================================================

315
src/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/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 );
}

139
src/math/cmath.h Normal file
View file

@ -0,0 +1,139 @@
#ifndef __CMATH_H
#define __CMATH_H
#include "tables.h"
#include "m_fixed.h"
#define USE_CUSTOM_MATH // we want reoreducably 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);
}
// 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;
inline double fastcosdeg(double v)
{
return fasttrig.cos(FLOAT2ANGLE(v));
}
inline double fastsindeg(double v)
{
return fasttrig.sin(FLOAT2ANGLE(v));
}
inline double fastcos(double v)
{
return fasttrig.cos(RAD2ANGLE(v));
}
inline double fastsin(double v)
{
return fasttrig.sin(RAD2ANGLE(v));
}
inline double sindeg(double v)
{
#ifdef USE_CUSTOM_MATH
return c_sin(v * (M_PI / 180.));
#else
return sin(v * (M_PI / 180.));
#endif
}
inline double cosdeg(double v)
{
#ifdef USE_CUSTOM_MATH
return c_cos(v * (M_PI / 180.));
#else
return cos(v * (M_PI / 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
#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
#endif
#endif

252
src/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/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/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);
}

102
src/math/fastsin.cpp Normal file
View file

@ -0,0 +1,102 @@
/*
** 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"
#include "m_fixed.h"
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 - 1)) == 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/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/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/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/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/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/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 );
}

15
src/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/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/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/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/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/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 );
}

2
src/nodebuild_events.cpp
View file

@ -219,7 +219,7 @@ void FEventTree::PrintTree (const FEvent *event) const
{
PrintTree(event->Left);
sprintf(buff, " Distance %g, vertex %d, seg %u\n",
sqrt(event->Distance/4294967296.0), event->Info.Vertex, (unsigned)event->Info.FrontSeg);
g_sqrt(event->Distance/4294967296.0), event->Info.Vertex, (unsigned)event->Info.FrontSeg);
Printf(PRINT_LOG, "%s", buff);
PrintTree(event->Right);
}

3
src/nodebuild_utility.cpp
View file

@ -50,6 +50,7 @@
#include "i_system.h"
#include "po_man.h"
#include "r_state.h"
#include "math/cmath.h"
static const int PO_LINE_START = 1;
static const int PO_LINE_EXPLICIT = 5;
@ -74,7 +75,7 @@ angle_t FNodeBuilder::PointToAngle (fixed_t x, fixed_t y)
// See https://gcc.gnu.org/wiki/Math_Optimization_Flags for details
long double ang = atan2l (double(y), double(x));
#else // !__APPLE__ || __llvm__
double ang = atan2 (double(y), double(x));
double ang = g_atan2 (double(y), double(x));
#endif // __APPLE__ && !__llvm__
// Convert to signed first since negative double to unsigned is undefined.
return angle_t(int(ang * rad2bam)) << 1;

4
src/p_acs.cpp
View file

@ -5335,10 +5335,10 @@ int DLevelScript::CallFunction(int argCount, int funcIndex, SDWORD *args)
return P_IsTIDUsed(args[0]);
case ACSF_Sqrt:
return xs_FloorToInt(sqrt(double(args[0])));
return xs_FloorToInt(g_sqrt(double(args[0])));
case ACSF_FixedSqrt:
return FLOAT2FIXED(sqrt(FIXED2DBL(args[0])));
return FLOAT2FIXED(g_sqrt(FIXED2DBL(args[0])));
case ACSF_VectorLength:
return FLOAT2FIXED(DVector2(FIXED2DBL(args[0]), FIXED2DBL(args[1])).Length());

2
src/p_buildmap.cpp
View file

@ -639,7 +639,7 @@ static void LoadWalls (walltype *walls, int numwalls, sectortype *bsec)
slope.wal2 = &walls[slope.wal->point2];
slope.dx = slope.wal2->x - slope.wal->x;
slope.dy = slope.wal2->y - slope.wal->y;
slope.i = long (sqrt ((double)(slope.dx*slope.dx+slope.dy*slope.dy))) << 5;
slope.i = long (g_sqrt ((double)(slope.dx*slope.dx+slope.dy*slope.dy))) << 5;
if (slope.i == 0)
{
continue;

2
src/p_effect.cpp
View file

@ -662,7 +662,7 @@ void P_DrawRailTrail(AActor *source, const DVector3 &start, const DVector3 &end,
dir = end - start;
lengthsquared = dir | dir;
length = sqrt(lengthsquared);
length = g_sqrt(lengthsquared);
steps = xs_FloorToInt(length / 3);
fullbright = !!(flags & RAF_FULLBRIGHT);

7
src/p_enemy.cpp
View file

@ -52,6 +52,7 @@
#include "teaminfo.h"
#include "p_spec.h"
#include "p_checkposition.h"
#include "math/cmath.h"
#include "gi.h"
@ -2896,8 +2897,8 @@ void A_Face (AActor *self, AActor *other, angle_t max_turn, angle_t max_pitch, a
target_z += z_add;
double dist_z = target_z - source_z;
double ddist = sqrt(dist.X*dist.X + dist.Y*dist.Y + dist_z*dist_z);
int other_pitch = (int)RAD2ANGLE(asin(dist_z / ddist));
double ddist = g_sqrt(dist.X*dist.X + dist.Y*dist.Y + dist_z*dist_z);
int other_pitch = (int)RAD2ANGLE(g_asin(dist_z / ddist));
if (max_pitch != 0)
{
@ -3009,7 +3010,7 @@ DEFINE_ACTION_FUNCTION(AActor, A_MonsterRail)
fixedvec2 pos = self->Vec2To(self->target);
DVector2 xydiff(pos.x, pos.y);
double zdiff = (self->target->Z() + (self->target->height>>1)) - (self->Z() + (self->height>>1) - self->floorclip);
self->pitch = int(atan2(zdiff, xydiff.Length()) * ANGLE_180 / -M_PI);
self->pitch = int(g_atan2(zdiff, xydiff.Length()) * ANGLE_180 / -M_PI);
}
// Let the aim trail behind the player

2
src/p_glnodes.cpp
View file

@ -1388,7 +1388,7 @@ static bool PointOnLine (int x, int y, int x1, int y1, int dx, int dy)
// Either the point is very near the line, or the segment defining
// the line is very short: Do a more expensive test to determine
// just how far from the line the point is.
double l = sqrt(d_dx*d_dx+d_dy*d_dy);
double l = g_sqrt(d_dx*d_dx+d_dy*d_dy);
double dist = fabs(s_num)/l;
if (dist < SIDE_EPSILON)
{

7
src/p_map.cpp
View file

@ -31,6 +31,7 @@
#include "m_random.h"
#include "i_system.h"
#include "c_dispatch.h"
#include "math/cmath.h"
#include "doomdef.h"
#include "p_local.h"
@ -3309,7 +3310,7 @@ bool FSlide::BounceWall(AActor *mo)
deltaangle >>= ANGLETOFINESHIFT;
movelen = fixed_t(sqrt(double(mo->velx)*mo->velx + double(mo->vely)*mo->vely));
movelen = fixed_t(g_sqrt(double(mo->velx)*mo->velx + double(mo->vely)*mo->vely));
movelen = FixedMul(movelen, mo->wallbouncefactor);
FBoundingBox box(mo->X(), mo->Y(), mo->radius);
@ -4589,7 +4590,7 @@ void P_TraceBleed(int damage, AActor *target, AActor *missile)
{
double aim;
aim = atan((double)missile->velz / (double)target->AproxDistance(missile));
aim = g_atan((double)missile->velz / (double)target->AproxDistance(missile));
pitch = -(int)(aim * ANGLE_180 / PI);
}
else
@ -5336,7 +5337,7 @@ void P_RadiusAttack(AActor *bombspot, AActor *bombsource, int bombdamage, int bo
else
{
len -= boxradius;
len = sqrt(len*len + dz*dz);
len = g_sqrt(len*len + dz*dz);
}
}
else

2
src/p_maputl.cpp
View file

@ -388,7 +388,7 @@ bool AActor::FixMapthingPos()
// Get the exact distance to the line
divline_t dll, dlv;
fixed_t linelen = (fixed_t)sqrt((double)ldef->dx*ldef->dx + (double)ldef->dy*ldef->dy);
fixed_t linelen = (fixed_t)g_sqrt((double)ldef->dx*ldef->dx + (double)ldef->dy*ldef->dy);
P_MakeDivline(ldef, &dll);

2
src/p_mobj.cpp
View file

@ -6775,6 +6775,6 @@ void PrintMiscActorInfo(AActor *query)
FIXED2DBL(query->floorz), FIXED2DBL(query->ceilingz));
Printf("\nSpeed= %f, velocity= x:%f, y:%f, z:%f, combined:%f.\n",
FIXED2DBL(query->Speed), FIXED2DBL(query->velx), FIXED2DBL(query->vely), FIXED2DBL(query->velz),
sqrt(pow(FIXED2DBL(query->velx), 2) + pow(FIXED2DBL(query->vely), 2) + pow(FIXED2DBL(query->velz), 2)));
g_sqrt(pow(FIXED2DBL(query->velx), 2) + pow(FIXED2DBL(query->vely), 2) + pow(FIXED2DBL(query->velz), 2)));
}
}

4
src/p_setup.cpp
View file

@ -1324,7 +1324,7 @@ void P_LoadSegs (MapData * map)
segangle >>= (ANGLETOFINESHIFT-16);
dx = (li->v1->x - li->v2->x)>>FRACBITS;
dy = (li->v1->y - li->v2->y)>>FRACBITS;
dis = ((int) sqrt((double)(dx*dx + dy*dy)))<<FRACBITS;
dis = ((int) g_sqrt((double)(dx*dx + dy*dy)))<<FRACBITS;
dx = finecosine[segangle];
dy = finesine[segangle];
if ((vnum2 > vnum1) && (vertchanged[vnum2] == 0))
@ -2024,7 +2024,7 @@ void P_FinishLoadingLineDef(line_t *ld, int alpha)
}
// [RH] Set some new sidedef properties
int len = (int)(sqrt (dx*dx + dy*dy) + 0.5f);
int len = (int)(g_sqrt (dx*dx + dy*dy) + 0.5f);
if (ld->sidedef[0] != NULL)
{

9
src/p_things.cpp
View file

@ -50,6 +50,7 @@
#include "d_player.h"
#include "r_utility.h"
#include "p_spec.h"
#include "math/cmath.h"
// Set of spawnable things for the Thing_Spawn and Thing_Projectile specials.
FClassMap SpawnableThings;
@ -283,9 +284,9 @@ bool P_Thing_Projectile (int tid, AActor *source, int type, const char *type_nam
double dist = aim.Length();
double targspeed = tvel.Length();
double ydotx = -aim | tvel;
double a = acos (clamp (ydotx / targspeed / dist, -1.0, 1.0));
double a = g_acos (clamp (ydotx / targspeed / dist, -1.0, 1.0));
double multiplier = double(pr_leadtarget.Random2())*0.1/255+1.1;
double sinb = -clamp (targspeed*multiplier * sin(a) / fspeed, -1.0, 1.0);
double sinb = -clamp (targspeed*multiplier * g_sin(a) / fspeed, -1.0, 1.0);
// Use the cross product of two of the triangle's sides to get a
// rotation vector.
@ -293,7 +294,7 @@ bool P_Thing_Projectile (int tid, AActor *source, int type, const char *type_nam
// The vector must be normalized.
rv.MakeUnit();
// Now combine the rotation vector with angle b to get a rotation matrix.
DMatrix3x3 rm(rv, cos(asin(sinb)), sinb);
DMatrix3x3 rm(rv, g_cos(g_asin(sinb)), sinb);
// And multiply the original aim vector with the matrix to get a
// new aim vector that leads the target.
DVector3 aimvec = rm * aim;
@ -328,7 +329,7 @@ nolead:
// Set the missile's speed to reflect the speed it was spawned at.
if (mobj->flags & MF_MISSILE)
{
mobj->Speed = fixed_t (sqrt (double(mobj->velx)*mobj->velx + double(mobj->vely)*mobj->vely + double(mobj->velz)*mobj->velz));
mobj->Speed = fixed_t (g_sqrt (double(mobj->velx)*mobj->velx + double(mobj->vely)*mobj->vely + double(mobj->velz)*mobj->velz));
}
// Hugger missiles don't have any vertical velocity
if (mobj->flags3 & (MF3_FLOORHUGGER|MF3_CEILINGHUGGER))

2
src/po_man.cpp
View file

@ -1778,7 +1778,7 @@ void PO_Init (void)
node_t *no = &nodes[i];
double fdx = (double)no->dx;
double fdy = (double)no->dy;
no->len = (float)sqrt(fdx * fdx + fdy * fdy);
no->len = (float)g_sqrt(fdx * fdx + fdy * fdy);
}
// mark all subsectors which have a seg belonging to a polyobj

9
src/portal.cpp
View file

@ -53,6 +53,7 @@
#include "p_maputl.h"
#include "p_spec.h"
#include "p_checkposition.h"
#include "math/cmath.h"
// simulation recurions maximum
CVAR(Int, sv_portal_recursions, 4, CVAR_ARCHIVE|CVAR_SERVERINFO)
@ -247,9 +248,9 @@ static void SetRotation(FLinePortal *port)
{
line_t *dst = port->mDestination;
line_t *line = port->mOrigin;
double angle = atan2(dst->dy, dst->dx) - atan2(line->dy, line->dx) + M_PI;
port->mSinRot = FLOAT2FIXED(sin(angle));
port->mCosRot = FLOAT2FIXED(cos(angle));
double angle = g_atan2(dst->dy, dst->dx) - g_atan2(line->dy, line->dx) + M_PI;
port->mSinRot = FLOAT2FIXED(g_sin(angle));
port->mCosRot = FLOAT2FIXED(g_cos(angle));
port->mAngleDiff = RAD2ANGLE(angle);
}
@ -665,7 +666,7 @@ void P_NormalizeVXVY(fixed_t& vx, fixed_t& vy)
{
double _vx = FIXED2DBL(vx);
double _vy = FIXED2DBL(vy);
double len = sqrt(_vx*_vx+_vy*_vy);
double len = g_sqrt(_vx*_vx+_vy*_vy);
vx = FLOAT2FIXED(_vx/len);
vy = FLOAT2FIXED(_vy/len);
}

11
src/r_utility.cpp
View file

@ -57,6 +57,7 @@
#include "r_utility.h"
#include "d_player.h"
#include "p_local.h"
#include "math/cmath.h"
// EXTERNAL DATA DECLARATIONS ----------------------------------------------
@ -254,7 +255,7 @@ angle_t R_PointToAngle2 (fixed_t x1, fixed_t y1, fixed_t x, fixed_t y)
else
{
// we have to use the slower but more precise floating point atan2 function here.
return xs_RoundToUInt(atan2(double(y), double(x)) * (ANGLE_180/M_PI));
return xs_RoundToUInt(g_atan2(double(y), double(x)) * (ANGLE_180/M_PI));
}
}
@ -273,7 +274,7 @@ void R_InitPointToAngle (void)
//
for (i = 0; i <= SLOPERANGE; i++)
{
f = atan2 ((double)i, (double)SLOPERANGE) / (6.28318530718 /* 2*pi */);
f = g_atan2 ((double)i, (double)SLOPERANGE) / (6.28318530718 /* 2*pi */);
tantoangle[i] = (angle_t)(0xffffffff*f);
}
}
@ -322,16 +323,16 @@ void R_InitTables (void)
const double pimul = PI*2/FINEANGLES;
// viewangle tangent table
finetangent[0] = (fixed_t)(FRACUNIT*tan ((0.5-FINEANGLES/4)*pimul)+0.5);
finetangent[0] = (fixed_t)(FRACUNIT*g_tan ((0.5-FINEANGLES/4)*pimul)+0.5);
for (i = 1; i < FINEANGLES/2; i++)
{
finetangent[i] = (fixed_t)(FRACUNIT*tan ((i-FINEANGLES/4)*pimul)+0.5);
finetangent[i] = (fixed_t)(FRACUNIT*g_tan ((i-FINEANGLES/4)*pimul)+0.5);
}
// finesine table
for (i = 0; i < FINEANGLES/4; i++)
{
finesine[i] = (fixed_t)(FRACUNIT * sin (i*pimul));
finesine[i] = (fixed_t)(FRACUNIT * g_sin (i*pimul));
}
for (i = 0; i < FINEANGLES/4; i++)
{

2
src/r_utility.h
View file

@ -60,7 +60,7 @@ angle_t R_PointToAngle2 (fixed_t x1, fixed_t y1, fixed_t x2, fixed_t y2);
inline angle_t R_PointToAngle (fixed_t x, fixed_t y) { return R_PointToAngle2 (viewx, viewy, x, y); }
inline angle_t R_PointToAnglePrecise (fixed_t viewx, fixed_t viewy, fixed_t x, fixed_t y)
{
return xs_RoundToUInt(atan2(double(y-viewy), double(x-viewx)) * (ANGLE_180/M_PI));
return xs_RoundToUInt(g_atan2(double(y-viewy), double(x-viewx)) * (ANGLE_180/M_PI));
}
// Used for interpolation waypoints.

3
src/thingdef/thingdef_codeptr.cpp
View file

@ -76,6 +76,7 @@
#include "d_player.h"
#include "p_maputl.h"
#include "p_spec.h"
#include "math/cmath.h"
AActor *SingleActorFromTID(int tid, AActor *defactor);
@ -1935,7 +1936,7 @@ DEFINE_ACTION_FUNCTION_PARAMS(AActor, A_CustomRailgun)
DVector2 xydiff(pos.x, pos.y);
double zdiff = (self->target->Z() + (self->target->height>>1)) -
(self->Z() + (self->height>>1) - self->floorclip);
self->pitch = int(atan2(zdiff, xydiff.Length()) * ANGLE_180 / -M_PI);
self->pitch = int(g_atan2(zdiff, xydiff.Length()) * ANGLE_180 / -M_PI);
}
// Let the aim trail behind the player
if (aim)

27
src/thingdef/thingdef_expression.cpp
View file

@ -53,6 +53,7 @@
#include "m_fixed.h"
#include "vmbuilder.h"
#include "v_text.h"
#include "math/cmath.h"
struct FLOP
{
@ -65,23 +66,23 @@ struct FLOP
// degrees to radians for those that work with angles.
static const FLOP FxFlops[] =
{
{ NAME_Exp, FLOP_EXP, [](double v) { return exp(v); } },
{ NAME_Log, FLOP_LOG, [](double v) { return log(v); } },
{ NAME_Log10, FLOP_LOG10, [](double v) { return log10(v); } },
{ NAME_Sqrt, FLOP_SQRT, [](double v) { return sqrt(v); } },
{ NAME_Exp, FLOP_EXP, [](double v) { return g_exp(v); } },
{ NAME_Log, FLOP_LOG, [](double v) { return g_log(v); } },
{ NAME_Log10, FLOP_LOG10, [](double v) { return g_log10(v); } },
{ NAME_Sqrt, FLOP_SQRT, [](double v) { return g_sqrt(v); } },
{ NAME_Ceil, FLOP_CEIL, [](double v) { return ceil(v); } },
{ NAME_Floor, FLOP_FLOOR, [](double v) { return floor(v); } },
{ NAME_ACos, FLOP_ACOS_DEG, [](double v) { return acos(v) * (180.0 / M_PI); } },
{ NAME_ASin, FLOP_ASIN_DEG, [](double v) { return asin(v) * (180.0 / M_PI); } },
{ NAME_ATan, FLOP_ATAN_DEG, [](double v) { return atan(v) * (180.0 / M_PI); } },
{ NAME_Cos, FLOP_COS_DEG, [](double v) { return cos(v * (M_PI / 180.0)); } },
{ NAME_Sin, FLOP_SIN_DEG, [](double v) { return sin(v * (M_PI / 180.0)); } },
{ NAME_Tan, FLOP_TAN_DEG, [](double v) { return tan(v * (M_PI / 180.0)); } },
{ NAME_ACos, FLOP_ACOS_DEG, [](double v) { return g_acos(v) * (180.0 / M_PI); } },
{ NAME_ASin, FLOP_ASIN_DEG, [](double v) { return g_asin(v) * (180.0 / M_PI); } },
{ NAME_ATan, FLOP_ATAN_DEG, [](double v) { return g_atan(v) * (180.0 / M_PI); } },
{ NAME_Cos, FLOP_COS_DEG, [](double v) { return g_cosdeg(v); } },
{ NAME_Sin, FLOP_SIN_DEG, [](double v) { return g_sindeg(v); } },
{ NAME_Tan, FLOP_TAN_DEG, [](double v) { return g_tan(v * (M_PI / 180.0)); } },
{ NAME_CosH, FLOP_COSH, [](double v) { return cosh(v); } },
{ NAME_SinH, FLOP_SINH, [](double v) { return sinh(v); } },
{ NAME_TanH, FLOP_TANH, [](double v) { return tanh(v); } },
{ NAME_CosH, FLOP_COSH, [](double v) { return g_cosh(v); } },
{ NAME_SinH, FLOP_SINH, [](double v) { return g_sinh(v); } },
{ NAME_TanH, FLOP_TANH, [](double v) { return g_tanh(v); } },
};
ExpEmit::ExpEmit(VMFunctionBuilder *build, int type)

76
src/vectors.h
View file

@ -43,6 +43,7 @@
#include <math.h>
#include <string.h>
#include "m_fixed.h"
#include "math/cmath.h"
#define EQUAL_EPSILON (1/65536.f)
@ -227,7 +228,7 @@ struct TVector2
// Vector length
vec_t Length() const
{
return (vec_t)sqrt (X*X + Y*Y);
return (vec_t)g_sqrt (X*X + Y*Y);
}
vec_t LengthSquared() const
@ -262,14 +263,23 @@ struct TVector2
// Returns the angle (in radians) that the ray (0,0)-(X,Y) faces
double Angle() const
{
return atan2 (X, Y);
return g_atan2 (X, Y);
}
// Returns a rotated vector. TAngle is in radians.
// Returns a rotated vector. angle is in radians.
TVector2 Rotated (double angle)
{
double cosval = cos (angle);
double sinval = sin (angle);
double cosval = g_cos (angle);
double sinval = g_sin (angle);
return TVector2(X*cosval - Y*sinval, Y*cosval + X*sinval);
}
// Returns a rotated vector. angle is in degrees.
template<class T>
TVector2 Rotated(TAngle<T> angle)
{
double cosval = angle.Cos();
double sinval = angle.Sin();
return TVector2(X*cosval - Y*sinval, Y*cosval + X*sinval);
}
@ -490,7 +500,7 @@ struct TVector3
// Vector length
double Length() const
{
return sqrt (X*X + Y*Y + Z*Z);
return g_sqrt (X*X + Y*Y + Z*Z);
}
double LengthSquared() const
@ -573,7 +583,7 @@ struct TMatrix3x3
// (The axis vector must be normalized.)
TMatrix3x3(const Vector3 &axis, double radians)
{
double c = cos(radians), s = sin(radians), t = 1 - c;
double c = g_cos(radians), s = g_sin(radians), t = 1 - c;
/* In comments: A more readable version of the matrix setup.
This was found in Diana Gruber's article "The Mathematics of the
3D Rotation Matrix" at <http://www.makegames.com/3drotation/> and is
@ -782,7 +792,8 @@ struct TAngle
return *this;
}
operator vec_t() const { return Degrees; }
// intentionally disabled so that common math functions cannot be accidentally called with a TAngle.
//operator vec_t() const { return Degrees; }
TAngle operator- () const
{
@ -971,10 +982,26 @@ struct TAngle
return Degrees * (M_PI / 180.0);
}
unsigned BAM() const
unsigned BAMs() const
{
return FLOAT2ANGLE(Degrees);
}
double Cos() const
{
return g_cosdeg(Degrees);
}
double Sin() const
{
return g_sindeg(Degrees);
}
double Tan() const
{
return g_tan(Degrees);
}
};
template<class T>
@ -989,24 +1016,6 @@ inline TAngle<T> ToDegrees (double rad)
return TAngle<T> (T(rad * (180.0 / M_PI)));
}
template<class T>
inline double cos (const TAngle<T> &deg)
{
return cos(ToRadians(deg));
}
template<class T>
inline double sin (const TAngle<T> &deg)
{
return sin(ToRadians(deg));
}
template<class T>
inline double tan (const TAngle<T> &deg)
{
return tan(ToRadians(deg));
}
template<class T>
inline TAngle<T> fabs (const TAngle<T> &deg)
{
@ -1016,13 +1025,13 @@ inline TAngle<T> fabs (const TAngle<T> &deg)
template<class T>
inline TAngle<T> vectoyaw (const TVector2<T> &vec)
{
return (vec_t)atan2(vec.Y, vec.X) * (180.0 / M_PI);
return (vec_t)g_atan2(vec.Y, vec.X) * (180.0 / M_PI);
}
template<class T>
inline TAngle<T> vectoyaw (const TVector3<T> &vec)
{
return (vec_t)atan2(vec.Y, vec.X) * (180.0 / M_PI);
return (vec_t)g_atan2(vec.Y, vec.X) * (180.0 / M_PI);
}
// Much of this is copied from TVector3. Is all that functionality really appropriate?
@ -1192,13 +1201,16 @@ struct TRotator
template<class T>
inline TVector3<T>::TVector3 (const TRotator<T> &rot)
: X(cos(rot.Pitch)*cos(rot.Yaw)), Y(cos(rot.Pitch)*sin(rot.Yaw)), Z(-sin(rot.Pitch))
{
double pcos = rot.Pitch.Cos();
X = pcos * rot.Yaw.Cos();
Y = pcos * rot.Yaw.Sin();
Z = rot.Pitch.Sin();
}
template<class T>
inline TVector2<T>::TVector2(const TRotator<T> &rot)
: X(cos(rot.Yaw)), Y(sin(rot.Yaw))
: X(rot.Yaw.Cos()), Y(rot.Yaw.Sin())
{
}
@ -1206,7 +1218,7 @@ inline TVector2<T>::TVector2(const TRotator<T> &rot)
template<class T>
inline TMatrix3x3<T>::TMatrix3x3(const TVector3<T> &axis, TAngle<T> degrees)
{
double c = cos(degrees), s = sin(degrees), t = 1 - c;
double c = degrees.Cos(), s = degrees.Sin(), t = 1 - c;
double sx = s*axis.X, sy = s*axis.Y, sz = s*axis.Z;
double tx, ty, txx, tyy, u, v;

1
src/zscript/vmexec.cpp
View file

@ -1,6 +1,7 @@
#include <math.h>
#include "vm.h"
#include "xs_Float.h"
#include "math/cmath.h"
#define IMPLEMENT_VMEXEC

40
src/zscript/vmexec.h
View file

@ -1269,7 +1269,7 @@ begin:
OP(LENV):
ASSERTF(a); ASSERTF(B+2);
reg.f[a] = sqrt(reg.f[B] * reg.f[B] + reg.f[B+1] * reg.f[B+1] + reg.f[B+2] * reg.f[B+2]);
reg.f[a] = g_sqrt(reg.f[B] * reg.f[B] + reg.f[B+1] * reg.f[B+1] + reg.f[B+2] * reg.f[B+2]);
NEXTOP;
OP(EQV_R):
@ -1392,30 +1392,30 @@ static double DoFLOP(int flop, double v)
{
case FLOP_ABS: return fabs(v);
case FLOP_NEG: return -v;
case FLOP_EXP: return exp(v);
case FLOP_LOG: return log(v);
case FLOP_LOG10: return log10(v);
case FLOP_SQRT: return sqrt(v);
case FLOP_EXP: return g_exp(v);
case FLOP_LOG: return g_log(v);
case FLOP_LOG10: return g_log10(v);
case FLOP_SQRT: return g_sqrt(v);
case FLOP_CEIL: return ceil(v);
case FLOP_FLOOR: return floor(v);
case FLOP_ACOS: return acos(v);
case FLOP_ASIN: return asin(v);
case FLOP_ATAN: return atan(v);
case FLOP_COS: return cos(v);
case FLOP_SIN: return sin(v);
case FLOP_TAN: return tan(v);
case FLOP_ACOS: return g_acos(v);
case FLOP_ASIN: return g_asin(v);
case FLOP_ATAN: return g_atan(v);
case FLOP_COS: return g_cos(v);
case FLOP_SIN: return g_sin(v);
case FLOP_TAN: return g_tan(v);
case FLOP_ACOS_DEG: return acos(v) * (180 / M_PI);
case FLOP_ASIN_DEG: return asin(v) * (180 / M_PI);
case FLOP_ATAN_DEG: return atan(v) * (180 / M_PI);
case FLOP_COS_DEG: return cos(v * (M_PI / 180));
case FLOP_SIN_DEG: return sin(v * (M_PI / 180));
case FLOP_TAN_DEG: return tan(v * (M_PI / 180));
case FLOP_ACOS_DEG: return g_acos(v) * (180 / M_PI);
case FLOP_ASIN_DEG: return g_asin(v) * (180 / M_PI);
case FLOP_ATAN_DEG: return g_atan(v) * (180 / M_PI);
case FLOP_COS_DEG: return g_cosdeg(v);
case FLOP_SIN_DEG: return g_sindeg(v);
case FLOP_TAN_DEG: return g_tan(v * (M_PI / 180));
case FLOP_COSH: return cosh(v);
case FLOP_SINH: return sinh(v);
case FLOP_TANH: return tanh(v);
case FLOP_COSH: return g_cosh(v);
case FLOP_SINH: return g_sinh(v);
case FLOP_TANH: return g_tanh(v);
}
assert(0);
return 0;