home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Simtel MSDOS - Coast to Coast
/
simteldosarchivecoasttocoast2.iso
/
c
/
math.c
< prev
next >
Wrap
C/C++ Source or Header
|
1994-03-04
|
12KB
|
520 lines
/*****************************************************************************
* *
* MATH.C Mathematical subroutines for use with Lattice or Microsoft- *
* compiler. *
* By Max R. D^Arsteler *
* *
*****************************************************************************
*/
/*****************************************************************************
* This file contains the combined source files of the most
* important mathe- matical formulas written for the Lattice C Compiler.
* It makes no use of the 8087. Some of the functions make use of the
* internal binary representation of double precision numbers. As the
* internal data representation may vary between different compilers,
* watch out for errors e.g. in the power function. The function values
* are calculated through polynominal or rational approximations. The
* results are for allmost all functions precise for the first 9 to 10
* decimal digits. For still higher precision or checking of values use
* functions from file mathf.c (not available here), which contains
* Taylor's expansions series. As this requires many iterations with
* multiplications and divisions, the calculation time is too long for use
* in regular programs. In order to use these functions, a program has to
* include the header file "math.h". Furthermore, either the compiled
* version (math.obj) or a library with the compiled functions has to be
* linked to the program. As exemples see included testprograms
* "tstsin.c" and "tstpow.c" To compile: "mcb tstsin <CR>", to link:"link
* c+%1+math.obj, %1, ,mc <CR>". I recommend to compile the source of
* each function seperately with your personal compiler and use the MS-DOS
* librarian to build a library of just the functions you need. Some of
* the more esotheric functions used on "UNIX" are such as Bessel
* functions or the gamma function are not implemented. Evidently, this
* is not a professional implementation.
*
* Rockville, 11/27/83 Max R. D^Arsteler
* 12405 Village Sq. Terrace
* Rockville Md. 20852 H 301 984-1168
******************************************************************************/
/* SQRT.C
* Calculates square root
* Precision: 11.05 E-10
* Range: 0 to e308
* Author: Max R. D^Arsteler
*/
double sqrt(x)
double x;
{
typedef union { /* Makes use of internal data */
double d; /* representation */
unsigned u[4];
} DBL;
double y, re, p, q;
unsigned ex; /* exponens*/
DBL *xp;
xp = (DBL *)&x;
ex = xp->u[3] & ~0100017; /* save exponens */
re = ex & 020 ? 1.4142135623730950488 : 1.0;
ex = ex - 037740 >> 1;
ex &= ~0100017;
xp->u[3] &= ~0177760; /* erase exponens and sign*/
xp->u[3] |= 037740; /* arrange for mantissa in range 0.5 to 1.0 */
if (xp->d < .7071067812) { /* multiply by sqrt(2) if mantissa < 0.7 */
xp->d *= 1.4142135623730950488;
if (re > 1.0) re = 1.189207115002721;
else re = 1.0 / 1.18920711500271;
}
p = .54525387389085 + /* Polynomial approximation */
xp->d * (13.65944682358639 + /* from: COMPUTER APPROXIMATIONS*/
xp->d * (27.090122712655 + /* Hart, J.F. et al. 1968 */
xp->d * 6.43107724139784));
q = 4.17003771413707 +
xp->d * (24.84175210765124 +
xp->d * (17.71411083016702 +
xp->d ));
xp->d = p / q;
xp->u[3] = xp->u[3] + ex & ~0100000;
xp->d *= re;
return(xp->d);
}
/* EXP2.C
* Calculates exponens to base 2 using polynomial approximations
* Precision: 10E-10
* Range:
* Author: Max R. D^Arsteler
*/
double exp2(x)
double x;
{
typedef union {
double d;
unsigned u[4];
} DBL;
double y, x2, p, q, re;
int ix;
DBL *xp, *yp;
xp = (DBL *)&x;
y = 0.0;
yp = (DBL *)&y;
if(xp->d > 1023.0 || xp->d < -1023.0) return (1E307);
ix = (int) xp->d;
yp->u[3] += ix + 1023 << 4;
yp->u[3] &= ~0100017;
if ((xp->d -= (double) ix) == 0.0) return(yp->d);
if (xp->d < 0.0) {
yp->u[3] -= 1 << 4;
yp->u[3] &= ~0100017;
xp->d++;
}
if (xp->d >= 0.5) { /* adjust to range 0-0.5 */
xp->d -= 0.5;
re = 1.41421356237309504880;
}
else re = 1.0;
x2 = xp->d * xp->d;
p = xp->d * (7.2152891511493 +
x2 * 0.0576900723731);
q = 20.8189237930062 +
x2 ;
xp->d = (q + p) / (q - p);
yp->d *= re * xp->d;
return(yp->d);
}
/* POW.C
* Calculates y^x
* Precision:
* Range: o to big for y, +/-1023 for x
* Author: Max R. D^Arsteler, 10/2/83
*/
double pow(y,x)
double y, x;
{
typedef union {
double d;
unsigned u[4];
} DBL;
double z, w, p, p2, q, re;
unsigned ex; /* exponens*/
int iz;
DBL *yp, *zp, *wp;
yp = (DBL *)&y;
if (yp->d <= 0.0) y = -y;
z = 0.0;
zp = (DBL *)&z;
zp->u[3] = yp->u[3] & ~0100017; /* save exponens */
iz = (zp->u[3] >> 4)-1023;
if ((yp->d - zp->d) == 0.0)
z = (double)iz;
else {
yp->u[3] -= ++iz << 4; /* arrange for range 0.5 to 0.99999999999 */
yp->d *= 1.4142135623730950488; /* shift for 1/sqrt(2) to sqrt(2) */
p = (yp->d - 1.0) / (yp->d + 1.0);
p2 = p * p;
z = p * (2.000000000046727 + /* Polynomial approximation */
p2 * (0.666666635059382 + /* from: COMPUTER APPROXIMATIONS*/
p2 * (0.4000059794795 + /* Hart, J.F. et al. 1968 */
p2 * (0.28525381498 +
p2 * 0.2376245609 ))));
z = z * 1.442695040888634 + (double)iz - 0.5;
}
z *= x;
w = 0.0;
wp = (DBL *)&w;
if(zp->d > 1023.0 || zp->d < -1023.0) return (1E307);
iz = (int) zp->d;
wp->u[3] += iz + 1023 << 4;
wp->u[3] &= ~0100017;
if ((zp->d -= (double) iz) == 0.0) return(wp->d);
while (zp->d < 0.0) {
wp->u[3] -= 1 << 4;
wp->u[3] &= ~0100017;
zp->d++;
}
if (zp->d >= 0.5) { /* adjust to range 0-0.5 */
zp->d -= 0.5;
re = 1.41421356237309504880;
}
else re = 1.0;
p2 = zp->d * zp->d;
p = zp->d * (7.2152891511493 +
p2 * 0.0576900723731);
q = 20.8189237930062 +
p2 ;
zp->d = re * wp->d * (q + p) / (q - p);
return(zp->d);
}
/* FMOD.C
* Returns the number f such that x = i*y + f. i is an integer, and
* 0 <= f < y.
*/
double fmod(x, y)
double x, y;
{
double zint, z;
int i;
z = x / y;
zint = 0.0;
while (z > 32768.0) {
zint += 32768.0;
z -= 32768;
}
while (z < -32768.0) {
zint -= 32768.0;
z += 32768.0;
}
i = (int) z;
zint += (double) i;
return( x - zint * y);
}
/*ASIN.C
*Calculates arcsin(x)
*Range: 0 <= x <= 1
*Precision: +/- .000,000,02
*Header: math.h
*Author: Max R. D^Arsteler
*/
extern double sqrt();
double asin(x)
double x;
{
double y;
int sign;
if (x > 1.0 || x < -1.0) exit(1);
sign = 0;
if (x < 0) {
sign = 1;
x = -x;
}
y = ((((((-.0012624911 * x
+ .0066700901) * x
- .0170881256) * x
+ .0308918810) * x
- .0501743046) * x
+ .0889789874) * x
- .2145988016) * x
+1.5707963050;
y = 1.57079632679 - sqrt(1.0 - x) * y;
if (sign) y = -y;
return(y);
}
/* SIN.C
* Calculates sin(x), angle x must be in rad.
* Range: -pi/2 <= x <= pi/2
* Precision: +/- .000,000,005
* Header: math.h
* Author: Max R. D^Arsteler
*/
double sin(x)
double x;
{
double xi, y, q, q2;
int sign;
xi = x; sign = 1;
while (xi < -1.57079632679489661923) xi += 6.28318530717958647692;
while (xi > 4.71238898038468985769) xi -= 6.28318530717958647692;
if (xi > 1.57079632679489661923) {
xi -= 3.141592653589793238462643;
sign = -1;
}
q = xi / 1.57079632679; q2 = q * q;
y = ((((.00015148419 * q2
- .00467376557) * q2
+ .07968967928) * q2
- .64596371106) * q2
+1.57079631847) * q;
return(sign < 0? -y : y);
}
/* LOG.C
* Calculates natural logarithmus
* Precision: 11.56 E-10
* Range: 0 to e308
* Author: Max R. D^Arsteler
*/
double log(x)
double x;
{
typedef union {
double d;
unsigned u[4];
} DBL;
double y, z, z2, p;
unsigned ex; /* exponens*/
int ix;
DBL *xp, *yp;
xp = (DBL *)&x;
if (xp->d <= 0.0) return(y);
y = 0.0;
yp = (DBL *)&y;
yp->u[3] = xp->u[3] & ~0100017; /* save exponens */
ix = (yp->u[3] >> 4)-1023;
if ((xp->d - yp->d) == 0.0) return( .693147180559945 * (double)ix);
xp->u[3] -= ++ix << 4; /* arrange for range 0.5 to 0.99999999999 */
xp->d *= 1.4142135623730950488; /* shift for 1/sqrt(2) to sqrt(2) */
z = (xp->d - 1.0) / (xp->d + 1.0);
z2 = z * z;
y = z * (2.000000000046727 + /* Polynomial approximation */
z2 * (0.666666635059382 + /* from: COMPUTER APPROXIMATIONS*/
z2 * (0.4000059794795 + /* Hart, J.F. et al. 1968 */
z2 * (0.28525381498 +
z2 * 0.2376245609 ))));
y = y + .693147180559945 * ((double)ix - 0.5);
return(yp->d);
}
/* ATAN.C
* Calculates arctan(x)
* Range: -infinite <= x <= infinite (Output -pi/2 to +pi/2)
* Precision: +/- .000,000,04
* Header: math.h
* Author: Max R. D^Arsteler 9/15/83
*/
double atan(x)
double x;
{
double xi, q, q2, y;
int sign;
xi = (x < 0. ? -x : x);
q = (xi - 1.0) / (xi + 1.0); q2 = q * q;
y = ((((((( - .0040540580 * q2
+ .0218612286) * q2
- .0559098861) * q2
+ .0964200441) * q2
- .1390853351) * q2
+ .1994653599) * q2
- .3332985605) * q2
+ .9999993329) * q + 0.785398163397;
return(x < 0. ? -y: y);
}
/* FLOOR.C
* Returns largest integer not greater than x
* Author: Max R. D^Arsteler 9/26/83
*/
double floor(x)
double x;
{
double y;
int ix;
y = 0.0;
while (x >= 32768.0) {
y += 32768.0;
x -= 32768.0;
}
while (x <= -32768.0) {
y -= 32768.0;
x += 32768.0;
}
if (x > 0.0) ix = (int) x;
else ix = (int)(x - 0.9999999999);
return( y + (double) ix);
}
/* CEIL.C
* Returns smallest integer not less than x
* Author: Max R. D^Arsteler 9/26/83
*/
double ceil(x)
double x;
{
double y;
int ix;
y = 0.0;
while (x >= 32768.0) {
y += 32768.0;
x -= 32768.0;
}
while (x <= -32768.0) {
y -= 32768.0;
x += 32768.0;
}
if (x > 0.0) ix = (int) (x + 0.999999999999999);
else ix = (int) x;
return( y + (double) ix);
}
/* EXP.C
* Calculates exponens of x to base e
* Range: +/- exp(88)
* Precision: +/- .000,000,000,1
* Author: Max R. D^Arsteler, 9/20/83
*/
double exp(xi)
double xi;
{
double y, ex, px, nn, ds, in;
if (xi > 88.0) return(1.7014117331926443e38);
if (xi < -88.0) return(0.0);
ex = 1.0;
while( xi > 1.0) {
ex *= 2.718281828459; /* const. e */
xi--;
}
while( xi < -1.0) {
ex *= .367879441171; /* 1/e */
xi++;
}
/* Slow, but more precise Taylor expansion series */
y = ds = 1.0; nn = 0.0;
while ((ds < 0.0 ? -ds : ds) > 0.000000000001) {
px = xi/++nn; /* above precision required */
ds *= px;
y += ds;
}
y *= ex;
/* Chebyshev polynomials: fast, but less precise then expected!
* xi = -xi;
* y = (((((.0000006906 * xi
* +.0000054302) * xi
* +.0001715620) * xi
* +.0025913712) * xi
* +.0312575832) * xi
* +.2499986842) * xi
* +1.0;
* y = ex / (y * y * y * y);
*/
return(y);
}
/* LOG10.C
* Approximation for logarithm of basis 10
* Range 0 < x < 1e+38
* Precision: +/- 0.000,000,1
* Header: math.h
* Author: Max R. D^Arsteler 9/15/83
*/
/* Method of Chebyshev polynomials */
/* C. Hastings, jr. 1955 */
double log10(x)
double x;
{
double xi, y, q, q2;
int ex;
if (x <= 0.0) return(0.); /* Error!! */
ex = 0.0; xi = x;
while (xi < 1.0 ) {
xi *= 10.0;
ex--;
}
while (xi > 10.0) {
xi *= 0.1;
ex++;
}
q = (xi - 3.16227766) / (xi + 3.16227766); q2 = q * q;
y = ((((.191337714 * q2
+ .094376476) * q2
+ .177522071) * q2
+ .289335524) * q2
+ .868591718) * q + .5;
y += (double) ex;
return(y);
}
/* PW10.C
* Calculates 10 power x
* Range: 0 <= x <= 1
* Precision: +/- 0.000,000,005
* Header: math.lib
* Author: Max R. D^Arsteler 9/15/83
*/
double pw10(x)
double x;
{
double xi,ex,y;
if (x > 38.0) return(1.7014117331926443e38);
if (x < -38.0) return(0.0);
xi = x; ex = 1.0;
while (xi > 1.0) {
ex *= 10.0;
xi--;
}
while (xi < 0.0) {
ex *= 0.1;
xi++;
}
y = ((((((.00093264267 * xi
+ .00255491796) * xi
+ .01742111988) * xi
+ .07295173666) * xi
+ .25439357484) * xi
+ .66273088429) * xi
+1.15129277603) * xi + 1.0;
y *= y;
return(y*ex);
}