/* 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 #include #include /* This file contains code for most of the math routines from */ /* 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 /* for forward references */ #include /* 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 */