/*
* math.c
*
* Mathematic funcs
*
* Most of these functions are NOT created by me, they are created by
* the guys behind MirrorMOO (for DGD) project.
* (C) Frank Schmidt, Jesus@NorseMUD
*
*/
#include <math.h>
/* return float squareroot */
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);
}
/* return sine or cosine of value */
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;
}
/* 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);
}
/* return cosine of argument */
static float cos(float x) {
return sin_cos(x, fabs(x) + PIO2, 1.0);
}
/* 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;
}
}
/* compute tangent of argument */
static float tan(float x) {
return tan_cot(x, 0);
}
/* 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);
}
private 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
/* 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
/* arctangent of argument */
static float atan(float x) {
return atan2(x, 1.0);
}
/* help-function to exp() */
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;
}
/* another help-function to exp() */
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;
}
/* compute exponential */
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);
}
/* compute logarithm */
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;
}
/* decimal logarithm */
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;
}