/*
* NAME: math.c
* DESCRIPTION: mathematics library
*/
# include <math.h>
# include <dgd/float.h>
/*
* NAME: abs()
* DESCRIPTION: return integer absolute value
*/
static
int abs(int num)
{ return (num < 0 ? -num : num); }
/*
* NAME: sqrt()
* DESCRIPTION: return floating point square root
*/
static
float sqrt(float x)
{
mixed *split;
float f, y, z;
int n;
if (x == 0.0)
return 0.0;
if (x < 0.0)
error("Square root of negative number");
split = frexp(x);
f = split[0];
n = split[1];
y = 0.4173102246 + 0.5901604053 * f;
z = y + f / y;
y = 0.25 * z + f / z;
z = 0.5 * (y + f / y);
if (n & 1)
{
++n;
z *= SQRT05;
}
return ldexp(z, n / 2);
}
/*
* NAME: sin_cos()
* DESCRIPTION: compute sine or cosine
*/
private
float sin_cos(float x, float y, float sgn)
{
float f, g, r, xn, *tmp;
tmp = modf(ONEOPI * y + 0.5);
xn = tmp[1];
tmp = modf(0.5 * xn);
if ((int) (2.0 * tmp[0]) & 1 == 1)
sgn = -sgn;
r = (x < 0.0 ? -x : x);
if (r != y)
xn -= 0.5;
f = sgn * ((r - xn * 3.1416015625) - xn * -8.90891020676153736e-6);
if (f < ROOTEPS && f > -ROOTEPS)
return f;
g = f * f;
r = (((((((0.272047909578888462e-14 * g +
-0.764291780689104677e-12) * g +
0.160589364903715891e-9) * g +
-0.250521067982745845e-7) * g +
0.275573192101527561e-5) * g +
-0.198412698412018405e-3) * g +
0.833333333333316503e-2) * g +
-0.166666666666666651e0) * g;
return f + f * r;
}
/*
* NAME: sin()
* DESCRIPTION: return sine of argument
*/
static
float sin(float x)
{
if (x < 0.0)
return sin_cos(x, -x, -1.0);
else
return sin_cos(x, x, 1.0);
}
/*
* NAME: cos()
* DESCRIPTION: return cosine of argument
*/
static
float cos(float x)
{
return sin_cos(x, fabs(x) + PIO2, 1.0);
}
/*
* NAME: tan_cot()
* DESCRIPTION: compute tangent or cotangent
*/
private
float tan_cot(float x, int flag)
{
float f, g, p, q, *tmp;
int n;
f = x < 0.0 ? -x : x;
if (x < 0.0)
tmp = modf(x * ONEOPIO2 - 0.5);
else
tmp = modf(x * ONEOPIO2 + 0.5);
g = tmp[1];
tmp = modf(0.5 * g);
p = tmp[1];
n = (int) (2.0 * tmp[0]);
f = (x - g * 1.57080078125) - g * -4.45445510338076868e-6;
if (f < ROOTEPS && f > -ROOTEPS)
{
p = f;
q = 1.0;
}
else
{
g = f * f;
p = ((-0.178617073422544267e-4 * g +
0.342488782358905900e-2) * g +
-0.133383500064219607e0) * g * f + f;
q = (((0.498194339937865123e-6 * g +
-0.311815319070100273e-3) * g +
0.256638322894401129e-1) * g +
-0.466716833397552942e0) * g + 1.0;
}
if (flag == 0)
{
if (n == 0)
return p / q;
else
return -q / p;
}
else
{
if (n == 0)
return q / p;
else
return -p / q;
}
}
/*
* NAME: tan()
* DESCRIPTION: compute tangent of argument
*/
static
float tan(float x)
{
return tan_cot(x, 0);
}
/*
* NAME: cotan()
* DESCRIPTION: compute cotangent of argument
*/
static
float cotan(float x)
{
if (x < (1.0 / FLT_MAX) && x > -(1.0 / FLT_MAX))
return x < 0.0 ? -FLT_MAX : FLT_MAX;
return tan_cot(x, 1);
}
/*
* NAME: copysign()
* DESCRIPTION: return first argument with sign of second
*/
static
float copysign(float x, float y)
{
return ((x < 0.0 && y < 0.0) || (x > 0.0 && y > 0.0)) ? x : -x;
}
# define small 1.0e-9
# define big 1.0e+18
/*
* NAME: atan2()
* DESCRIPTION: return arctangent of y/x
*/
static
float atan2(float y, float x)
{
float t, z, signy, signx, hi, lo;
int k, m;
signy = copysign(1.0, y);
signx = copysign(1.0, x);
if (x == 1.0)
t = y = copysign(y, 1.0);
else
{
if (y == 0.0)
return (signx == 1.0) ? y : copysign(PI, signy);
if (x == 0.0)
return copysign(PIO2, signy);
x = copysign(x, 1.0);
y = copysign(y, 1.0);
if ((m = (k = frexp(y)[1]) - frexp(x)[1]) > 60)
t = big + big;
else if (m < -80)
t = y / x;
else
{
t = y / x;
y = ldexp(y, -k);
x = ldexp(x, -k);
}
}
/* begin argument reduction */
if (t < 2.4375)
{
/* truncate 4(t+1/16) to integer for branching */
k = (int) (4.0 * (t + 0.0625));
switch (k)
{
case 0: /* t is in [0, 7/16] */
case 1:
if (t < small)
return copysign((signx > 0.0) ? t : PI - t, signy);
hi = lo = 0.0;
break;
case 2: /* t is in [7/16, 11/16] */
hi = 4.6364760900080609352e-1;
lo = 4.6249969567426939759e-18;
z = x + x;
t = ((y + y) - x) / (z + y);
break;
case 3: /* t is in [11/16, 19/16] */
case 4:
hi = PIO4;
lo = 0.0;
t = (y - x) / (x + y);
break;
default: /* t is in [19/16, 39/16] */
hi = 9.8279372324732905408e-1;
lo = -2.4407677060164810007e-17;
z = y - x;
y = y + y + y;
t = x + x;
t = ((z + z) - x) / (t + y);
break;
}
}
else /* t >= 2.4375 */
{
hi = PIO2;
lo = 0.0;
if (t <= big) /* t is in [2.4375, big] */
t = -x / y;
else /* t is in [big, INF] */
t = 0.0;
}
/* compute atan(t) for t in [-0.4375, 0.4375] */
z = t * t;
z = t * (z * (3.3333333333333942106e-1 + z *
(-1.9999999999979536924e-1 + z *
(1.4285714278004377209e-1 + z *
(-1.1111110579344973814e-1 + z *
(9.0908906105474668324e-2 + z *
(-7.6919217767468239799e-2 + z *
(6.6614695906082474486e-2 + z *
(-5.8358371008508623523e-2 + z *
(4.9850617156082015213e-2 + z *
(-3.6700606902093604877e-2 + z *
1.6438029044759730479e-2)))))))))));
z = lo - z;
z += t;
z += hi;
return copysign((signx > 0.0) ? z : PI - z, signy);
}
# undef small
# undef big
/*
* NAME: atan()
* DESCRIPTION: return arctangent of argument
*/
static
float atan(float x)
{
return atan2(x, 1.0);
}
/* The following functions were borrowed from The Pattern (author: Dworkin) */
private
float poly(float x, float *coef)
{
int i;
float y;
for (i = sizeof(coef) - 1, y = coef[i]; i > 0; )
y = y * x + coef[--i];
return y;
}
private
float poly2(float x, float *coef1, float *coef2)
{
int i;
float y, z;
for (i = sizeof(coef1) - 1, y = coef1[i]; i > 0; )
y = y * x + coef1[--i];
for (i = sizeof(coef2) - 1, z = x + coef2[i]; i > 0; )
z = z * x + coef2[--i];
return y / z;
}
static
float exp(float x)
{
float n, xx, px;
if (x > MAXLOG)
error("Result too large");
if (x < MINLOG)
return 0.0;
/* Express e**x = e**g 2**n
* = e**g e**( n loge(2) )
* = e**( g + n loge(2) )
*/
n = floor(x * LOG2E + 0.5); /* floor() truncates toward -infinity. */
x -= n * 6.9335937500000000000E-1;
x += n * 2.1219444005469058277E-4;
/* 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 * poly(xx, ({ 1.00000000000000000000E0,
3.02996887658430129200E-2,
1.26183092834458542160E-4 }) );
x = px / (poly(xx, ({ 2.00000000000000000005E0,
2.27266044198352679519E-1,
2.52453653553222894311E-3,
3.00227947279887615146E-6 }) ) - px);
return ldexp(ldexp(x, 1) + 1.0, (int) n);
}
/* 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
/* 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
static
float log(float x)
{
mixed *fe;
int e;
float y, z;
/* Test for domain */
if (x <= 0.0)
{
if (x == 0.0)
error("Result too large");
else
error("Math argument");
}
/* separate mantissa from exponent */
fe = frexp(x);
x = fe[0];
e = fe[1];
if (e > 2 || e < -2)
{
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if (x < SQRTH)
{
/* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5;
y = ldexp(z, -1) + 0.5;
}
else
{
/* 2 (x-1)/(x+1) */
z = x - 0.5;
z -= 0.5;
y = ldexp(x, -1) + 0.5;
}
x = z / y;
/* rational form */
z = x * x;
z = x + x * (z * poly2(z,
({ -6.41409952958715622951E1,
1.63866645699558079767E1,
-7.89580278884799154124E-1 }),
({ -7.69691943550460008604E2,
3.12093766372244180303E2,
-3.56722798256324312549E1 }) ));
}
else
{
/* 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 * poly2(x,
({ 1.98892446572874072159E1,
5.67349287391754285487E1,
6.06127134467767258030E1,
2.97877425097986925891E1,
6.56312093769992875930E0,
4.98531067254050724270E-1,
4.58482948458143443514E-5 }),
({ 5.96677339718622216300E1,
2.14955586696422947765E2,
3.07254189979530058263E2,
2.20664384982121929218E2,
8.27410449222435217021E1,
1.50314182634250003249E1 }) ));
y = y - ldexp(z, -1); /* y - 0.5 * z */
z = x + y;
}
/* recombine with exponent term */
if (e != 0)
{
y = (float) e;
z = z - y * 2.121944400546905827679e-4;
z = z + y * 0.693359375;
}
return z;
}
/* 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
static
float log10(float x)
{
mixed *fe;
int e;
float w, y, z;
/* Test for domain */
if (x <= 0.0)
{
if (x == 0.0)
error("Result too large");
else
error("Math argument");
}
/* separate mantissa from exponent */
fe = frexp(x);
x = fe[0];
e = fe[1];
/* 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 * poly2(x,
({ 1.98892446572874072159E1,
5.67349287391754285487E1,
6.06127134467767258030E1,
2.97877425097986925891E1,
6.56312093769992875930E0,
4.98531067254050724270E-1,
4.58482948458143443514E-5 }),
({ 5.96677339718622216300E1,
2.14955586696422947765E2,
3.07254189979530058263E2,
2.20664384982121929218E2,
8.27410449222435217021E1,
1.50314182634250003249E1 }) ));
y = y - ldexp(z, -1); /* y - 0.5 * x**2 */
/* multiply log of fraction by log10(e)
* and base 2 exponent by log10(2)
*/
/* accumulate terms in order of size */
z = (x + y) * 7.00731903251827651129E-4; /* L10EB */
z += y * 4.3359375E-1; /* L10EA */
z += x * 4.3359375E-1; /* L10EA */
w = (float) e;
z += w * 2.48745663981195213739E-4; /* L102B */
z += w * 3.0078125E-1; /* L102A */
return z;
}
# if 0
/* 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, 30 2000 2.1e-17 5.2e-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.1: December, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
float sqrt(float x)
{
float w;
mixed *fe;
int e;
if (x <= 0.0) {
if (x < 0.0) {
error("Domain error");
}
return 0.0;
}
w = x;
/* separate exponent and significand */
fe = frexp(x);
x = fe[0];
e = fe[1];
/* approximate square root of number between 0.5 and 1
* relative error of approximation = 7.47e-3
*/
x = 4.173075996388649989089E-1 + 5.9016206709064458299663E-1 * x;
/* adjust for odd powers of 2 */
if (e & 1) {
x *= 1.41421356237309504880;
}
/* re-insert exponent */
x = ldexp(x, e / 2);
/* Newton iterations: */
x = 0.5 * (x + w / x);
x = 0.5 * (x + w / x);
x = 0.5 * (x + w / x);
return x;
}
# endif