/* Copyright 1996 Acorn Computers Ltd
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
/* math.c: ANSI draft (X3J11 May 86) library code, section D.5 */
/* Copyright (C) Codemist Ltd, 1988 */
/* version 0.04b */
/* Nov 87: fix bug in ibm frexp(-ve arg). */
/*
* This version of the code takes the view that whenever there is an
* error a NaN should be handed back (as well as errno getting set). The
* value HUGE_VAL is used, which is not actually a NaN but which will
* often lead to exponent overflow pretty soon if it is used. ACN is
* unclear if this is sensible, and has had a program fall over when
* atan2(0.0, 0.0) handed back HUGE_VAL rather than some less vicious
* value (e.g. 0.0). He can imagine people who expect pow(0.0, 0.0) to
* be 1.0 (or maybe 0.0, but certainly not HUGE_VAL), and who expect
* sin(x) to be <= 1.0 in absolute value regardless of anything. Thus
* the current state is OK if we are being strict, but mey be unfriendly
* in some cases? Thoughts and comments, anybody?
*/
#include "hostsys.h"
#include <limits.h>
#include <float.h>
#include <errno.h>
/* This file contains code for most of the math routines from <math.h> */
/* On the ARM some of these routines are implemented as floating point */
/* opcodes and as such appear in startup.s */
#ifndef NO_FLOATING_POINT
#include <math.h> /* for forward references */
#include <fenv.h>
/* On the ARM, this has moved into the library's static data area */
/* so that it still works with the Shared C Library module. */
/* const double __huge_val = 1.79769313486231571e+308; */
#ifndef DO_NOT_SUPPORT_UNNORMALIZED_NUMBERS
double frexp(double d, int *lvn)
{
/* This version works even if d starts off as an unnormalized number in */
/* the IEEE sense. But in that special case it will be mighty slow! */
/* By that we mean at most 52 iterations for the smallest number. */
fp_number d1;
int n;
if (d==0.0)
{ *lvn = 0;
return d;
}
d1.d = d;
if ((n = d1.i.x - 0x3fe) == -0x3fe)
{ int flag;
/* Here d1 has zero in its exponent field - this means that the mantissa */
/* is un-normalized. I have to shift it left (at least one bit) until a */
/* suitable nonzero bit appears to go in the implicit-bit place in the */
/* fractional result. For each bit shifted I need to adjust the final */
/* exponent that will be returned. */
/* I have already tested to see if d was zero so the following loop MUST */
/* terminate. */
n++;
do
{ flag = d1.i.mhi & 0x00080000;
d1.i.mhi = (d1.i.mhi << 1) | (d1.i.mlo >> 31);
d1.i.mlo = d1.i.mlo << 1;
n--;
} while (flag==0);
}
else if (n == 0x401)
{
/* Here d1 has 0x7ff in its exponent field - it's a NaN or infinity */
*lvn = (d1.i.mhi || d1.i.mlo) ? FP_ILOGBNAN : INT_MAX;
return d1.d;
}
*lvn = n;
d1.i.x = 0x3fe;
return d1.d;
}
float frexpf(float s, int *lvn)
{
fp_number_single s1;
int n;
if (s==0.0F)
{ *lvn = 0;
return s;
}
s1.s = s;
if ((n = s1.i.x - 0x7e) == -0x7e)
{ int flag;
n++;
do
{ flag = s1.i.m & 0x00400000;
s1.i.m = s1.i.m << 1;
n--;
} while (flag==0);
}
else if (n == 0x81)
{
/* Here d1 has 0xff in its exponent field - it's a NaN or infinity */
*lvn = (s1.i.m) ? FP_ILOGBNAN : INT_MAX;
return s1.s;
}
*lvn = n;
s1.i.x = 0x7e;
return s1.s;
}
#else /* DO_NOT_SUPPORT_UNNORMALIZED_NUMBERS */
double frexp(double d, int *lvn)
{
fp_number d1;
if (d==0.0)
{ *lvn = 0;
return d;
}
d1.d = d;
*lvn = d1.i.x - 0x3fe;
d1.i.x = 0x3fe;
return d1.d;
}
float frexpf(float s, int *lvn)
{
fp_number_single s1;
int n;
if (s==0.0F)
{ *lvn = 0;
return s;
}
s1.s = s;
*lvn = s1.i.x - 0x7e;
s1.i.x = 0x7e;
return s1.s;
}
#endif /* DO_NOT_SUPPORT_UNNORMALIZED_NUMBERS */
double scalbln(double d, long int n)
{
fp_number x;
long int exponent;
if (n == 0 || d == 0.0) return d;
x.d = d;
exponent = x.i.x;
if (exponent == 0x7ff) return d; /* NaN or infinite */
#ifndef DO_NOT_SUPPORT_UNNORMALIZED_NUMBERS
if (exponent == 0) /* starting subnormal */
{ int flag;
exponent++;
do
{ flag = x.i.mhi & 0x00080000;
x.i.mhi = (x.i.mhi << 1) | (x.i.mlo >> 31);
x.i.mlo = x.i.mlo << 1;
exponent--;
} while (flag==0);
}
#endif
if (n > 0x1000 || exponent + n >= 0x7ff) /* overflow */
{ errno = ERANGE;
return x.i.s ? -HUGE_VAL : HUGE_VAL;
}
if (n < -0x1000 || exponent + n <= -53) /* total underflow */
{ x.i.x = x.i.mhi = x.i.mlo = 0;
errno = ERANGE;
return x.d;
}
if (exponent + n <= 0) /* subnormal result */
{ unsigned round = 0, hi;
n = 1 - (exponent + n);
hi = 0x00100000 | x.i.mhi;
if (n >= 32)
{ round = x.i.mlo;
x.i.mlo = hi;
hi = 0;
n -= 32;
}
if (n > 0)
{ round = (x.i.mlo << (32-n)) | (round != 0);
x.i.mlo = (hi << (32-n)) | (x.i.mlo >> n);
hi = hi >> n;
}
x.i.mhi = hi;
x.i.x = 0;
/* perform round-to-nearest; to even in tie cases */
if (round > 0x80000000 ||
(round == 0x80000000 && (x.i.mlo & 1)))
if (++x.i.mlo == 0)
if (++x.i.mhi == 0)
x.i.x = 1;
if (round != 0) errno = ERANGE; /* only if inexact */
return x.d;
}
/* normal result */
x.i.x = (int) (exponent + n);
return x.d;
}
float scalblnf(float s, long int n)
{
fp_number_single x;
long int exponent;
if (n == 0 || s == 0.0F) return s;
x.s = s;
exponent = x.i.x;
if (exponent == 0xff) return s; /* NaN or infinite */
#ifndef DO_NOT_SUPPORT_UNNORMALIZED_NUMBERS
if (exponent == 0) /* starting subnormal */
{ int flag;
exponent++;
do
{ flag = x.i.m & 0x00400000;
x.i.m = x.i.m << 1;
exponent--;
} while (flag==0);
}
#endif
if (n > 0x1000 || exponent + n >= 0xff) /* overflow */
{ errno = ERANGE;
return x.i.s ? -HUGE_VALF : HUGE_VALF;
}
if (n < -0x1000 || exponent + n <= -24) /* total underflow */
{ x.i.x = x.i.m = 0;
errno = ERANGE;
return x.s;
}
if (exponent + n <= 0) /* subnormal result */
{ unsigned round, m;
n = 1 - (exponent + n);
m = 0x00800000 | x.i.m;
round = m << (32-n);
x.i.m = m >> n;
x.i.x = 0;
/* perform round-to-nearest; to even in tie cases */
if (round > 0x80000000 ||
(round == 0x80000000 && (x.i.m & 1)))
if (++x.i.m == 0)
x.i.x = 1;
if (round != 0) errno = ERANGE; /* only if inexact */
return x.s;
}
/* normal result */
x.i.x = (int) (exponent + n);
return x.s;
}
double scalbn(double x, int n)
{
return scalbln(x, n);
}
float scalbnf(float x, int n)
{
return scalblnf(x, n);
}
double ldexp(double x, int n)
{
return scalbn(x, n);
}
float ldexpf(float x, int n)
{
return scalbnf(x, n);
}
int ilogb(double d)
{
fp_number x;
int exponent;
if (d == 0.0)
{ errno = ERANGE;
return FP_ILOGB0;
}
x.d = d;
exponent = x.i.x;
#ifndef DO_NOT_SUPPORT_UNNORMALIZED_NUMBERS
if (exponent == 0)
{ int flag;
exponent = 1;
do
{ flag = x.i.mhi & 0x00080000;
x.i.mhi = (x.i.mhi << 1) | (x.i.mlo >> 31);
x.i.mlo = x.i.mlo << 1;
exponent--;
} while (flag==0);
}
else if (exponent == 0x7ff)
{ if (x.i.mhi || x.i.mlo)
return FP_ILOGBNAN;
else
return INT_MAX;
}
#endif
return exponent - 0x3ff;
}
int ilogbf(float s)
{
fp_number_single x;
int exponent;
if (s == 0.0F)
{ errno = ERANGE;
return FP_ILOGB0;
}
x.s = s;
exponent = x.i.x;
#ifndef DO_NOT_SUPPORT_UNNORMALIZED_NUMBERS
if (exponent == 0)
{ int flag;
exponent = 1;
do
{ flag = x.i.m & 0x00400000;
x.i.m = x.i.m << 1;
exponent--;
} while (flag==0);
}
else if (exponent == 0xff)
{ if (x.i.m)
return FP_ILOGBNAN;
else
return INT_MAX;
}
#endif
return exponent - 0x7f;
}
double logb(double d)
{
int x = ilogb(d);
switch (x)
{
case FP_ILOGB0: return -HUGE_VAL;
case FP_ILOGBNAN: return d;
case INT_MAX: return INFINITY;
default: return x;
}
}
float logbf(float s)
{
int x = ilogbf(s);
switch (x)
{
case FP_ILOGB0: return -HUGE_VALF;
case FP_ILOGBNAN: return s;
case INT_MAX: return INFINITY;
default: return x;
}
}
#define _exp_arg_limit 709.78271289338397
/* machine independent code - but beware the 1e-20's */
#define _pi_ 3.14159265358979323846
#define _pi_3 1.04719755119659774615
#define _pi_2 1.57079632679489661923
#define _pi_4 0.78539816339744830962
#define _pi_6 0.52359877559829887038
#define _sqrt_half 0.70710678118654752440
#ifndef HOST_HAS_TRIG
/* The ARM has a load of trig functions etc available as opcodes, so */
/* has no need for these software versions shown here. */
#undef sin
#undef cos
#undef atan
static double _sincos(double x, double y, int sign, int coscase)
{
int n;
double xn, f, g, r;
if (y >= 1.0e9) /* fail if argument is overlarge */
{ errno = EDOM;
return -HUGE_VAL;
}
n = (int) ((y + _pi_2) / _pi_);
xn = (double) n;
if ((n & 1) != 0) sign = -sign;
if (coscase) xn = xn - 0.5;
/* Note that these days C is REQUIRED to observe the brackets used here. */
/* The compiler is NOT allowed to re-associate the additions. */
#ifdef IBMFLOAT
{ double x1 = (double)(int)x;
/* observe that the range check on y assures me that (int)x is OK. */
double x2 = x - x1;
f = ((x1 - xn*3.1416015625) + x2) + xn*8.908910206761537356617e-6;
}
#else
f = (x - xn*3.1416015625) + xn*8.908910206761537356617e-6;
#endif
/* I expect that the absolute value of f is less than pi/2 here */
if (fabs(f) >= 1.e-10)
#define _sincos_r1 -0.16666666666666665052
#define _sincos_r2 0.83333333333331650315e-2
#define _sincos_r3 -0.19841269841201840457e-3
#define _sincos_r4 0.27557319210152756119e-5
#define _sincos_r5 -0.25052106798274584544e-7
#define _sincos_r6 0.16058936490371589114e-9
#define _sincos_r7 -0.76429178068910467734e-12
#define _sincos_r8 0.27204790957888846175e-14
{ g = f*f;
r = ((((((((_sincos_r8) * g + _sincos_r7) * g + _sincos_r6) * g +
_sincos_r5) * g + _sincos_r4) * g + _sincos_r3) * g +
_sincos_r2) * g + _sincos_r1) * g;
f += f*r;
};
if (sign < 0) return -f;
else return f;
}
double sin(double x)
{
if (x < 0.0) return _sincos(-x, -x, -1, 0);
else return _sincos(x, x, 1, 0);
}
double cos(double x)
{
if (x < 0.0) return _sincos(-x, _pi_2 - x, 1, 1);
else return _sincos(x, _pi_2 + x, 1, 1);
}
/* NB: the following value is a proper limit provided denormalised */
/* values are not being supported. It would need to be changed if they */
/* were to start existing. */
#define _exp_negative_arg_limit -708.39641853226408
double exp(double x)
{
int n;
double xn, g, z, gp, q, r;
if (x > _exp_arg_limit)
{ errno = ERANGE;
return HUGE_VAL;
}
if (x < _exp_negative_arg_limit) return 0.0;
if (fabs(x) < 1.e-20) return 1.0;
/* In C the cast (int)x truncates towards zero. Here I want to round. */
n = (int)((x >= 0 ? 0.5 : -0.5) + 1.44266950408889634074 * x);
xn = (double)n;
g = (x - xn * 0.693359375) - xn * (-2.1219444005469058277e-4);
z = g * g;
#define _exp_p0 0.249999999999999993
#define _exp_p1 0.694360001511792852e-2
#define _exp_p2 0.165203300268279130e-4
#define _exp_q0 0.500000000000000000
#define _exp_q1 0.555538666969001188e-1
#define _exp_q2 0.495862884905441294e-3
gp = ((_exp_p2 * z + _exp_p1) * z + _exp_p0) * g;
q = (_exp_q2 * z + _exp_q1) * z + _exp_q0;
r = 0.5 + gp / (q - gp);
return ldexp(r, n + 1);
}
double log(double x)
{
double f, znum, zden, z, w, r, xn;
int n;
if (x <= 0.0)
{ if (x==0.0)
{ errno = ERANGE;
return -HUGE_VAL;
} else
{ errno = EDOM;
return -HUGE_VAL;
}
}
f = frexp(x, &n);
if (f > _sqrt_half)
{ znum = (f - 0.5) - 0.5;
zden = f * 0.5 + 0.5;
}
else
{ n -= 1;
znum = f - 0.5;
zden = znum * 0.5 + 0.5;
}
z = znum / zden;
w = z * z;
#define _log_a0 -0.64124943423745581147e2
#define _log_a1 0.16383943563021534222e2
#define _log_a2 -0.78956112887491257267e0
#define _log_b0 -0.76949932108494879777e3
#define _log_b1 0.31203222091924532844e3
#define _log_b2 -0.35667977739034646171e2
r = w * ( ((_log_a2 * w + _log_a1) * w + _log_a0) /
(((w + _log_b2) * w + _log_b1) * w + _log_b0) );
r = z + z * r;
xn = (double)n;
return (xn*(-2.121944400546905827679e-4) + r) + xn*(355.0/512.0);
}
double log10(double x)
{
return log(x) * 0.43429448190325182765; /* log10(e) */
}
double sqrt(double x)
{
fp_number f;
double y0;
int n;
if (x <= 0.0)
{ if (x < 0.0) errno = EDOM;
return -HUGE_VAL;
}
f.d = x;
n = f.i.x - 0x3fe;
f.i.x = 0x3fe;
{ double fd = f.d;
y0 = 0.41731 + 0.59016 * fd;
y0 = 0.5 * (y0 + fd/y0);
y0 = 0.5 * (y0 + fd/y0);
y0 = 0.5 * (y0 + fd/y0);
}
if (n & 1)
{
y0 *= _sqrt_half;
n += 1;
}
n /= 2;
f.d = y0;
f.i.x += n;
return f.d;
}
double _tancot(double x, int iflag)
{
int n;
double f, g, xnum, xden, y, xn;
y = fabs(x);
if (y >= 1.0e9) /* fail if argument is overlarge */
{ errno = EDOM;
return -HUGE_VAL;
}
n = (int) ((2.0 * y + _pi_2) / _pi_);
if (x < 0) n = - n;
xn = (double) n;
f = (x - xn*1.57080078125) + xn*4.454455103380768678308e-6;
if (fabs(f) > 1.e-10)
{ g = f * f;
#define _tan_p1 -0.13338350006421960681
#define _tan_p2 0.34248878235890589960e-2
#define _tan_p3 -0.17861707342254426711e-4
#define _tan_q0 1.00000000000000000000
#define _tan_q1 -0.46671683339755294240
#define _tan_q2 0.25663832289440112864e-1
#define _tan_q3 -0.31181531907010027307e-3
#define _tan_q4 0.49819433993786512270e-6
xnum = ((_tan_p3*g + _tan_p2)*g + _tan_p1)*g*f + f;
xden = (((_tan_q4*g + _tan_q3)*g + _tan_q2)*g + _tan_q1)*g + _tan_q0;
}
else
{ xnum = f;
xden = 1.0;
}
/* It is probable that overflow can never occur here, since floating */
/* point values fall about or more 1.e-16 apart near singularities */
/* of the tangent function. */
if (iflag==0)
{ if ((n & 1) == 0) return xnum / xden;
else return - xden / xnum;
}
else
{ if ((n & 1) == 0) return xden / xnum;
else return - xnum / xden;
}
}
double tan(double x)
{
return _tancot(x, 0);
}
double _cot(double x) /* Not specified by ANSI hence the funny name */
{
if (fabs(x) < 1.0/HUGE_VAL)
{ errno = ERANGE;
if (x < 0.0) return -HUGE_VAL;
else return HUGE_VAL;
}
return _tancot(x, 1);
}
double atan(double x)
{
int n;
double f;
const static double a[4] = { 0.0, _pi_6, _pi_2, _pi_3 };
f = fabs(x);
if (f > 1.0)
{ f = 1.0 / f;
n = 2;
}
else n = 0;
#define _two_minus_root_three 0.26794919243112270647
#define _sqrt_three 1.73205080756887729353
#define _sqrt_three_minus_one 0.73205080756887729353
if (f > _two_minus_root_three)
{ f = (((_sqrt_three_minus_one*f - 0.5) - 0.5) + f) / (_sqrt_three + f);
n++;
}
if (fabs(f) > 1.e-10)
{ double g = f * f;
double r;
#define _atan_p0 -0.13688768894191926929e2
#define _atan_p1 -0.20505855195861651981e2
#define _atan_p2 -0.84946240351320683534e1
#define _atan_p3 -0.83758299368150059274
#define _atan_q0 0.41066306682575781263e2
#define _atan_q1 0.86157349597130242515e2
#define _atan_q2 0.59578436142597344465e2
#define _atan_q3 0.15024001160028576121e2
r = ((((_atan_p3*g + _atan_p2)*g + _atan_p1)*g + _atan_p0)*g) /
((((g + _atan_q3)*g + _atan_q2)*g + _atan_q1)*g + _atan_q0);
f = f + f * r;
}
if (n > 1) f = -f;
f = f + a[n];
if (x < 0) return -f;
else return f;
}
double _asinacos(double x, int flag)
{
int i;
double y, g, r;
const static double a[2] = {0.0, _pi_4 };
const static double b[2] = {_pi_2, _pi_4 };
y = fabs(x);
if (y < 1.e-10) i = flag;
else
{ if (y > 0.5)
{ i = 1 - flag;
if (y > 1.0)
{ errno = EDOM;
return -HUGE_VAL;
}
g = ((0.5 - y) + 0.5) * 0.5;
y = -2.0 * sqrt(g);
}
else
{ i = flag;
g = y * y;
}
#define _asin_p1 -0.27368494524164255994e2
#define _asin_p2 0.57208227877891731407e2
#define _asin_p3 -0.39688862997504877339e2
#define _asin_p4 0.10152522233806463645e2
#define _asin_p5 -0.69674573447350646411e0
#define _asin_q0 -0.16421096714498560795e3
#define _asin_q1 0.41714430248260412556e3
#define _asin_q2 -0.38186303361750149284e3
#define _asin_q3 0.15095270841030604719e3
#define _asin_q4 -0.23823859153670238830e2
r = (((((_asin_p5*g + _asin_p4)*g + _asin_p3)*g +
_asin_p2)*g + _asin_p1)*g) /
(((((g + _asin_q4)*g + _asin_q3)*g +
_asin_q2)*g + _asin_q1)*g + _asin_q0);
y = y + y*r;
}
if (flag==0)
{ y = (a[i] + y) + a[i];
if (x<0) y = -y;
}
else if (x < 0) y = (b[i] + y) + b[i];
else y = (a[i] - y) + a[i];
return y;
}
double asin(double x)
{
return _asinacos(x, 0);
}
double acos(double x)
{
return _asinacos(x, 1);
}
double pow(double x, double y)
{
int sign = 0, m, p, i, mdash, pdash;
double g, r, z, v, u1, u2;
/* The table a1[] contains properly rounded values for 2**(i/16), and */
/* a2[] contains differences between the true values of 2**(i/16) and */
/* the a1[] values for odd i. */
const static double a1[17] = {
/* It is painfully important that the following 17 floating point */
/* numbers are read in to yield the quantities shown on the right. */
1.0, /* 3ff00000:00000000 */
0.9576032806985737, /* 3feea4af:a2a490da */
0.91700404320467121, /* 3fed5818:dcfba487 */
0.87812608018664973, /* 3fec199b:dd85529c */
0.8408964152537145, /* 3feae89f:995ad3ad */
0.80524516597462714, /* 3fe9c491:82a3f090 */
0.77110541270397037, /* 3fe8ace5:422aa0db */
0.73841307296974967, /* 3fe7a114:73eb0187 */
0.70710678118654757, /* 3fe6a09e:667f3bcd */
0.67712777346844632, /* 3fe5ab07:dd485429 */
0.64841977732550482, /* 3fe4bfda:d5362a27 */
0.620928906036742, /* 3fe3dea6:4c123422 */
0.59460355750136051, /* 3fe306fe:0a31b715 */
0.56939431737834578, /* 3fe2387a:6e756238 */
0.54525386633262884, /* 3fe172b8:3c7d517b */
0.52213689121370688, /* 3fe0b558:6cf9890f */
0.5 /* 3fe00000:00000000 */
};
const static double a2[8] = {
-5.3099730280915798e-17,
1.4800703477186887e-17,
1.2353596284702225e-17,
-1.7419972784343138e-17,
3.8504741898901863e-17,
2.3290137959059465e-17,
4.4563878092065126e-17,
4.2759448527536718e-17
};
if (y == 1.0) return x;
if (x <= 0.0)
{ int ny;
if (x==0.0)
{ if (y <= 0.0)
{ errno = EDOM;
if (y==0.0) return 1.0;
return -HUGE_VAL;
}
return 0.0;
}
if (y < (double)INT_MIN || y > (double)INT_MAX ||
(double)(ny = (int)y) != y)
{ errno = EDOM;
return -HUGE_VAL;
}
/* Here y is an integer and x is negative. */
x = -x;
sign = (ny & 1);
}
if (y == 2.0 && x < 1.e20) return x*x; /* special case. */
g = frexp(x, &m);
p = 0;
if (g <= a1[8]) p = 8;
if (g <= a1[p+4]) p += 4;
if (g <= a1[p+2]) p += 2;
z = ((g - a1[p+1]) - a2[p/2]) / (0.5*g + 0.5*a1[p+1]);
/* Expect abs(z) <= 0.044 here */
v = z * z;
#define _pow_p1 0.83333333333333211405e-1
#define _pow_p2 0.12500000000503799174e-1
#define _pow_p3 0.22321421285924258967e-2
#define _pow_p4 0.43445775672163119635e-3
r = (((_pow_p4*v + _pow_p3)*v + _pow_p2)*v + _pow_p1)*v*z;
#define _pow_k 0.44269504088896340736
r = r + _pow_k * r;
u2 = (r + z * _pow_k) + z;
#define _reduce(v) ((double)(int)(16.0*(v))*0.0625)
u1 = (double)(16*m-p-1) * 0.0625;
{ double y1 = _reduce(y);
double y2 = y - y1;
double w = u2*y + u1*y2;
double w1 = _reduce(w);
double w2 = w - w1;
int iw1;
w = w1 + u1*y1;
w1 = _reduce(w);
w2 = w2 + (w - w1);
w = _reduce(w2);
iw1 = (int)(16.0*(w1+w));
w2 = w2 - w;
/* The following values have been determined experimentally, buth their */
/* values are not very critical. */
# define _negative_pow_limit -16352
if (iw1 < _negative_pow_limit)
{ errno = ERANGE; /* Underflow certain */
return 0.0;
}
if (w2 > 0.0)
{ iw1 += 1;
w2 -= 0.0625;
}
if (iw1 < 0) i = 0;
else i = 1;
mdash = iw1/16 + i;
pdash = 16*mdash - iw1;
#define _pow_q1 0.69314718055994529629
#define _pow_q2 0.24022650695909537056
#define _pow_q3 0.55504108664085595326e-1
#define _pow_q4 0.96181290595172416964e-2
#define _pow_q5 0.13333541313585784703e-2
#define _pow_q6 0.15400290440989764601e-3
#define _pow_q7 0.14928852680595608186e-4
z = ((((((_pow_q7*w2 + _pow_q6)*w2 + _pow_q5)*w2 +
_pow_q4)*w2 + _pow_q3)*w2 + _pow_q2)*w2 + _pow_q1)*w2;
z = a1[pdash] + a1[pdash]*z;
z = frexp(z, &m);
mdash += m;
if (mdash >= 0x7ff-0x3fe)
{ errno = ERANGE;
if (sign) r = -HUGE_VAL;
else r = HUGE_VAL;
}
else if (mdash <= -0x3fe)
{ errno = ERANGE;
r = 0.0;
}
else
{ r = ldexp(z, mdash);
if (sign) r = -r;
}
}
return r;
}
#endif /* HOST_HAS_TRIG */
double atan2(double y, double x)
{
if (isunordered(y, x)) return y+x;
if (y==0.0 || (isinf(x) && isfinite(y)))
return signbit(x) ? copysign(_pi_, y) : y;
if (isinf(x)) /* y must also be infinite, from above */
return copysign(x > 0 ? _pi_4 : 3*_pi_4, y);
if (fabs(x) < fabs(y)) /* handles infinite y, finite x */
{ if (fabs(x / y)<1.0e-20)
{ if (y<0.0) return - _pi_2;
else return _pi_2;
}
}
y = atan(y / x);
if (x<0.0)
{ if (y>0.0) y -= _pi_;
else y += _pi_;
}
return y;
}
float atan2f(float y, float x)
{
if (isunordered(y, x)) return y+x;
if (y==0.0F || (isinf(x) && isfinite(y)))
return signbit(x) ? copysignf((float)_pi_, y) : y;
if (isinf(x)) /* y must also be infinite, from above */
return copysignf(x > 0 ? (float)_pi_4 : (float)(3*_pi_4), y);
if (fabsf(x) < fabsf(y)) /* handles infinite y, finite x */
{ if (fabsf(x / y)<1.0e-12F)
{ if (y<0.0F) return (float) - _pi_2;
else return (float) _pi_2;
}
}
y = atanf(y / x);
if (x<0.0F)
{ if (y>0.0F) y -= (float) _pi_;
else y += (float) _pi_;
}
return y;
}
#ifndef HOST_HAS_TRIG
double hypot(double x, double y)
{
#pragma STDC FENV_ACCESS ON
int excepts;
double p, q, r;
fenv_t env;
p = fabs(x);
q = fabs(y);
if (isunordered(p, q))
{
if (p == INFINITY || q == INFINITY)
return INFINITY;
else
return p + q;
}
if (isless(p, q))
{
r = p; p = q; q = r;
}
if (p == 0 || q == INFINITY)
return q;
feholdexcept(&env);
r = q / p;
p = p * sqrt(1 + r * r);
excepts = fetestexcept(FE_OVERFLOW);
fesetenv(&env);
if (excepts)
{
errno = ERANGE;
return HUGE_VAL;
}
return p;
}
/* All over again, in floats */
float hypotf(float x, float y)
{
#pragma STDC FENV_ACCESS ON
int excepts;
float p, q; double r;
fenv_t env;
float sqrtf(float);
p = fabsf(x);
q = fabsf(y);
if (isunordered(p, q))
{
if (p == INFINITY || q == INFINITY)
return INFINITY;
else
return p + q;
}
if (isless(p, q))
{
r = p; p = q; q = r;
}
if (p == 0 || q == INFINITY)
return q;
feholdexcept(&env);
r = (double) q / p;
p = (float) (p * sqrt(1 + r * r));
excepts = fetestexcept(FE_OVERFLOW);
fesetenv(&env);
if (excepts)
{
errno = ERANGE;
return HUGE_VAL;
}
return p;
}
#endif
/* Use the macros, which expand to inline ABS instructions */
double (fabs)(double x) { return fabs(x); }
float (fabsf)(float x) { return fabsf(x); }
double sinh(double x)
{
int sign;
double y;
if (x==0) return x;
if (isless(x,0.0)) y = -x, sign = 1; else y = x, sign = 0;
if (isgreater(y,1.0))
{
/* _sinh_lnv is REQUIRED to read in as a number with the lower part of */
/* its floating point representation zero. */
#define _sinh_lnv 0.69316101074218750000 /* ln(v) */
#define _sinh_vm2 0.24999308500451499336 /* 1/v^2 */
#define _sinh_v2m1 0.13830277879601902638e-4 /* (v/2) - 1 */
double w = y - _sinh_lnv, z, r;
if (w>_exp_arg_limit)
{ if (isinf(x)) return x;
errno = ERANGE;
if (sign) return -HUGE_VAL;
else return HUGE_VAL;
}
z = exp(w); /* should not overflow! */
if (z < 1.0e10) z = z - _sinh_vm2/z;
r = z + _sinh_v2m1 * z;
if (sign) return -r;
else return r;
}
else if (!isgreater(y,1.0e-10)) return x;
else
{
#define _sinh_p0 -0.35181283430177117881e6
#define _sinh_p1 -0.11563521196851768270e5
#define _sinh_p2 -0.16375798202630751372e3
#define _sinh_p3 -0.78966127417357099479e0
#define _sinh_q0 -0.21108770058106271242e7
#define _sinh_q1 0.36162723109421836460e5
#define _sinh_q2 -0.27773523119650701667e3
#define _sinh_q3 1.0
double g = x*x;
double r;
/* Use a (minimax) rational approximation. See Cody & Waite. */
r = ((((_sinh_p3*g + _sinh_p2)*g + _sinh_p1)*g + _sinh_p0)*g) /
(((g + _sinh_q2)*g + _sinh_q1)*g + _sinh_q0);
return x + x*r;
}
}
double cosh(double x)
{
x = fabs(x);
if (isgreater(x,1.0))
{ x = x - _sinh_lnv;
if (x>_exp_arg_limit)
{ if (isinf(x)) return x;
errno = ERANGE;
return HUGE_VAL;
}
x = exp(x); /* the range check above ensures that this does */
/* not overflow. */
if (x < 1.0e10) x = x + _sinh_vm2/x;
/* This very final line might JUST overflow even though the call */
/* to exp is safe and even though _exp_arg_limit is conservative */
return x + _sinh_v2m1 * x;
}
/* This second part is satisfactory, even though it is simple! */
x = exp(x);
return 0.5*(x + 1/x);
}
double tanh(double x)
{
/* The first two exits avoid premature overflow as well as needless use */
/* of the exp() function. */
int sign;
if (isgreater(x,27.0)) return 1.0; /* here exp(2x) dominates 1.0 */
else if (isless(x,-27.0)) return -1.0;
if (isless(x,0.0)) x = -x, sign = 1;
else sign = 0;
if (isgreater(x,0.549306144334054846)) /* ln(3)/2 is crossover point */
{ x = exp(2.0*x);
x = 1.0 - 2.0/(x + 1.0);
if (sign) return -x;
else return x;
}
#define _tanh_p0 -0.16134119023996228053e4
#define _tanh_p1 -0.99225929672236083313e2
#define _tanh_p2 -0.96437492777225469787e0
#define _tanh_q0 0.48402357071988688686e4
#define _tanh_q1 0.22337720718962312926e4
#define _tanh_q2 0.11274474380534949335e3
#define _tanh_q3 1.0
if (isgreater(x,1.0e-10))
{ double y = x*x;
/* minimax rational approximation */
y = (((_tanh_p2*y + _tanh_p1)*y + _tanh_p0)*y) /
(((y + _tanh_q2)*y + _tanh_q1)*y + _tanh_q0);
x = x + x*y;
}
if (sign) return -x;
else return x;
}
/* Simple forms for now */
float coshf(float x)
{
return (float)cosh(x);
}
float sinhf(float x)
{
return (float)sinh(x);
}
float tanhf(float x)
{
return (float)tanh(x);
}
double asinh(double x)
{
if (isless(x, 0.0))
return -asinh(-x);
if (x == 0) return x; // asinh(�0) returns �0
return log(x+hypot(x,1.0));
}
double acosh(double x)
{
if (isless(x, 1.0))
{ errno = EDOM;
return NAN;
}
// acosh(0) = 2.0*log(sqrt(1)+sqrt(0)) = 2.0*log(1) = +0
return 2.0*log(sqrt((x+1)*0.5)+sqrt((x-1)*0.5));
}
double atanh(double x)
{
if (isless(x, 0.0))
return -atanh(-x);
if (x == 0) return x; // atanh(�0) returns �0
if (isgreater(x, 1.0))
{ errno = EDOM;
return NAN;
}
if (x == 1) // atanh(�1) returns �inf
{ errno = ERANGE;
return HUGE_VAL;
}
return 0.5*log((1+x)/(1-x));
}
float asinhf(float x)
{
if (isless(x, 0.0F))
return -asinhf(-x);
if (x == 0) return x; // asinhf(�0) returns �0
return (float) log(x+hypot(x,1.0));
}
float acoshf(float x)
{
if (isless(x, 1.0F))
{ errno = EDOM;
return NAN;
}
// acoshf(0) = 2.0*log(sqrt(1)+sqrt(0)) = 2.0*log(1) = +0
return (float) (2.0*log(sqrt((x+1.0)*0.5)+sqrt((x-1.0)*0.5)));
}
float atanhf(float x)
{
if (isless(x, 0.0F))
return -atanhf(-x);
if (x == 0) return x; // atanhf(�0) returns �0
if (isgreater(x, 1.0))
{ errno = EDOM;
return NAN;
}
if (x == 1) // atanhf(�1) returns �inf
{ errno = ERANGE;
return HUGE_VALF;
}
return (float) (0.5*log((1.0+x)/(1.0-x)));
}
double fmod(double x, double y)
{
/* floating point remainder of (x/y) for integral quotient. Remainder */
/* has same sign as x. */
double q, r;
if (y==0.0)
{
errno = EDOM;
return -HUGE_VAL;
}
if (x==0.0) return 0.0;
if (y < 0.0) y = -y;
r = modf(x/y, &q);
r = x - q * y;
/* The next few lines are an ultra-cautious scheme to ensure that the */
/* result is less than fabs(y) in value and that it has the sign of x. */
if (x > 0.0)
{ while (r >= y) r -= y;
while (r < 0.0) r += y;
}
else
{ while (r <= -y) r += y;
while (r > 0.0) r -= y;
}
return r;
}
float fmodf(float x, float y)
{
/* floating point remainder of (x/y) for integral quotient. Remainder */
/* has same sign as x. */
float q, r;
if (y==0.0F)
{
errno = EDOM;
return -HUGE_VALF;
}
if (x==0.0F) return x;
if (y < 0.0F) y = -y;
r = modff(x/y, &q);
r = x - q * y;
/* The next few lines are an ultra-cautious scheme to ensure that the */
/* result is less than fabs(y) in value and that it has the sign of x. */
if (x > 0.0F)
{ while (r >= y) r -= y;
while (r < 0.0F) r += y;
}
else
{ while (r <= -y) r += y;
while (r > 0.0F) r -= y;
}
return r;
}
#if 0
double (floor)(double d)
{
/* round x down to an integer towards minus infinity. */
fp_number x;
int exponent, mask, exact;
if (d == 0.0) return 0.0;
x.d = d; /* pun on union type */
if ((exponent = x.i.x - 0x3ff) < 0)
{ if (x.i.s) return -1.0;
else return 0.0;
}
else if (exponent >= 52) return x.d;
if (exponent >= 20)
{ mask = ((unsigned) 0xffffffff) >> (exponent - 20);
exact = x.i.mlo & mask;
x.i.mlo &= ~mask;
}
else
{ mask = 0xfffff >> exponent;
exact = (x.i.mhi & mask) | x.i.mlo;
x.i.mhi &= ~mask;
x.i.mlo = 0;
}
if (exact!=0 && x.i.s) return x.d - 1.0;
else return x.d;
}
double (ceil)(double d)
{
/* round x up to an integer towards plus infinity. */
fp_number x;
int exponent, mask, exact;
if (d == 0.0) return d;
x.d = d; /* pun on union type */
if ((exponent = x.i.x - 0x3ff) < 0)
{ if (x.i.s) return 0.0;
else return 1.0;
}
else if (exponent >= 52) return x.d;
if (exponent >= 20)
{ mask = ((unsigned) 0xffffffff) >> (exponent - 20);
exact = x.i.mlo & mask;
x.i.mlo &= ~mask;
}
else
{ mask = 0xfffff >> exponent;
exact = (x.i.mhi & mask) | x.i.mlo;
x.i.mhi &= ~mask;
x.i.mlo = 0;
}
if (exact!=0 && x.i.s==0) return x.d + 1.0;
else return x.d;
}
#endif
double (ceil)(double x) { return ceil(x); }
float (ceilf)(float x) { return ceilf(x); }
double (floor)(double x) { return floor(x); }
float (floorf)(float x) { return floorf(x); }
double (rint)(double x) { return rint(x); }
float (rintf)(float x) { return rintf(x); }
long int (lrint)(double x) { return lrint(x); }
long int (lrintf)(float x) { return lrintf(x); }
double (trunc)(double x) { return trunc(x); }
float (truncf)(float x) { return truncf(x); }
double (atan)(double x) { return atan(x); }
float (atanf)(float x) { return atanf(x); }
double (sin)(double x) { return sin(x); }
float (sinf)(float x) { return sinf(x); }
double (cos)(double x) { return cos(x); }
float (cosf)(float x) { return cosf(x); }
double modf(double value, double *iptr)
{
/* splits value into integral part & fraction (both same sign) */
fp_number x;
int exponent, mask;
if (value == 0.0)
{ *iptr = value;
return value;
}
x.d = value;
if ((exponent = x.i.x - 0x3ff) < 0)
{ *iptr = copysign(0.0, value);
return value;
}
else if (exponent >= 52)
{ *iptr = value;
return copysign(0.0, value);
}
if (exponent >= 20)
{ mask = ((unsigned) 0xffffffff) >> (exponent - 20);
x.i.mlo &= ~mask;
}
else
{ mask = 0xfffff >> exponent;
x.i.mhi &= ~mask;
x.i.mlo = 0;
}
*iptr = x.d;
return value - x.d;
}
float modff(float value, float *iptr)
{
/* splits value into integral part & fraction (both same sign) */
fp_number_single x;
int exponent, mask;
if (value == 0.0F)
{ *iptr = value;
return value;
}
x.s = value;
if ((exponent = x.i.x - 0xff) < 0)
{ *iptr = copysignf(0.0F, value);
return value;
}
else if (exponent >= 23)
{ *iptr = value;
return copysignf(0.0F, value);
}
mask = 0x7fffff >> exponent;
x.i.m &= ~mask;
*iptr = x.s;
return value - x.s;
}
double nan(const char *s)
{
return (double) NAN;
}
float nanf(const char *s)
{
return NAN;
}
double fdim(double x, double y)
{
if (islessequal(x, y))
return 0.0;
else
return x - y; /* Will return NaN for NaN input */
}
float fdimf(float x, float y)
{
if (islessequal(x, y))
return 0.0F;
else
return x - y; /* Will return NaN for NaN input */
}
#define _sqrt2pi 2.50662827463100050242
/*
* Gamma functions computed using Lanczos' approximation.
* Double version uses coefficients computed as per Spouge (1994)
* to achieve (theoretically) < 1 ulp error.
* Float version uses original Lanczos coefficients.
*
* Lanczos' approximation:
*
* G(x+1) = (x+a)^(x+.5)* e^-(x+a) * sqrt(2*pi) *
* [ c0 + c1/(x+1) + c2/(x+3).. + cn/(x+n) ] (x > 0)
*
*
* Spouge's coefficients:
* c0 = 1
* c[k] = (-1)^(k-1) * e^(a-k) * (a-k)^(k-0.5) ( 1 <= k < ceil(a) )
* ------------------------------------
* sqrt(2*pi) * (k-1)!
*
* giving relative error < sqrt(a) * (2*pi)^-(a+0.5) / (a+x)
*
* Useful relations: gamma(x+1) = x*gamma(x)
* gamma(1-x) = pi / (gamma(x)*sin(pi*x))
* gamma(n+1) = n!
*/
#define _lgamma_c1 166599025.853949811570 // actually c[k]*sqrt(2*pi)
#define _lgamma_c2 -981939810.195007931170
#define _lgamma_c3 2548964102.46316408700
#define _lgamma_c4 -3836248618.43839895348
#define _lgamma_c5 3709080184.61731236844
#define _lgamma_c6 -2412708472.49486138749
#define _lgamma_c7 1075449464.75190680642
#define _lgamma_c8 -328534715.011179056348
#define _lgamma_c9 67752870.4715251633277
#define _lgamma_c10 -9145761.20044444915581
#define _lgamma_c11 768402.637723269577278
#define _lgamma_c12 -37175.9448951564986832
#define _lgamma_c13 917.944248521710658895
#define _lgamma_c14 -9.51510944794615044564
#define _lgamma_c15 0.0294177178100457006509
#define _lgamma_c16 -1.37185031730621246722e-5
#define _lgamma_c17 1.72316912091954830013e-10
#define _lgamma_c18 -2.50065513100139951901e-20
static inline double _gamma_sum(double x)
{
/* Do it like this to avoid using client static data for an array */
return _sqrt2pi
+ _lgamma_c1 / (x + 1)
+ _lgamma_c2 / (x + 2)
+ _lgamma_c3 / (x + 3)
+ _lgamma_c4 / (x + 4)
+ _lgamma_c5 / (x + 5)
+ _lgamma_c6 / (x + 6)
+ _lgamma_c7 / (x + 7)
+ _lgamma_c8 / (x + 8)
+ _lgamma_c9 / (x + 9)
+ _lgamma_c10 / (x + 10)
+ _lgamma_c11 / (x + 11)
+ _lgamma_c12 / (x + 12)
+ _lgamma_c13 / (x + 13)
+ _lgamma_c14 / (x + 14)
+ _lgamma_c15 / (x + 15)
+ _lgamma_c16 / (x + 16)
+ _lgamma_c17 / (x + 17)
+ _lgamma_c18 / (x + 18);
}
#define _lgammaf_c0 1.000000000190015
#define _lgammaf_c1 76.18009172947146
#define _lgammaf_c2 -86.50532032941677
#define _lgammaf_c3 24.01409824083091
#define _lgammaf_c4 -1.231739572450155
#define _lgammaf_c5 1.208650973866179e-3
#define _lgammaf_c6 -5.395239384953e-6
static inline double _gammaf_sum(double x)
{
return _lgammaf_c0
+ _lgammaf_c1 / (x + 1)
+ _lgammaf_c2 / (x + 2)
+ _lgammaf_c3 / (x + 3)
+ _lgammaf_c4 / (x + 4)
+ _lgammaf_c5 / (x + 5)
+ _lgammaf_c6 / (x + 6);
}
double lgamma(double x)
{
if (isinf(x)) return INFINITY;
if (isnan(x)) return x;
if (x < 1)
{
if (floor(x) == x)
{ errno = ERANGE;
return HUGE_VAL;
}
return log(fabs((1-x)*_pi_/sin(_pi_*x)))-lgamma(2-x);
}
if (x == 1 || x == 2) return +0;
return (x+0.5)*log(x+18.5) - (x+18.5) + log(_gamma_sum(x) / x);
}
float lgammaf(float x)
{
if (isinf(x)) return INFINITY;
if (isnan(x)) return x;
if (x < 1)
{
if (floorf(x) == x)
{ errno = ERANGE;
return HUGE_VALF;
}
return (float) log(fabs((1.0-x)*_pi_/sin(_pi_*x)))-lgammaf(2-x);
}
if (x == 1 || x == 2) return +0;
return (float) ((x+0.5)*log(x+5.5) - (x+5.5) + log(_sqrt2pi * _gammaf_sum(x) / x));
}
double tgamma(double x)
{
if (isinf(x))
{ if (x > 0)
return x;
else
{ errno = ERANGE;
return NAN;
}
}
if (isnan(x)) return x;
if (x == 0)
{ errno = ERANGE;
return copysign(HUGE_VAL, x);
}
if (floor(x) == x)
{
if (x < 0)
{ errno = EDOM;
return NAN;
}
else if (x <= 171)
{ double r = 1;
for (int n = (int) x - 1; n > 1; n--)
r *= n;
return r;
}
}
if (x < 1)
return (1-x)*_pi_/(sin(_pi_*x)*tgamma(2-x));
return exp(lgamma(x));
}
float tgammaf(float x)
{
if (isinf(x))
{ if (x > 0)
return x;
else
{ errno = ERANGE;
return NAN;
}
}
if (isnan(x)) return x;
if (x == 0)
{ errno = ERANGE;
return copysignf(HUGE_VALF, x);
}
if (floorf(x) == x)
{
if (x < 0)
{ errno = EDOM;
return NAN;
}
else if (x <= 35)
{
float r = 1;
for (int n = (int) x - 1; n > 1; n--)
r *= n;
return r;
}
}
if (x < 1)
return (float) ((1.0-x)*_pi_/(sin(_pi_*x)*tgammaf(2-x)));
return expf(lgammaf(x));
}
/* Error function algorithms derived from "Numerical recipes in C" */
#define _log_sqrt_pi 0.5723649429247000870717137 // == lgamma(0.5)
/* Series calculation for incomplete gamma function P(0.5,x); good for */
/* x <= 1.5. */
static double gser05(double x, double epsilon)
{
double sum,del,ap;
ap = 0.5;
del = sum = 2;
for (;;)
{ ++ap;
del *= x/ap;
sum += del;
if (fabs(del) < fabs(sum)*epsilon)
return sum*exp(-x+0.5*log(x)-_log_sqrt_pi);
}
}
/* Continued fraction calculation for incomplete gamma function */
/* 1 - P(0.5,x); good for x >= 1.5. */
static double gcf05(double x, double epsilon)
{
double an,b,c,d,del,h;
#define FPMIN 1e-300
b = x+0.5;
c = 1/FPMIN;
d = 1/b;
h = d;
for (int i=1; ; i++)
{ an = i*(0.5-i);
b += 2;
d = an*d+b;
if (fabs(d) < FPMIN) d = FPMIN;
c = b + an/c;
if (fabs(c) < FPMIN) c = FPMIN;
d = 1/d;
del = d*c;
h *= del;
if (fabs(del-1) < epsilon) break;
}
return exp(-x+0.5*log(x)-_log_sqrt_pi)*h;
}
static double gammp05(double x, double epsilon)
{
if (isunordered(x, 1.5))
return x;
else if (isless(x, 1.5))
return gser05(x,epsilon);
else
return 1 - gcf05(x,epsilon);
}
static double gammq05(double x, double epsilon)
{
if (isunordered(x, 1.5))
return x;
else if (isless(x, 1.5))
return 1 - gser05(x,epsilon);
else
return gcf05(x,epsilon);
}
double erf(double x)
{
if (x == 0) return x;
if (isless(x, 0.0)) return -erf(-x);
if (isgreater(x, 1e100)) return 1.0;
return gammp05(x*x, 3*DBL_EPSILON);
}
double erfc(double x)
{
if (isgreater(fabs(x), 1e100)) return isless(x, 0.0) ? 2 : 0;
return isless(x, 0.0) ? 1+gammp05(x*x, 3*DBL_EPSILON)
: gammq05(x*x, 3*DBL_EPSILON);
}
float erff(float x)
{
if (x == 0) return x;
if (isless(x, 0.0)) return -erff(-x);
if (isgreater(x, 1e15F)) return 1.0F;
return (float) gammp05((double) x*x, 3.0*FLT_EPSILON);
}
#define _erfcf_c0 -1.26551223
#define _erfcf_c1 1.00002368
#define _erfcf_c2 0.37409196
#define _erfcf_c3 0.09678418
#define _erfcf_c4 -0.18628806
#define _erfcf_c5 0.27886807
#define _erfcf_c6 -1.13520398
#define _erfcf_c7 1.48851586
#define _erfcf_c8 -0.82215223
#define _erfcf_c9 0.17087277
float erfcf(float x)
{
double t,r;
if (isgreater(fabsf(x), 1e15F)) return isless(x, 0.0F) ? 2 : 0;
t = 1/(1+0.5*fabsf(x));
r = (((((((((_erfcf_c9) * t + _erfcf_c8) * t + _erfcf_c7) * t +
_erfcf_c6) * t + _erfcf_c5) * t + _erfcf_c4) * t +
_erfcf_c3) * t + _erfcf_c2) * t + _erfcf_c1) * t +
_erfcf_c0;
r = t*exp(r - (double)x*x);
return (float) (isless(x, 0.0F) ? 2-r : r);
}
#endif /* NO_FLOATING_POINT */
/* end of math.c */