133 lines
2.8 KiB
C
Executable File
133 lines
2.8 KiB
C
Executable File
/* HPmath.c
|
|
*
|
|
* Copyright (C) 1994 by TV Games, Inc.
|
|
* All Rights Reserved
|
|
*
|
|
* These routines are for use only when a high degree of precision
|
|
* are nessesary.
|
|
*/
|
|
|
|
|
|
|
|
#define BITS 23 /* There are 23 bits in the mantissa */
|
|
#define MAXX 88.72283906 /* ln(HUGE_VAL) */
|
|
#define MAXH 89.41598624 /* ln(HUGE_VAL) + ln(2) */
|
|
#define TWO23 8388608 /* 2 ^ BITS */
|
|
#define XBIG 8.664339757 /* (BITS/2 + 1) * ln(2) */
|
|
|
|
|
|
/* macros used in sin and cos */
|
|
|
|
#define INVSPI 0.31830988618379067154
|
|
#define HALFPI 1.57079632679489661923
|
|
|
|
#define C1 3.140625
|
|
#define C2 9.67653589793e-4
|
|
|
|
#define R1 -0.1666665668e+0
|
|
#define R2 0.8333025139e-2
|
|
#define R3 -0.1980741872e-3
|
|
#define R4 0.2601903036e-5
|
|
|
|
|
|
|
|
/* HPsin() - High Precision sine
|
|
*
|
|
* Based on the algorithm from "Software Manual for the Elementary
|
|
* Functions", Cody and Waite, Prentice Hall 1980, chapter 8.
|
|
*
|
|
* N = round(x / PI)
|
|
* f = x - N * PI
|
|
* g = f * f
|
|
* R = polynomial expansion
|
|
*
|
|
* result = f + f * R
|
|
*
|
|
* if x < 0, result = - result
|
|
* if N is even, result = - result
|
|
*
|
|
* This will return the wrong result for x >= MAXINT * PI
|
|
*/
|
|
double HPsin(double x)
|
|
{
|
|
double d,y,xn,f,g,rg;
|
|
float sgn = (x < 0) ? -1.0 : 1.0;
|
|
int n;
|
|
|
|
x = fabs(x);
|
|
n = (int) ((x * INVSPI) + 0.5);
|
|
xn = (double) n;
|
|
|
|
/*
|
|
* if n is odd, negate the sign
|
|
*/
|
|
if (n % 2) sgn = -sgn;
|
|
|
|
/*
|
|
* f = x - xn * PI (but mathematically more stable)
|
|
*/
|
|
f = (x - xn * C1) - xn * C2;
|
|
|
|
/*
|
|
* determine polynomial expression
|
|
*/
|
|
g = f * f;
|
|
|
|
rg = (((R4 * g + R3) * g + R2) * g + R1) * g;
|
|
|
|
return (sgn * (f + f * rg));
|
|
}
|
|
|
|
|
|
|
|
/* HPcos() - High Precision Cosine
|
|
*
|
|
* Based on the algorithm from "Software Manual for the Elementary
|
|
* Functions", Cody and Waite, Prentice Hall 1980, chapter 8.
|
|
*
|
|
* N = round(x / PI + 1/2) - 0.5
|
|
* f = x - N * PI
|
|
* g = f * f
|
|
* R = polynomial expression
|
|
*
|
|
* result = f + f * R
|
|
* if N is even, result = - result
|
|
*
|
|
* This will return the wrong result for x >= MAXINT * PI
|
|
*/
|
|
double HPcos(double x)
|
|
{
|
|
float sgn; /* the sign of the result */
|
|
double xn,f,g,rg;
|
|
int n;
|
|
|
|
/*
|
|
* cos(x) = cos(-x)
|
|
*/
|
|
x = fabs(x);
|
|
|
|
/*
|
|
* n = round(x/PI + 1/2) (can be rounded this way, since positive number)
|
|
*/
|
|
n = (int) (((x + HALFPI) * INVSPI) + 0.5);
|
|
xn = (double) n - 0.5;
|
|
|
|
/*
|
|
* if n is odd, negate the sign
|
|
*/
|
|
sgn = (n % 2) ? -1.0 : 1.0;
|
|
|
|
/*
|
|
* f = x - xn * PI (but more mathematically stable)
|
|
*/
|
|
f = (x - xn * C1) - xn * C2;
|
|
|
|
/*
|
|
* determine polynomial expression
|
|
*/
|
|
g = f * f;
|
|
|
|
rg = (((R4 * g + R3) * g + R2) * g + R1) * g;
|
|
return (sgn * (f + f * rg));
|
|
}
|