From 161d03231a21fc3c310447395ee40d3ace518308 Mon Sep 17 00:00:00 2001 From: Christoph Oelckers Date: Fri, 11 Mar 2016 15:45:47 +0100 Subject: [PATCH] - added custom math routines for reliability. --- CMakeLists.txt | 2 +- src/CMakeLists.txt | 16 ++ src/b_move.cpp | 3 +- src/fragglescript/t_func.cpp | 23 +- src/g_heretic/a_hereticweaps.cpp | 2 +- src/g_hexen/a_flechette.cpp | 2 +- src/g_shared/a_camera.cpp | 3 +- src/g_shared/a_decals.cpp | 2 +- src/g_shared/a_movingcamera.cpp | 8 +- src/g_shared/a_quake.cpp | 2 +- src/math/asin.c | 315 +++++++++++++++++++++ src/math/atan.c | 393 +++++++++++++++++++++++++++ src/math/cmath.h | 139 ++++++++++ src/math/const.c | 252 +++++++++++++++++ src/math/cosh.c | 83 ++++++ src/math/exp.c | 182 +++++++++++++ src/math/fastsin.cpp | 102 +++++++ src/math/isnan.c | 237 ++++++++++++++++ src/math/log.c | 341 +++++++++++++++++++++++ src/math/log10.c | 250 +++++++++++++++++ src/math/mconf.h | 199 ++++++++++++++ src/math/mtherr.c | 102 +++++++ src/math/polevl.c | 97 +++++++ src/math/readme.txt | 15 + src/math/sin.c | 387 ++++++++++++++++++++++++++ src/math/sinh.c | 148 ++++++++++ src/math/sqrt.c | 178 ++++++++++++ src/math/tan.c | 304 +++++++++++++++++++++ src/math/tanh.c | 141 ++++++++++ src/nodebuild_events.cpp | 2 +- src/nodebuild_utility.cpp | 3 +- src/p_acs.cpp | 4 +- src/p_buildmap.cpp | 2 +- src/p_effect.cpp | 2 +- src/p_enemy.cpp | 7 +- src/p_glnodes.cpp | 2 +- src/p_map.cpp | 7 +- src/p_maputl.cpp | 2 +- src/p_mobj.cpp | 2 +- src/p_setup.cpp | 4 +- src/p_things.cpp | 9 +- src/po_man.cpp | 2 +- src/portal.cpp | 9 +- src/r_utility.cpp | 11 +- src/r_utility.h | 2 +- src/thingdef/thingdef_codeptr.cpp | 3 +- src/thingdef/thingdef_expression.cpp | 27 +- src/vectors.h | 76 +++--- src/zscript/vmexec.cpp | 1 + src/zscript/vmexec.h | 40 +-- 50 files changed, 4025 insertions(+), 120 deletions(-) create mode 100644 src/math/asin.c create mode 100644 src/math/atan.c create mode 100644 src/math/cmath.h create mode 100644 src/math/const.c create mode 100644 src/math/cosh.c create mode 100644 src/math/exp.c create mode 100644 src/math/fastsin.cpp create mode 100644 src/math/isnan.c create mode 100644 src/math/log.c create mode 100644 src/math/log10.c create mode 100644 src/math/mconf.h create mode 100644 src/math/mtherr.c create mode 100644 src/math/polevl.c create mode 100644 src/math/readme.txt create mode 100644 src/math/sin.c create mode 100644 src/math/sinh.c create mode 100644 src/math/sqrt.c create mode 100644 src/math/tan.c create mode 100644 src/math/tanh.c diff --git a/CMakeLists.txt b/CMakeLists.txt index 2e0bf8b77..bbe4e92ba 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 94ac6fb32..609706837 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 diff --git a/src/b_move.cpp b/src/b_move.cpp index 6138ba3ad..3a4cf3c3a 100644 --- a/src/b_move.cpp +++ b/src/b_move.cpp @@ -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); } diff --git a/src/fragglescript/t_func.cpp b/src/fragglescript/t_func.cpp index e034ab7ad..13598e269 100644 --- a/src/fragglescript/t_func.cpp +++ b/src/fragglescript/t_func.cpp @@ -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"); diff --git a/src/g_heretic/a_hereticweaps.cpp b/src/g_heretic/a_hereticweaps.cpp index d85faaeb5..7a7d21f76 100644 --- a/src/g_heretic/a_hereticweaps.cpp +++ b/src/g_heretic/a_hereticweaps.cpp @@ -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); diff --git a/src/g_hexen/a_flechette.cpp b/src/g_hexen/a_flechette.cpp index 56ad33342..ef7d092f3 100644 --- a/src/g_hexen/a_flechette.cpp +++ b/src/g_hexen/a_flechette.cpp @@ -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]); diff --git a/src/g_shared/a_camera.cpp b/src/g_shared/a_camera.cpp index 7febf80d1..fee44a3eb 100644 --- a/src/g_shared/a_camera.cpp +++ b/src/g_shared/a_camera.cpp @@ -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) { diff --git a/src/g_shared/a_decals.cpp b/src/g_shared/a_decals.cpp index d37250887..ab0d50933 100644 --- a/src/g_shared/a_decals.cpp +++ b/src/g_shared/a_decals.cpp @@ -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) diff --git a/src/g_shared/a_movingcamera.cpp b/src/g_shared/a_movingcamera.cpp index b67e7f8e3..e186f8592 100644 --- a/src/g_shared/a_movingcamera.cpp +++ b/src/g_shared/a_movingcamera.cpp @@ -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); } diff --git a/src/g_shared/a_quake.cpp b/src/g_shared/a_quake.cpp index 205d43ea2..42d2c9df5 100644 --- a/src/g_shared/a_quake.cpp +++ b/src/g_shared/a_quake.cpp @@ -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))); } //========================================================================== diff --git a/src/math/asin.c b/src/math/asin.c new file mode 100644 index 000000000..50419597f --- /dev/null +++ b/src/math/asin.c @@ -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; +} diff --git a/src/math/atan.c b/src/math/atan.c new file mode 100644 index 000000000..3c5d4394c --- /dev/null +++ b/src/math/atan.c @@ -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 ); +} diff --git a/src/math/cmath.h b/src/math/cmath.h new file mode 100644 index 000000000..b161e2106 --- /dev/null +++ b/src/math/cmath.h @@ -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 \ No newline at end of file diff --git a/src/math/const.c b/src/math/const.c new file mode 100644 index 000000000..2d9ff7b33 --- /dev/null +++ b/src/math/const.c @@ -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 diff --git a/src/math/cosh.c b/src/math/cosh.c new file mode 100644 index 000000000..d69f5b344 --- /dev/null +++ b/src/math/cosh.c @@ -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 ); +} diff --git a/src/math/exp.c b/src/math/exp.c new file mode 100644 index 000000000..7d2b1e88f --- /dev/null +++ b/src/math/exp.c @@ -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); +} diff --git a/src/math/fastsin.cpp b/src/math/fastsin.cpp new file mode 100644 index 000000000..1c25955b0 --- /dev/null +++ b/src/math/fastsin.cpp @@ -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 +#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); + } +} + diff --git a/src/math/isnan.c b/src/math/isnan.c new file mode 100644 index 000000000..b5341e650 --- /dev/null +++ b/src/math/isnan.c @@ -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 +} diff --git a/src/math/log.c b/src/math/log.c new file mode 100644 index 000000000..34a8faf02 --- /dev/null +++ b/src/math/log.c @@ -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 ); +} diff --git a/src/math/log10.c b/src/math/log10.c new file mode 100644 index 000000000..df67d30b9 --- /dev/null +++ b/src/math/log10.c @@ -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 ); +} diff --git a/src/math/mconf.h b/src/math/mconf.h new file mode 100644 index 000000000..6e8e982ff --- /dev/null +++ b/src/math/mconf.h @@ -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 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; diff --git a/src/math/mtherr.c b/src/math/mtherr.c new file mode 100644 index 000000000..0d94cd14d --- /dev/null +++ b/src/math/mtherr.c @@ -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 +#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 ); +} diff --git a/src/math/polevl.c b/src/math/polevl.c new file mode 100644 index 000000000..4d050fbfc --- /dev/null +++ b/src/math/polevl.c @@ -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 ); +} diff --git a/src/math/readme.txt b/src/math/readme.txt new file mode 100644 index 000000000..965335ee6 --- /dev/null +++ b/src/math/readme.txt @@ -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 + \ No newline at end of file diff --git a/src/math/sin.c b/src/math/sin.c new file mode 100644 index 000000000..3eba0e22f --- /dev/null +++ b/src/math/sin.c @@ -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 ); +} diff --git a/src/math/sinh.c b/src/math/sinh.c new file mode 100644 index 000000000..b1d2a7500 --- /dev/null +++ b/src/math/sinh.c @@ -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)) ); +} diff --git a/src/math/sqrt.c b/src/math/sqrt.c new file mode 100644 index 000000000..22ebdc669 --- /dev/null +++ b/src/math/sqrt.c @@ -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); +} diff --git a/src/math/tan.c b/src/math/tan.c new file mode 100644 index 000000000..8ef8af55c --- /dev/null +++ b/src/math/tan.c @@ -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 ); +} diff --git a/src/math/tanh.c b/src/math/tanh.c new file mode 100644 index 000000000..7590eca5c --- /dev/null +++ b/src/math/tanh.c @@ -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 ); +} diff --git a/src/nodebuild_events.cpp b/src/nodebuild_events.cpp index 345a6b6be..84e187081 100644 --- a/src/nodebuild_events.cpp +++ b/src/nodebuild_events.cpp @@ -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); } diff --git a/src/nodebuild_utility.cpp b/src/nodebuild_utility.cpp index a2edd19b9..d3d560f4c 100644 --- a/src/nodebuild_utility.cpp +++ b/src/nodebuild_utility.cpp @@ -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; diff --git a/src/p_acs.cpp b/src/p_acs.cpp index 3aefcb46a..bf62a1f0a 100644 --- a/src/p_acs.cpp +++ b/src/p_acs.cpp @@ -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()); diff --git a/src/p_buildmap.cpp b/src/p_buildmap.cpp index 3da00fbaf..a866f852f 100644 --- a/src/p_buildmap.cpp +++ b/src/p_buildmap.cpp @@ -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; diff --git a/src/p_effect.cpp b/src/p_effect.cpp index 63355144d..c745ce546 100644 --- a/src/p_effect.cpp +++ b/src/p_effect.cpp @@ -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); diff --git a/src/p_enemy.cpp b/src/p_enemy.cpp index 2c1339921..e8dd62708 100644 --- a/src/p_enemy.cpp +++ b/src/p_enemy.cpp @@ -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 diff --git a/src/p_glnodes.cpp b/src/p_glnodes.cpp index ed4b1d345..84951a3a3 100644 --- a/src/p_glnodes.cpp +++ b/src/p_glnodes.cpp @@ -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) { diff --git a/src/p_map.cpp b/src/p_map.cpp index f021c3cc6..42a19bc3f 100644 --- a/src/p_map.cpp +++ b/src/p_map.cpp @@ -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 diff --git a/src/p_maputl.cpp b/src/p_maputl.cpp index 526626b22..31038d9db 100644 --- a/src/p_maputl.cpp +++ b/src/p_maputl.cpp @@ -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); diff --git a/src/p_mobj.cpp b/src/p_mobj.cpp index 73214763e..3dfb90693 100644 --- a/src/p_mobj.cpp +++ b/src/p_mobj.cpp @@ -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))); } } diff --git a/src/p_setup.cpp b/src/p_setup.cpp index 5232ab37f..4470d6fa1 100644 --- a/src/p_setup.cpp +++ b/src/p_setup.cpp @@ -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)))< 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) { diff --git a/src/p_things.cpp b/src/p_things.cpp index eecf57ebc..3ab4811fc 100644 --- a/src/p_things.cpp +++ b/src/p_things.cpp @@ -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)) diff --git a/src/po_man.cpp b/src/po_man.cpp index 64f1fde61..9e597719e 100644 --- a/src/po_man.cpp +++ b/src/po_man.cpp @@ -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 diff --git a/src/portal.cpp b/src/portal.cpp index 63dfeafe6..518ae15ee 100644 --- a/src/portal.cpp +++ b/src/portal.cpp @@ -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); } diff --git a/src/r_utility.cpp b/src/r_utility.cpp index f4bd85a9a..336b4cc2c 100644 --- a/src/r_utility.cpp +++ b/src/r_utility.cpp @@ -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++) { diff --git a/src/r_utility.h b/src/r_utility.h index 5c953c7ab..d12d8f5ae 100644 --- a/src/r_utility.h +++ b/src/r_utility.h @@ -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. diff --git a/src/thingdef/thingdef_codeptr.cpp b/src/thingdef/thingdef_codeptr.cpp index b9ad61e6f..761273666 100644 --- a/src/thingdef/thingdef_codeptr.cpp +++ b/src/thingdef/thingdef_codeptr.cpp @@ -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) diff --git a/src/thingdef/thingdef_expression.cpp b/src/thingdef/thingdef_expression.cpp index 7ecc3c3c4..441e53371 100644 --- a/src/thingdef/thingdef_expression.cpp +++ b/src/thingdef/thingdef_expression.cpp @@ -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) diff --git a/src/vectors.h b/src/vectors.h index ecd37e5a7..0fb8c0905 100644 --- a/src/vectors.h +++ b/src/vectors.h @@ -43,6 +43,7 @@ #include #include #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 + TVector2 Rotated(TAngle 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 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 @@ -989,24 +1016,6 @@ inline TAngle ToDegrees (double rad) return TAngle (T(rad * (180.0 / M_PI))); } -template -inline double cos (const TAngle °) -{ - return cos(ToRadians(deg)); -} - -template -inline double sin (const TAngle °) -{ - return sin(ToRadians(deg)); -} - -template -inline double tan (const TAngle °) -{ - return tan(ToRadians(deg)); -} - template inline TAngle fabs (const TAngle °) { @@ -1016,13 +1025,13 @@ inline TAngle fabs (const TAngle °) template inline TAngle vectoyaw (const TVector2 &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 inline TAngle vectoyaw (const TVector3 &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 inline TVector3::TVector3 (const TRotator &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 inline TVector2::TVector2(const TRotator &rot) - : X(cos(rot.Yaw)), Y(sin(rot.Yaw)) + : X(rot.Yaw.Cos()), Y(rot.Yaw.Sin()) { } @@ -1206,7 +1218,7 @@ inline TVector2::TVector2(const TRotator &rot) template inline TMatrix3x3::TMatrix3x3(const TVector3 &axis, TAngle 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; diff --git a/src/zscript/vmexec.cpp b/src/zscript/vmexec.cpp index ebce8c5f2..070aeb30a 100644 --- a/src/zscript/vmexec.cpp +++ b/src/zscript/vmexec.cpp @@ -1,6 +1,7 @@ #include #include "vm.h" #include "xs_Float.h" +#include "math/cmath.h" #define IMPLEMENT_VMEXEC diff --git a/src/zscript/vmexec.h b/src/zscript/vmexec.h index eb54e7b76..954a8cde8 100644 --- a/src/zscript/vmexec.h +++ b/src/zscript/vmexec.h @@ -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;