Index: head/lib/msun/ld128/e_powl.c =================================================================== --- head/lib/msun/ld128/e_powl.c (revision 336361) +++ head/lib/msun/ld128/e_powl.c (revision 336362) @@ -1,443 +1,443 @@ /*- * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* * Copyright (c) 2008 Stephen L. Moshier * * Permission to use, copy, modify, and distribute this software for any * purpose with or without fee is hereby granted, provided that the above * copyright notice and this permission notice appear in all copies. * * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ /* powl(x,y) return x**y * * n * Method: Let x = 2 * (1+f) * 1. Compute and return log2(x) in two pieces: * log2(x) = w1 + w2, * where w1 has 113-53 = 60 bit trailing zeros. * 2. Perform y*log2(x) = n+y' by simulating muti-precision * arithmetic, where |y'|<=0.5. * 3. Return x**y = 2**n*exp(y'*log2) * * Special cases: * 1. (anything) ** 0 is 1 * 2. (anything) ** 1 is itself * 3. (anything) ** NAN is NAN * 4. NAN ** (anything except 0) is NAN * 5. +-(|x| > 1) ** +INF is +INF * 6. +-(|x| > 1) ** -INF is +0 * 7. +-(|x| < 1) ** +INF is +0 * 8. +-(|x| < 1) ** -INF is +INF * 9. +-1 ** +-INF is NAN * 10. +0 ** (+anything except 0, NAN) is +0 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 * 12. +0 ** (-anything except 0, NAN) is +INF * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF * 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) * 15. +INF ** (+anything except 0,NAN) is +INF * 16. +INF ** (-anything except 0,NAN) is +0 * 17. -INF ** (anything) = -0 ** (-anything) * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) * 19. (-anything except 0 and inf) ** (non-integer) is NAN * */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" static const long double bp[] = { 1.0L, 1.5L, }; /* log_2(1.5) */ static const long double dp_h[] = { 0.0, 5.8496250072115607565592654282227158546448E-1L }; /* Low part of log_2(1.5) */ static const long double dp_l[] = { 0.0, 1.0579781240112554492329533686862998106046E-16L }; static const long double zero = 0.0L, one = 1.0L, two = 2.0L, two113 = 1.0384593717069655257060992658440192E34L, huge = 1.0e3000L, tiny = 1.0e-3000L; /* 3/2 log x = 3 z + z^3 + z^3 (z^2 R(z^2)) z = (x-1)/(x+1) 1 <= x <= 1.25 Peak relative error 2.3e-37 */ static const long double LN[] = { -3.0779177200290054398792536829702930623200E1L, 6.5135778082209159921251824580292116201640E1L, -4.6312921812152436921591152809994014413540E1L, 1.2510208195629420304615674658258363295208E1L, -9.9266909031921425609179910128531667336670E-1L }; static const long double LD[] = { -5.129862866715009066465422805058933131960E1L, 1.452015077564081884387441590064272782044E2L, -1.524043275549860505277434040464085593165E2L, 7.236063513651544224319663428634139768808E1L, -1.494198912340228235853027849917095580053E1L /* 1.0E0 */ }; /* exp(x) = 1 + x - x / (1 - 2 / (x - x^2 R(x^2))) 0 <= x <= 0.5 Peak relative error 5.7e-38 */ static const long double PN[] = { 5.081801691915377692446852383385968225675E8L, 9.360895299872484512023336636427675327355E6L, 4.213701282274196030811629773097579432957E4L, 5.201006511142748908655720086041570288182E1L, 9.088368420359444263703202925095675982530E-3L, }; static const long double PD[] = { 3.049081015149226615468111430031590411682E9L, 1.069833887183886839966085436512368982758E8L, 8.259257717868875207333991924545445705394E5L, 1.872583833284143212651746812884298360922E3L, /* 1.0E0 */ }; static const long double /* ln 2 */ lg2 = 6.9314718055994530941723212145817656807550E-1L, lg2_h = 6.9314718055994528622676398299518041312695E-1L, lg2_l = 2.3190468138462996154948554638754786504121E-17L, ovt = 8.0085662595372944372e-0017L, /* 2/(3*log(2)) */ cp = 9.6179669392597560490661645400126142495110E-1L, cp_h = 9.6179669392597555432899980587535537779331E-1L, cp_l = 5.0577616648125906047157785230014751039424E-17L; long double powl(long double x, long double y) { long double z, ax, z_h, z_l, p_h, p_l; long double yy1, t1, t2, r, s, t, u, v, w; long double s2, s_h, s_l, t_h, t_l; int32_t i, j, k, yisint, n; u_int32_t ix, iy; int32_t hx, hy; ieee_quad_shape_type o, p, q; p.value = x; hx = p.parts32.mswhi; ix = hx & 0x7fffffff; q.value = y; hy = q.parts32.mswhi; iy = hy & 0x7fffffff; /* y==zero: x**0 = 1 */ if ((iy | q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0) return one; /* 1.0**y = 1; -1.0**+-Inf = 1 */ if (x == one) return one; if (x == -1.0L && iy == 0x7fff0000 && (q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0) return one; /* +-NaN return x+y */ if ((ix > 0x7fff0000) || ((ix == 0x7fff0000) && ((p.parts32.mswlo | p.parts32.lswhi | p.parts32.lswlo) != 0)) || (iy > 0x7fff0000) || ((iy == 0x7fff0000) && ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) != 0))) - return x + y; + return nan_mix(x, y); /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer * yisint = 1 ... y is an odd int * yisint = 2 ... y is an even int */ yisint = 0; if (hx < 0) { if (iy >= 0x40700000) /* 2^113 */ yisint = 2; /* even integer y */ else if (iy >= 0x3fff0000) /* 1.0 */ { if (floorl (y) == y) { z = 0.5 * y; if (floorl (z) == z) yisint = 2; else yisint = 1; } } } /* special value of y */ if ((q.parts32.mswlo | q.parts32.lswhi | q.parts32.lswlo) == 0) { if (iy == 0x7fff0000) /* y is +-inf */ { if (((ix - 0x3fff0000) | p.parts32.mswlo | p.parts32.lswhi | p.parts32.lswlo) == 0) return y - y; /* +-1**inf is NaN */ else if (ix >= 0x3fff0000) /* (|x|>1)**+-inf = inf,0 */ return (hy >= 0) ? y : zero; else /* (|x|<1)**-,+inf = inf,0 */ return (hy < 0) ? -y : zero; } if (iy == 0x3fff0000) { /* y is +-1 */ if (hy < 0) return one / x; else return x; } if (hy == 0x40000000) return x * x; /* y is 2 */ if (hy == 0x3ffe0000) { /* y is 0.5 */ if (hx >= 0) /* x >= +0 */ return sqrtl (x); } } ax = fabsl (x); /* special value of x */ if ((p.parts32.mswlo | p.parts32.lswhi | p.parts32.lswlo) == 0) { if (ix == 0x7fff0000 || ix == 0 || ix == 0x3fff0000) { z = ax; /*x is +-0,+-inf,+-1 */ if (hy < 0) z = one / z; /* z = (1/|x|) */ if (hx < 0) { if (((ix - 0x3fff0000) | yisint) == 0) { z = (z - z) / (z - z); /* (-1)**non-int is NaN */ } else if (yisint == 1) z = -z; /* (x<0)**odd = -(|x|**odd) */ } return z; } } /* (x<0)**(non-int) is NaN */ if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0) return (x - x) / (x - x); /* |y| is huge. 2^-16495 = 1/2 of smallest representable value. If (1 - 1/131072)^y underflows, y > 1.4986e9 */ if (iy > 0x401d654b) { /* if (1 - 2^-113)^y underflows, y > 1.1873e38 */ if (iy > 0x407d654b) { if (ix <= 0x3ffeffff) return (hy < 0) ? huge * huge : tiny * tiny; if (ix >= 0x3fff0000) return (hy > 0) ? huge * huge : tiny * tiny; } /* over/underflow if x is not close to one */ if (ix < 0x3ffeffff) return (hy < 0) ? huge * huge : tiny * tiny; if (ix > 0x3fff0000) return (hy > 0) ? huge * huge : tiny * tiny; } n = 0; /* take care subnormal number */ if (ix < 0x00010000) { ax *= two113; n -= 113; o.value = ax; ix = o.parts32.mswhi; } n += ((ix) >> 16) - 0x3fff; j = ix & 0x0000ffff; /* determine interval */ ix = j | 0x3fff0000; /* normalize ix */ if (j <= 0x3988) k = 0; /* |x|> 31) - 1) | (yisint - 1)) == 0) s = -one; /* (-ve)**(odd int) */ /* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */ yy1 = y; o.value = yy1; o.parts32.lswlo = 0; o.parts32.lswhi &= 0xf8000000; yy1 = o.value; p_l = (y - yy1) * t1 + y * t2; p_h = yy1 * t1; z = p_l + p_h; o.value = z; j = o.parts32.mswhi; if (j >= 0x400d0000) /* z >= 16384 */ { /* if z > 16384 */ if (((j - 0x400d0000) | o.parts32.mswlo | o.parts32.lswhi | o.parts32.lswlo) != 0) return s * huge * huge; /* overflow */ else { if (p_l + ovt > z - p_h) return s * huge * huge; /* overflow */ } } else if ((j & 0x7fffffff) >= 0x400d01b9) /* z <= -16495 */ { /* z < -16495 */ if (((j - 0xc00d01bc) | o.parts32.mswlo | o.parts32.lswhi | o.parts32.lswlo) != 0) return s * tiny * tiny; /* underflow */ else { if (p_l <= z - p_h) return s * tiny * tiny; /* underflow */ } } /* compute 2**(p_h+p_l) */ i = j & 0x7fffffff; k = (i >> 16) - 0x3fff; n = 0; if (i > 0x3ffe0000) { /* if |z| > 0.5, set n = [z+0.5] */ n = floorl (z + 0.5L); t = n; p_h -= t; } t = p_l + p_h; o.value = t; o.parts32.lswlo = 0; o.parts32.lswhi &= 0xf8000000; t = o.value; u = t * lg2_h; v = (p_l - (t - p_h)) * lg2 + t * lg2_l; z = u + v; w = v - (z - u); /* exp(z) */ t = z * z; u = PN[0] + t * (PN[1] + t * (PN[2] + t * (PN[3] + t * PN[4]))); v = PD[0] + t * (PD[1] + t * (PD[2] + t * (PD[3] + t))); t1 = z - t * u / v; r = (z * t1) / (t1 - two) - (w + z * w); z = one - (r - z); o.value = z; j = o.parts32.mswhi; j += (n << 16); if ((j >> 16) <= 0) z = scalbnl (z, n); /* subnormal output */ else { o.parts32.mswhi = j; z = o.value; } return s * z; } Index: head/lib/msun/ld80/e_powl.c =================================================================== --- head/lib/msun/ld80/e_powl.c (revision 336361) +++ head/lib/msun/ld80/e_powl.c (revision 336362) @@ -1,616 +1,616 @@ /*- * Copyright (c) 2008 Stephen L. Moshier * * Permission to use, copy, modify, and distribute this software for any * purpose with or without fee is hereby granted, provided that the above * copyright notice and this permission notice appear in all copies. * * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ /* powl.c * * Power function, long double precision * * * * SYNOPSIS: * * long double x, y, z, powl(); * * z = powl( x, y ); * * * * DESCRIPTION: * * Computes x raised to the yth power. Analytically, * * x**y = exp( y log(x) ). * * Following Cody and Waite, this program uses a lookup table * of 2**-i/32 and pseudo extended precision arithmetic to * obtain several extra bits of accuracy in both the logarithm * and the exponential. * * * * ACCURACY: * * The relative error of pow(x,y) can be estimated * by y dl ln(2), where dl is the absolute error of * the internally computed base 2 logarithm. At the ends * of the approximation interval the logarithm equal 1/32 * and its relative error is about 1 lsb = 1.1e-19. Hence * the predicted relative error in the result is 2.3e-21 y . * * Relative error: * arithmetic domain # trials peak rms * * IEEE +-1000 40000 2.8e-18 3.7e-19 * .001 < x < 1000, with log(x) uniformly distributed. * -1000 < y < 1000, y uniformly distributed. * * IEEE 0,8700 60000 6.5e-18 1.0e-18 * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed. * * * ERROR MESSAGES: * * message condition value returned * pow overflow x**y > MAXNUM INFINITY * pow underflow x**y < 1/MAXNUM 0.0 * pow domain x<0 and y noninteger 0.0 * */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" /* Table size */ #define NXT 32 /* log2(Table size) */ #define LNXT 5 /* log(1+x) = x - .5x^2 + x^3 * P(z)/Q(z) * on the domain 2^(-1/32) - 1 <= x <= 2^(1/32) - 1 */ static long double P[] = { 8.3319510773868690346226E-4L, 4.9000050881978028599627E-1L, 1.7500123722550302671919E0L, 1.4000100839971580279335E0L, }; static long double Q[] = { /* 1.0000000000000000000000E0L,*/ 5.2500282295834889175431E0L, 8.4000598057587009834666E0L, 4.2000302519914740834728E0L, }; /* A[i] = 2^(-i/32), rounded to IEEE long double precision. * If i is even, A[i] + B[i/2] gives additional accuracy. */ static long double A[33] = { 1.0000000000000000000000E0L, 9.7857206208770013448287E-1L, 9.5760328069857364691013E-1L, 9.3708381705514995065011E-1L, 9.1700404320467123175367E-1L, 8.9735453750155359320742E-1L, 8.7812608018664974155474E-1L, 8.5930964906123895780165E-1L, 8.4089641525371454301892E-1L, 8.2287773907698242225554E-1L, 8.0524516597462715409607E-1L, 7.8799042255394324325455E-1L, 7.7110541270397041179298E-1L, 7.5458221379671136985669E-1L, 7.3841307296974965571198E-1L, 7.2259040348852331001267E-1L, 7.0710678118654752438189E-1L, 6.9195494098191597746178E-1L, 6.7712777346844636413344E-1L, 6.6261832157987064729696E-1L, 6.4841977732550483296079E-1L, 6.3452547859586661129850E-1L, 6.2092890603674202431705E-1L, 6.0762367999023443907803E-1L, 5.9460355750136053334378E-1L, 5.8186242938878875689693E-1L, 5.6939431737834582684856E-1L, 5.5719337129794626814472E-1L, 5.4525386633262882960438E-1L, 5.3357020033841180906486E-1L, 5.2213689121370692017331E-1L, 5.1094857432705833910408E-1L, 5.0000000000000000000000E-1L, }; static long double B[17] = { 0.0000000000000000000000E0L, 2.6176170809902549338711E-20L, -1.0126791927256478897086E-20L, 1.3438228172316276937655E-21L, 1.2207982955417546912101E-20L, -6.3084814358060867200133E-21L, 1.3164426894366316434230E-20L, -1.8527916071632873716786E-20L, 1.8950325588932570796551E-20L, 1.5564775779538780478155E-20L, 6.0859793637556860974380E-21L, -2.0208749253662532228949E-20L, 1.4966292219224761844552E-20L, 3.3540909728056476875639E-21L, -8.6987564101742849540743E-22L, -1.2327176863327626135542E-20L, 0.0000000000000000000000E0L, }; /* 2^x = 1 + x P(x), * on the interval -1/32 <= x <= 0 */ static long double R[] = { 1.5089970579127659901157E-5L, 1.5402715328927013076125E-4L, 1.3333556028915671091390E-3L, 9.6181291046036762031786E-3L, 5.5504108664798463044015E-2L, 2.4022650695910062854352E-1L, 6.9314718055994530931447E-1L, }; #define douba(k) A[k] #define doubb(k) B[k] #define MEXP (NXT*16384.0L) /* The following if denormal numbers are supported, else -MEXP: */ #define MNEXP (-NXT*(16384.0L+64.0L)) /* log2(e) - 1 */ #define LOG2EA 0.44269504088896340735992L #define F W #define Fa Wa #define Fb Wb #define G W #define Ga Wa #define Gb u #define H W #define Ha Wb #define Hb Wb static const long double MAXLOGL = 1.1356523406294143949492E4L; static const long double MINLOGL = -1.13994985314888605586758E4L; static const long double LOGE2L = 6.9314718055994530941723E-1L; static volatile long double z; static long double w, W, Wa, Wb, ya, yb, u; static const long double huge = 0x1p10000L; #if 0 /* XXX Prevent gcc from erroneously constant folding this. */ static const long double twom10000 = 0x1p-10000L; #else static volatile long double twom10000 = 0x1p-10000L; #endif static long double reducl( long double ); static long double powil ( long double, int ); long double powl(long double x, long double y) { /* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */ int i, nflg, iyflg, yoddint; long e; if( y == 0.0L ) return( 1.0L ); if( x == 1.0L ) return( 1.0L ); if( isnan(x) ) - return( x ); + return ( nan_mix(x, y) ); if( isnan(y) ) - return( y ); + return ( nan_mix(x, y) ); if( y == 1.0L ) return( x ); if( !isfinite(y) && x == -1.0L ) return( 1.0L ); if( y >= LDBL_MAX ) { if( x > 1.0L ) return( INFINITY ); if( x > 0.0L && x < 1.0L ) return( 0.0L ); if( x < -1.0L ) return( INFINITY ); if( x > -1.0L && x < 0.0L ) return( 0.0L ); } if( y <= -LDBL_MAX ) { if( x > 1.0L ) return( 0.0L ); if( x > 0.0L && x < 1.0L ) return( INFINITY ); if( x < -1.0L ) return( 0.0L ); if( x > -1.0L && x < 0.0L ) return( INFINITY ); } if( x >= LDBL_MAX ) { if( y > 0.0L ) return( INFINITY ); return( 0.0L ); } w = floorl(y); /* Set iyflg to 1 if y is an integer. */ iyflg = 0; if( w == y ) iyflg = 1; /* Test for odd integer y. */ yoddint = 0; if( iyflg ) { ya = fabsl(y); ya = floorl(0.5L * ya); yb = 0.5L * fabsl(w); if( ya != yb ) yoddint = 1; } if( x <= -LDBL_MAX ) { if( y > 0.0L ) { if( yoddint ) return( -INFINITY ); return( INFINITY ); } if( y < 0.0L ) { if( yoddint ) return( -0.0L ); return( 0.0 ); } } nflg = 0; /* flag = 1 if x<0 raised to integer power */ if( x <= 0.0L ) { if( x == 0.0L ) { if( y < 0.0 ) { if( signbit(x) && yoddint ) return( -INFINITY ); return( INFINITY ); } if( y > 0.0 ) { if( signbit(x) && yoddint ) return( -0.0L ); return( 0.0 ); } if( y == 0.0L ) return( 1.0L ); /* 0**0 */ else return( 0.0L ); /* 0**y */ } else { if( iyflg == 0 ) return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */ nflg = 1; } } /* Integer power of an integer. */ if( iyflg ) { i = w; w = floorl(x); if( (w == x) && (fabsl(y) < 32768.0) ) { w = powil( x, (int) y ); return( w ); } } if( nflg ) x = fabsl(x); /* separate significand from exponent */ x = frexpl( x, &i ); e = i; /* find significand in antilog table A[] */ i = 1; if( x <= douba(17) ) i = 17; if( x <= douba(i+8) ) i += 8; if( x <= douba(i+4) ) i += 4; if( x <= douba(i+2) ) i += 2; if( x >= douba(1) ) i = -1; i += 1; /* Find (x - A[i])/A[i] * in order to compute log(x/A[i]): * * log(x) = log( a x/a ) = log(a) + log(x/a) * * log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a */ x -= douba(i); x -= doubb(i/2); x /= douba(i); /* rational approximation for log(1+v): * * log(1+v) = v - v**2/2 + v**3 P(v) / Q(v) */ z = x*x; w = x * ( z * __polevll( x, P, 3 ) / __p1evll( x, Q, 3 ) ); w = w - ldexpl( z, -1 ); /* w - 0.5 * z */ /* Convert to base 2 logarithm: * multiply by log2(e) = 1 + LOG2EA */ z = LOG2EA * w; z += w; z += LOG2EA * x; z += x; /* Compute exponent term of the base 2 logarithm. */ w = -i; w = ldexpl( w, -LNXT ); /* divide by NXT */ w += e; /* Now base 2 log of x is w + z. */ /* Multiply base 2 log by y, in extended precision. */ /* separate y into large part ya * and small part yb less than 1/NXT */ ya = reducl(y); yb = y - ya; /* (w+z)(ya+yb) * = w*ya + w*yb + z*y */ F = z * y + w * yb; Fa = reducl(F); Fb = F - Fa; G = Fa + w * ya; Ga = reducl(G); Gb = G - Ga; H = Fb + Gb; Ha = reducl(H); w = ldexpl( Ga+Ha, LNXT ); /* Test the power of 2 for overflow */ if( w > MEXP ) return (huge * huge); /* overflow */ if( w < MNEXP ) return (twom10000 * twom10000); /* underflow */ e = w; Hb = H - Ha; if( Hb > 0.0L ) { e += 1; Hb -= (1.0L/NXT); /*0.0625L;*/ } /* Now the product y * log2(x) = Hb + e/NXT. * * Compute base 2 exponential of Hb, * where -0.0625 <= Hb <= 0. */ z = Hb * __polevll( Hb, R, 6 ); /* z = 2**Hb - 1 */ /* Express e/NXT as an integer plus a negative number of (1/NXT)ths. * Find lookup table entry for the fractional power of 2. */ if( e < 0 ) i = 0; else i = 1; i = e/NXT + i; e = NXT*i - e; w = douba( e ); z = w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */ z = z + w; z = ldexpl( z, i ); /* multiply by integer power of 2 */ if( nflg ) { /* For negative x, * find out if the integer exponent * is odd or even. */ w = ldexpl( y, -1 ); w = floorl(w); w = ldexpl( w, 1 ); if( w != y ) z = -z; /* odd exponent */ } return( z ); } /* Find a multiple of 1/NXT that is within 1/NXT of x. */ static long double reducl(long double x) { long double t; t = ldexpl( x, LNXT ); t = floorl( t ); t = ldexpl( t, -LNXT ); return(t); } /* powil.c * * Real raised to integer power, long double precision * * * * SYNOPSIS: * * long double x, y, powil(); * int n; * * y = powil( x, n ); * * * * DESCRIPTION: * * Returns argument x raised to the nth power. * The routine efficiently decomposes n as a sum of powers of * two. The desired power is a product of two-to-the-kth * powers of x. Thus to compute the 32767 power of x requires * 28 multiplications instead of 32767 multiplications. * * * * ACCURACY: * * * Relative error: * arithmetic x domain n domain # trials peak rms * IEEE .001,1000 -1022,1023 50000 4.3e-17 7.8e-18 * IEEE 1,2 -1022,1023 20000 3.9e-17 7.6e-18 * IEEE .99,1.01 0,8700 10000 3.6e-16 7.2e-17 * * Returns MAXNUM on overflow, zero on underflow. * */ static long double powil(long double x, int nn) { long double ww, y; long double s; int n, e, sign, asign, lx; if( x == 0.0L ) { if( nn == 0 ) return( 1.0L ); else if( nn < 0 ) return( LDBL_MAX ); else return( 0.0L ); } if( nn == 0 ) return( 1.0L ); if( x < 0.0L ) { asign = -1; x = -x; } else asign = 0; if( nn < 0 ) { sign = -1; n = -nn; } else { sign = 1; n = nn; } /* Overflow detection */ /* Calculate approximate logarithm of answer */ s = x; s = frexpl( s, &lx ); e = (lx - 1)*n; if( (e == 0) || (e > 64) || (e < -64) ) { s = (s - 7.0710678118654752e-1L) / (s + 7.0710678118654752e-1L); s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L; } else { s = LOGE2L * e; } if( s > MAXLOGL ) return (huge * huge); /* overflow */ if( s < MINLOGL ) return (twom10000 * twom10000); /* underflow */ /* Handle tiny denormal answer, but with less accuracy * since roundoff error in 1.0/x will be amplified. * The precise demarcation should be the gradual underflow threshold. */ if( s < (-MAXLOGL+2.0L) ) { x = 1.0L/x; sign = -sign; } /* First bit of the power */ if( n & 1 ) y = x; else { y = 1.0L; asign = 0; } ww = x; n >>= 1; while( n ) { ww = ww * ww; /* arg to the 2-to-the-kth power */ if( n & 1 ) /* if that bit is set, then include in product */ y *= ww; n >>= 1; } if( asign ) y = -y; /* odd power of negative number */ if( sign < 0 ) y = 1.0L/y; return(y); } Index: head/lib/msun/src/catrig.c =================================================================== --- head/lib/msun/src/catrig.c (revision 336361) +++ head/lib/msun/src/catrig.c (revision 336362) @@ -1,655 +1,655 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2012 Stephen Montgomery-Smith * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ #include __FBSDID("$FreeBSD$"); #include #include #include "math.h" #include "math_private.h" #undef isinf #define isinf(x) (fabs(x) == INFINITY) #undef isnan #define isnan(x) ((x) != (x)) #define raise_inexact() do { volatile float junk __unused = 1 + tiny; } while(0) #undef signbit #define signbit(x) (__builtin_signbit(x)) /* We need that DBL_EPSILON^2/128 is larger than FOUR_SQRT_MIN. */ static const double A_crossover = 10, /* Hull et al suggest 1.5, but 10 works better */ B_crossover = 0.6417, /* suggested by Hull et al */ FOUR_SQRT_MIN = 0x1p-509, /* >= 4 * sqrt(DBL_MIN) */ QUARTER_SQRT_MAX = 0x1p509, /* <= sqrt(DBL_MAX) / 4 */ m_e = 2.7182818284590452e0, /* 0x15bf0a8b145769.0p-51 */ m_ln2 = 6.9314718055994531e-1, /* 0x162e42fefa39ef.0p-53 */ pio2_hi = 1.5707963267948966e0, /* 0x1921fb54442d18.0p-52 */ RECIP_EPSILON = 1 / DBL_EPSILON, SQRT_3_EPSILON = 2.5809568279517849e-8, /* 0x1bb67ae8584caa.0p-78 */ SQRT_6_EPSILON = 3.6500241499888571e-8, /* 0x13988e1409212e.0p-77 */ SQRT_MIN = 0x1p-511; /* >= sqrt(DBL_MIN) */ static const volatile double pio2_lo = 6.1232339957367659e-17; /* 0x11a62633145c07.0p-106 */ static const volatile float tiny = 0x1p-100; static double complex clog_for_large_values(double complex z); /* * Testing indicates that all these functions are accurate up to 4 ULP. * The functions casin(h) and cacos(h) are about 2.5 times slower than asinh. * The functions catan(h) are a little under 2 times slower than atanh. * * The code for casinh, casin, cacos, and cacosh comes first. The code is * rather complicated, and the four functions are highly interdependent. * * The code for catanh and catan comes at the end. It is much simpler than * the other functions, and the code for these can be disconnected from the * rest of the code. */ /* * ================================ * | casinh, casin, cacos, cacosh | * ================================ */ /* * The algorithm is very close to that in "Implementing the complex arcsine * and arccosine functions using exception handling" by T. E. Hull, Thomas F. * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335, * http://dl.acm.org/citation.cfm?id=275324. * * Throughout we use the convention z = x + I*y. * * casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B) * where * A = (|z+I| + |z-I|) / 2 * B = (|z+I| - |z-I|) / 2 = y/A * * These formulas become numerically unstable: * (a) for Re(casinh(z)) when z is close to the line segment [-I, I] (that * is, Re(casinh(z)) is close to 0); * (b) for Im(casinh(z)) when z is close to either of the intervals * [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is * close to PI/2). * * These numerical problems are overcome by defining * f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2 * Then if A < A_crossover, we use * log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1))) * A-1 = f(x, 1+y) + f(x, 1-y) * and if B > B_crossover, we use * asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y))) * A-y = f(x, y+1) + f(x, y-1) * where without loss of generality we have assumed that x and y are * non-negative. * * Much of the difficulty comes because the intermediate computations may * produce overflows or underflows. This is dealt with in the paper by Hull * et al by using exception handling. We do this by detecting when * computations risk underflow or overflow. The hardest part is handling the * underflows when computing f(a, b). * * Note that the function f(a, b) does not appear explicitly in the paper by * Hull et al, but the idea may be found on pages 308 and 309. Introducing the * function f(a, b) allows us to concentrate many of the clever tricks in this * paper into one function. */ /* * Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2. * Pass hypot(a, b) as the third argument. */ static inline double f(double a, double b, double hypot_a_b) { if (b < 0) return ((hypot_a_b - b) / 2); if (b == 0) return (a / 2); return (a * a / (hypot_a_b + b) / 2); } /* * All the hard work is contained in this function. * x and y are assumed positive or zero, and less than RECIP_EPSILON. * Upon return: * rx = Re(casinh(z)) = -Im(cacos(y + I*x)). * B_is_usable is set to 1 if the value of B is usable. * If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y. * If returning sqrt_A2my2 has potential to result in an underflow, it is * rescaled, and new_y is similarly rescaled. */ static inline void do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B, double *sqrt_A2my2, double *new_y) { double R, S, A; /* A, B, R, and S are as in Hull et al. */ double Am1, Amy; /* A-1, A-y. */ R = hypot(x, y + 1); /* |z+I| */ S = hypot(x, y - 1); /* |z-I| */ /* A = (|z+I| + |z-I|) / 2 */ A = (R + S) / 2; /* * Mathematically A >= 1. There is a small chance that this will not * be so because of rounding errors. So we will make certain it is * so. */ if (A < 1) A = 1; if (A < A_crossover) { /* * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y). * rx = log1p(Am1 + sqrt(Am1*(A+1))) */ if (y == 1 && x < DBL_EPSILON * DBL_EPSILON / 128) { /* * fp is of order x^2, and fm = x/2. * A = 1 (inexactly). */ *rx = sqrt(x); } else if (x >= DBL_EPSILON * fabs(y - 1)) { /* * Underflow will not occur because * x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN */ Am1 = f(x, 1 + y, R) + f(x, 1 - y, S); *rx = log1p(Am1 + sqrt(Am1 * (A + 1))); } else if (y < 1) { /* * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and * A = 1 (inexactly). */ *rx = x / sqrt((1 - y) * (1 + y)); } else { /* if (y > 1) */ /* * A-1 = y-1 (inexactly). */ *rx = log1p((y - 1) + sqrt((y - 1) * (y + 1))); } } else { *rx = log(A + sqrt(A * A - 1)); } *new_y = y; if (y < FOUR_SQRT_MIN) { /* * Avoid a possible underflow caused by y/A. For casinh this * would be legitimate, but will be picked up by invoking atan2 * later on. For cacos this would not be legitimate. */ *B_is_usable = 0; *sqrt_A2my2 = A * (2 / DBL_EPSILON); *new_y = y * (2 / DBL_EPSILON); return; } /* B = (|z+I| - |z-I|) / 2 = y/A */ *B = y / A; *B_is_usable = 1; if (*B > B_crossover) { *B_is_usable = 0; /* * Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1). * sqrt_A2my2 = sqrt(Amy*(A+y)) */ if (y == 1 && x < DBL_EPSILON / 128) { /* * fp is of order x^2, and fm = x/2. * A = 1 (inexactly). */ *sqrt_A2my2 = sqrt(x) * sqrt((A + y) / 2); } else if (x >= DBL_EPSILON * fabs(y - 1)) { /* * Underflow will not occur because * x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN * and * x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN */ Amy = f(x, y + 1, R) + f(x, y - 1, S); *sqrt_A2my2 = sqrt(Amy * (A + y)); } else if (y > 1) { /* * fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and * A = y (inexactly). * * y < RECIP_EPSILON. So the following * scaling should avoid any underflow problems. */ *sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y / sqrt((y + 1) * (y - 1)); *new_y = y * (4 / DBL_EPSILON / DBL_EPSILON); } else { /* if (y < 1) */ /* * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and * A = 1 (inexactly). */ *sqrt_A2my2 = sqrt((1 - y) * (1 + y)); } } } /* * casinh(z) = z + O(z^3) as z -> 0 * * casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2) as z -> infinity * The above formula works for the imaginary part as well, because * Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3) * as z -> infinity, uniformly in y */ double complex casinh(double complex z) { double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; int B_is_usable; double complex w; x = creal(z); y = cimag(z); ax = fabs(x); ay = fabs(y); if (isnan(x) || isnan(y)) { /* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */ if (isinf(x)) return (CMPLX(x, y + y)); /* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */ if (isinf(y)) return (CMPLX(y, x + x)); /* casinh(NaN + I*0) = NaN + I*0 */ if (y == 0) return (CMPLX(x + x, y)); /* * All other cases involving NaN return NaN + I*NaN. * C99 leaves it optional whether to raise invalid if one of * the arguments is not NaN, so we opt not to raise it. */ - return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLX(nan_mix(x, y), nan_mix(x, y))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { /* clog...() will raise inexact unless x or y is infinite. */ if (signbit(x) == 0) w = clog_for_large_values(z) + m_ln2; else w = clog_for_large_values(-z) + m_ln2; return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y))); } /* Avoid spuriously raising inexact for z = 0. */ if (x == 0 && y == 0) return (z); /* All remaining cases are inexact. */ raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) return (z); do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); if (B_is_usable) ry = asin(B); else ry = atan2(new_y, sqrt_A2my2); return (CMPLX(copysign(rx, x), copysign(ry, y))); } /* * casin(z) = reverse(casinh(reverse(z))) * where reverse(x + I*y) = y + I*x = I*conj(z). */ double complex casin(double complex z) { double complex w = casinh(CMPLX(cimag(z), creal(z))); return (CMPLX(cimag(w), creal(w))); } /* * cacos(z) = PI/2 - casin(z) * but do the computation carefully so cacos(z) is accurate when z is * close to 1. * * cacos(z) = PI/2 - z + O(z^3) as z -> 0 * * cacos(z) = -sign(y)*I*clog(z) + O(1/z^2) as z -> infinity * The above formula works for the real part as well, because * Re(cacos(z)) = atan2(fabs(y), x) + O(y/z^3) * as z -> infinity, uniformly in y */ double complex cacos(double complex z) { double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; int sx, sy; int B_is_usable; double complex w; x = creal(z); y = cimag(z); sx = signbit(x); sy = signbit(y); ax = fabs(x); ay = fabs(y); if (isnan(x) || isnan(y)) { /* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */ if (isinf(x)) return (CMPLX(y + y, -INFINITY)); /* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */ if (isinf(y)) return (CMPLX(x + x, -y)); /* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */ if (x == 0) return (CMPLX(pio2_hi + pio2_lo, y + y)); /* * All other cases involving NaN return NaN + I*NaN. * C99 leaves it optional whether to raise invalid if one of * the arguments is not NaN, so we opt not to raise it. */ - return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLX(nan_mix(x, y), nan_mix(x, y))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { /* clog...() will raise inexact unless x or y is infinite. */ w = clog_for_large_values(z); rx = fabs(cimag(w)); ry = creal(w) + m_ln2; if (sy == 0) ry = -ry; return (CMPLX(rx, ry)); } /* Avoid spuriously raising inexact for z = 1. */ if (x == 1 && y == 0) return (CMPLX(0, -y)); /* All remaining cases are inexact. */ raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) return (CMPLX(pio2_hi - (x - pio2_lo), -y)); do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); if (B_is_usable) { if (sx == 0) rx = acos(B); else rx = acos(-B); } else { if (sx == 0) rx = atan2(sqrt_A2mx2, new_x); else rx = atan2(sqrt_A2mx2, -new_x); } if (sy == 0) ry = -ry; return (CMPLX(rx, ry)); } /* * cacosh(z) = I*cacos(z) or -I*cacos(z) * where the sign is chosen so Re(cacosh(z)) >= 0. */ double complex cacosh(double complex z) { double complex w; double rx, ry; w = cacos(z); rx = creal(w); ry = cimag(w); /* cacosh(NaN + I*NaN) = NaN + I*NaN */ if (isnan(rx) && isnan(ry)) return (CMPLX(ry, rx)); /* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */ /* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */ if (isnan(rx)) return (CMPLX(fabs(ry), rx)); /* cacosh(0 + I*NaN) = NaN + I*NaN */ if (isnan(ry)) return (CMPLX(ry, ry)); return (CMPLX(fabs(ry), copysign(rx, cimag(z)))); } /* * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON. */ static double complex clog_for_large_values(double complex z) { double x, y; double ax, ay, t; x = creal(z); y = cimag(z); ax = fabs(x); ay = fabs(y); if (ax < ay) { t = ax; ax = ay; ay = t; } /* * Avoid overflow in hypot() when x and y are both very large. * Divide x and y by E, and then add 1 to the logarithm. This * depends on E being larger than sqrt(2), since the return value of * hypot cannot overflow if neither argument is greater in magnitude * than 1/sqrt(2) of the maximum value of the return type. Likewise * this determines the necessary threshold for using this method * (however, actually use 1/2 instead as it is simpler). * * Dividing by E causes an insignificant loss of accuracy; however * this method is still poor since it is uneccessarily slow. */ if (ax > DBL_MAX / 2) return (CMPLX(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x))); /* * Avoid overflow when x or y is large. Avoid underflow when x or * y is small. */ if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) return (CMPLX(log(hypot(x, y)), atan2(y, x))); return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x))); } /* * ================= * | catanh, catan | * ================= */ /* * sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow). * Assumes x*x and y*y will not overflow. * Assumes x and y are finite. * Assumes y is non-negative. * Assumes fabs(x) >= DBL_EPSILON. */ static inline double sum_squares(double x, double y) { /* Avoid underflow when y is small. */ if (y < SQRT_MIN) return (x * x); return (x * x + y * y); } /* * real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y). * Assumes x and y are not NaN, and one of x and y is larger than * RECIP_EPSILON. We avoid unwarranted underflow. It is important to not use * the code creal(1/z), because the imaginary part may produce an unwanted * underflow. * This is only called in a context where inexact is always raised before * the call, so no effort is made to avoid or force inexact. */ static inline double real_part_reciprocal(double x, double y) { double scale; uint32_t hx, hy; int32_t ix, iy; /* * This code is inspired by the C99 document n1124.pdf, Section G.5.1, * example 2. */ GET_HIGH_WORD(hx, x); ix = hx & 0x7ff00000; GET_HIGH_WORD(hy, y); iy = hy & 0x7ff00000; #define BIAS (DBL_MAX_EXP - 1) /* XXX more guard digits are useful iff there is extra precision. */ #define CUTOFF (DBL_MANT_DIG / 2 + 1) /* just half or 1 guard digit */ if (ix - iy >= CUTOFF << 20 || isinf(x)) return (1 / x); /* +-Inf -> +-0 is special */ if (iy - ix >= CUTOFF << 20) return (x / y / y); /* should avoid double div, but hard */ if (ix <= (BIAS + DBL_MAX_EXP / 2 - CUTOFF) << 20) return (x / (x * x + y * y)); scale = 1; SET_HIGH_WORD(scale, 0x7ff00000 - ix); /* 2**(1-ilogb(x)) */ x *= scale; y *= scale; return (x / (x * x + y * y) * scale); } /* * catanh(z) = log((1+z)/(1-z)) / 2 * = log1p(4*x / |z-1|^2) / 4 * + I * atan2(2*y, (1-x)*(1+x)-y*y) / 2 * * catanh(z) = z + O(z^3) as z -> 0 * * catanh(z) = 1/z + sign(y)*I*PI/2 + O(1/z^3) as z -> infinity * The above formula works for the real part as well, because * Re(catanh(z)) = x/|z|^2 + O(x/z^4) * as z -> infinity, uniformly in x */ double complex catanh(double complex z) { double x, y, ax, ay, rx, ry; x = creal(z); y = cimag(z); ax = fabs(x); ay = fabs(y); /* This helps handle many cases. */ if (y == 0 && ax <= 1) return (CMPLX(atanh(x), y)); /* To ensure the same accuracy as atan(), and to filter out z = 0. */ if (x == 0) return (CMPLX(x, atan(y))); if (isnan(x) || isnan(y)) { /* catanh(+-Inf + I*NaN) = +-0 + I*NaN */ if (isinf(x)) return (CMPLX(copysign(0, x), y + y)); /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */ if (isinf(y)) return (CMPLX(copysign(0, x), copysign(pio2_hi + pio2_lo, y))); /* * All other cases involving NaN return NaN + I*NaN. * C99 leaves it optional whether to raise invalid if one of * the arguments is not NaN, so we opt not to raise it. */ - return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLX(nan_mix(x, y), nan_mix(x, y))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) return (CMPLX(real_part_reciprocal(x, y), copysign(pio2_hi + pio2_lo, y))); if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { /* * z = 0 was filtered out above. All other cases must raise * inexact, but this is the only case that needs to do it * explicitly. */ raise_inexact(); return (z); } if (ax == 1 && ay < DBL_EPSILON) rx = (m_ln2 - log(ay)) / 2; else rx = log1p(4 * ax / sum_squares(ax - 1, ay)) / 4; if (ax == 1) ry = atan2(2, -ay) / 2; else if (ay < DBL_EPSILON) ry = atan2(2 * ay, (1 - ax) * (1 + ax)) / 2; else ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; return (CMPLX(copysign(rx, x), copysign(ry, y))); } /* * catan(z) = reverse(catanh(reverse(z))) * where reverse(x + I*y) = y + I*x = I*conj(z). */ double complex catan(double complex z) { double complex w = catanh(CMPLX(cimag(z), creal(z))); return (CMPLX(cimag(w), creal(w))); } #if LDBL_MANT_DIG == 53 __weak_reference(cacosh, cacoshl); __weak_reference(cacos, cacosl); __weak_reference(casinh, casinhl); __weak_reference(casin, casinl); __weak_reference(catanh, catanhl); __weak_reference(catan, catanl); #endif Index: head/lib/msun/src/catrigf.c =================================================================== --- head/lib/msun/src/catrigf.c (revision 336361) +++ head/lib/msun/src/catrigf.c (revision 336362) @@ -1,395 +1,395 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2012 Stephen Montgomery-Smith * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ /* * The algorithm is very close to that in "Implementing the complex arcsine * and arccosine functions using exception handling" by T. E. Hull, Thomas F. * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335, * http://dl.acm.org/citation.cfm?id=275324. * * See catrig.c for complete comments. * * XXX comments were removed automatically, and even short ones on the right * of statements were removed (all of them), contrary to normal style. Only * a few comments on the right of declarations remain. */ #include __FBSDID("$FreeBSD$"); #include #include #include "math.h" #include "math_private.h" #undef isinf #define isinf(x) (fabsf(x) == INFINITY) #undef isnan #define isnan(x) ((x) != (x)) #define raise_inexact() do { volatile float junk __unused = 1 + tiny; } while(0) #undef signbit #define signbit(x) (__builtin_signbitf(x)) static const float A_crossover = 10, B_crossover = 0.6417, FOUR_SQRT_MIN = 0x1p-61, QUARTER_SQRT_MAX = 0x1p61, m_e = 2.7182818285e0, /* 0xadf854.0p-22 */ m_ln2 = 6.9314718056e-1, /* 0xb17218.0p-24 */ pio2_hi = 1.5707962513e0, /* 0xc90fda.0p-23 */ RECIP_EPSILON = 1 / FLT_EPSILON, SQRT_3_EPSILON = 5.9801995673e-4, /* 0x9cc471.0p-34 */ SQRT_6_EPSILON = 8.4572793338e-4, /* 0xddb3d7.0p-34 */ SQRT_MIN = 0x1p-63; static const volatile float pio2_lo = 7.5497899549e-8, /* 0xa22169.0p-47 */ tiny = 0x1p-100; static float complex clog_for_large_values(float complex z); static inline float f(float a, float b, float hypot_a_b) { if (b < 0) return ((hypot_a_b - b) / 2); if (b == 0) return (a / 2); return (a * a / (hypot_a_b + b) / 2); } static inline void do_hard_work(float x, float y, float *rx, int *B_is_usable, float *B, float *sqrt_A2my2, float *new_y) { float R, S, A; float Am1, Amy; R = hypotf(x, y + 1); S = hypotf(x, y - 1); A = (R + S) / 2; if (A < 1) A = 1; if (A < A_crossover) { if (y == 1 && x < FLT_EPSILON * FLT_EPSILON / 128) { *rx = sqrtf(x); } else if (x >= FLT_EPSILON * fabsf(y - 1)) { Am1 = f(x, 1 + y, R) + f(x, 1 - y, S); *rx = log1pf(Am1 + sqrtf(Am1 * (A + 1))); } else if (y < 1) { *rx = x / sqrtf((1 - y) * (1 + y)); } else { *rx = log1pf((y - 1) + sqrtf((y - 1) * (y + 1))); } } else { *rx = logf(A + sqrtf(A * A - 1)); } *new_y = y; if (y < FOUR_SQRT_MIN) { *B_is_usable = 0; *sqrt_A2my2 = A * (2 / FLT_EPSILON); *new_y = y * (2 / FLT_EPSILON); return; } *B = y / A; *B_is_usable = 1; if (*B > B_crossover) { *B_is_usable = 0; if (y == 1 && x < FLT_EPSILON / 128) { *sqrt_A2my2 = sqrtf(x) * sqrtf((A + y) / 2); } else if (x >= FLT_EPSILON * fabsf(y - 1)) { Amy = f(x, y + 1, R) + f(x, y - 1, S); *sqrt_A2my2 = sqrtf(Amy * (A + y)); } else if (y > 1) { *sqrt_A2my2 = x * (4 / FLT_EPSILON / FLT_EPSILON) * y / sqrtf((y + 1) * (y - 1)); *new_y = y * (4 / FLT_EPSILON / FLT_EPSILON); } else { *sqrt_A2my2 = sqrtf((1 - y) * (1 + y)); } } } float complex casinhf(float complex z) { float x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; int B_is_usable; float complex w; x = crealf(z); y = cimagf(z); ax = fabsf(x); ay = fabsf(y); if (isnan(x) || isnan(y)) { if (isinf(x)) return (CMPLXF(x, y + y)); if (isinf(y)) return (CMPLXF(y, x + x)); if (y == 0) return (CMPLXF(x + x, y)); - return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXF(nan_mix(x, y), nan_mix(x, y))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { if (signbit(x) == 0) w = clog_for_large_values(z) + m_ln2; else w = clog_for_large_values(-z) + m_ln2; return (CMPLXF(copysignf(crealf(w), x), copysignf(cimagf(w), y))); } if (x == 0 && y == 0) return (z); raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) return (z); do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); if (B_is_usable) ry = asinf(B); else ry = atan2f(new_y, sqrt_A2my2); return (CMPLXF(copysignf(rx, x), copysignf(ry, y))); } float complex casinf(float complex z) { float complex w = casinhf(CMPLXF(cimagf(z), crealf(z))); return (CMPLXF(cimagf(w), crealf(w))); } float complex cacosf(float complex z) { float x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; int sx, sy; int B_is_usable; float complex w; x = crealf(z); y = cimagf(z); sx = signbit(x); sy = signbit(y); ax = fabsf(x); ay = fabsf(y); if (isnan(x) || isnan(y)) { if (isinf(x)) return (CMPLXF(y + y, -INFINITY)); if (isinf(y)) return (CMPLXF(x + x, -y)); if (x == 0) return (CMPLXF(pio2_hi + pio2_lo, y + y)); - return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXF(nan_mix(x, y), nan_mix(x, y))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { w = clog_for_large_values(z); rx = fabsf(cimagf(w)); ry = crealf(w) + m_ln2; if (sy == 0) ry = -ry; return (CMPLXF(rx, ry)); } if (x == 1 && y == 0) return (CMPLXF(0, -y)); raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) return (CMPLXF(pio2_hi - (x - pio2_lo), -y)); do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); if (B_is_usable) { if (sx == 0) rx = acosf(B); else rx = acosf(-B); } else { if (sx == 0) rx = atan2f(sqrt_A2mx2, new_x); else rx = atan2f(sqrt_A2mx2, -new_x); } if (sy == 0) ry = -ry; return (CMPLXF(rx, ry)); } float complex cacoshf(float complex z) { float complex w; float rx, ry; w = cacosf(z); rx = crealf(w); ry = cimagf(w); if (isnan(rx) && isnan(ry)) return (CMPLXF(ry, rx)); if (isnan(rx)) return (CMPLXF(fabsf(ry), rx)); if (isnan(ry)) return (CMPLXF(ry, ry)); return (CMPLXF(fabsf(ry), copysignf(rx, cimagf(z)))); } static float complex clog_for_large_values(float complex z) { float x, y; float ax, ay, t; x = crealf(z); y = cimagf(z); ax = fabsf(x); ay = fabsf(y); if (ax < ay) { t = ax; ax = ay; ay = t; } if (ax > FLT_MAX / 2) return (CMPLXF(logf(hypotf(x / m_e, y / m_e)) + 1, atan2f(y, x))); if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) return (CMPLXF(logf(hypotf(x, y)), atan2f(y, x))); return (CMPLXF(logf(ax * ax + ay * ay) / 2, atan2f(y, x))); } static inline float sum_squares(float x, float y) { if (y < SQRT_MIN) return (x * x); return (x * x + y * y); } static inline float real_part_reciprocal(float x, float y) { float scale; uint32_t hx, hy; int32_t ix, iy; GET_FLOAT_WORD(hx, x); ix = hx & 0x7f800000; GET_FLOAT_WORD(hy, y); iy = hy & 0x7f800000; #define BIAS (FLT_MAX_EXP - 1) #define CUTOFF (FLT_MANT_DIG / 2 + 1) if (ix - iy >= CUTOFF << 23 || isinf(x)) return (1 / x); if (iy - ix >= CUTOFF << 23) return (x / y / y); if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23) return (x / (x * x + y * y)); SET_FLOAT_WORD(scale, 0x7f800000 - ix); x *= scale; y *= scale; return (x / (x * x + y * y) * scale); } float complex catanhf(float complex z) { float x, y, ax, ay, rx, ry; x = crealf(z); y = cimagf(z); ax = fabsf(x); ay = fabsf(y); if (y == 0 && ax <= 1) return (CMPLXF(atanhf(x), y)); if (x == 0) return (CMPLXF(x, atanf(y))); if (isnan(x) || isnan(y)) { if (isinf(x)) return (CMPLXF(copysignf(0, x), y + y)); if (isinf(y)) return (CMPLXF(copysignf(0, x), copysignf(pio2_hi + pio2_lo, y))); - return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXF(nan_mix(x, y), nan_mix(x, y))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) return (CMPLXF(real_part_reciprocal(x, y), copysignf(pio2_hi + pio2_lo, y))); if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { raise_inexact(); return (z); } if (ax == 1 && ay < FLT_EPSILON) rx = (m_ln2 - logf(ay)) / 2; else rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4; if (ax == 1) ry = atan2f(2, -ay) / 2; else if (ay < FLT_EPSILON) ry = atan2f(2 * ay, (1 - ax) * (1 + ax)) / 2; else ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; return (CMPLXF(copysignf(rx, x), copysignf(ry, y))); } float complex catanf(float complex z) { float complex w = catanhf(CMPLXF(cimagf(z), crealf(z))); return (CMPLXF(cimagf(w), crealf(w))); } Index: head/lib/msun/src/catrigl.c =================================================================== --- head/lib/msun/src/catrigl.c (revision 336361) +++ head/lib/msun/src/catrigl.c (revision 336362) @@ -1,417 +1,417 @@ /*- * Copyright (c) 2012 Stephen Montgomery-Smith * Copyright (c) 2017 Mahdi Mokhtari * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ /* * The algorithm is very close to that in "Implementing the complex arcsine * and arccosine functions using exception handling" by T. E. Hull, Thomas F. * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335, * http://dl.acm.org/citation.cfm?id=275324. * * See catrig.c for complete comments. * * XXX comments were removed automatically, and even short ones on the right * of statements were removed (all of them), contrary to normal style. Only * a few comments on the right of declarations remain. */ #include __FBSDID("$FreeBSD$"); #include #include #include "invtrig.h" #include "math.h" #include "math_private.h" #undef isinf #define isinf(x) (fabsl(x) == INFINITY) #undef isnan #define isnan(x) ((x) != (x)) #define raise_inexact() do { volatile float junk __unused = 1 + tiny; } while(0) #undef signbit #define signbit(x) (__builtin_signbitl(x)) #if LDBL_MAX_EXP != 0x4000 #error "Unsupported long double format" #endif static const long double A_crossover = 10, B_crossover = 0.6417, FOUR_SQRT_MIN = 0x1p-8189L, HALF_MAX = 0x1p16383L, QUARTER_SQRT_MAX = 0x1p8189L, RECIP_EPSILON = 1 / LDBL_EPSILON, SQRT_MIN = 0x1p-8191L; #if LDBL_MANT_DIG == 64 static const union IEEEl2bits um_e = LD80C(0xadf85458a2bb4a9b, 1, 2.71828182845904523536e+0L), um_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L); #define m_e um_e.e #define m_ln2 um_ln2.e static const long double /* The next 2 literals for non-i386. Misrounding them on i386 is harmless. */ SQRT_3_EPSILON = 5.70316273435758915310e-10, /* 0x9cc470a0490973e8.0p-94 */ SQRT_6_EPSILON = 8.06549008734932771664e-10; /* 0xddb3d742c265539e.0p-94 */ #elif LDBL_MANT_DIG == 113 static const long double m_e = 2.71828182845904523536028747135266250e0L, /* 0x15bf0a8b1457695355fb8ac404e7a.0p-111 */ m_ln2 = 6.93147180559945309417232121458176568e-1L, /* 0x162e42fefa39ef35793c7673007e6.0p-113 */ SQRT_3_EPSILON = 2.40370335797945490975336727199878124e-17, /* 0x1bb67ae8584caa73b25742d7078b8.0p-168 */ SQRT_6_EPSILON = 3.39934988877629587239082586223300391e-17; /* 0x13988e1409212e7d0321914321a55.0p-167 */ #else #error "Unsupported long double format" #endif static const volatile float tiny = 0x1p-100; static long double complex clog_for_large_values(long double complex z); static inline long double f(long double a, long double b, long double hypot_a_b) { if (b < 0) return ((hypot_a_b - b) / 2); if (b == 0) return (a / 2); return (a * a / (hypot_a_b + b) / 2); } static inline void do_hard_work(long double x, long double y, long double *rx, int *B_is_usable, long double *B, long double *sqrt_A2my2, long double *new_y) { long double R, S, A; long double Am1, Amy; R = hypotl(x, y + 1); S = hypotl(x, y - 1); A = (R + S) / 2; if (A < 1) A = 1; if (A < A_crossover) { if (y == 1 && x < LDBL_EPSILON * LDBL_EPSILON / 128) { *rx = sqrtl(x); } else if (x >= LDBL_EPSILON * fabsl(y - 1)) { Am1 = f(x, 1 + y, R) + f(x, 1 - y, S); *rx = log1pl(Am1 + sqrtl(Am1 * (A + 1))); } else if (y < 1) { *rx = x / sqrtl((1 - y) * (1 + y)); } else { *rx = log1pl((y - 1) + sqrtl((y - 1) * (y + 1))); } } else { *rx = logl(A + sqrtl(A * A - 1)); } *new_y = y; if (y < FOUR_SQRT_MIN) { *B_is_usable = 0; *sqrt_A2my2 = A * (2 / LDBL_EPSILON); *new_y = y * (2 / LDBL_EPSILON); return; } *B = y / A; *B_is_usable = 1; if (*B > B_crossover) { *B_is_usable = 0; if (y == 1 && x < LDBL_EPSILON / 128) { *sqrt_A2my2 = sqrtl(x) * sqrtl((A + y) / 2); } else if (x >= LDBL_EPSILON * fabsl(y - 1)) { Amy = f(x, y + 1, R) + f(x, y - 1, S); *sqrt_A2my2 = sqrtl(Amy * (A + y)); } else if (y > 1) { *sqrt_A2my2 = x * (4 / LDBL_EPSILON / LDBL_EPSILON) * y / sqrtl((y + 1) * (y - 1)); *new_y = y * (4 / LDBL_EPSILON / LDBL_EPSILON); } else { *sqrt_A2my2 = sqrtl((1 - y) * (1 + y)); } } } long double complex casinhl(long double complex z) { long double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; int B_is_usable; long double complex w; x = creall(z); y = cimagl(z); ax = fabsl(x); ay = fabsl(y); if (isnan(x) || isnan(y)) { if (isinf(x)) return (CMPLXL(x, y + y)); if (isinf(y)) return (CMPLXL(y, x + x)); if (y == 0) return (CMPLXL(x + x, y)); - return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXL(nan_mix(x, y), nan_mix(x, y))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { if (signbit(x) == 0) w = clog_for_large_values(z) + m_ln2; else w = clog_for_large_values(-z) + m_ln2; return (CMPLXL(copysignl(creall(w), x), copysignl(cimagl(w), y))); } if (x == 0 && y == 0) return (z); raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) return (z); do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); if (B_is_usable) ry = asinl(B); else ry = atan2l(new_y, sqrt_A2my2); return (CMPLXL(copysignl(rx, x), copysignl(ry, y))); } long double complex casinl(long double complex z) { long double complex w; w = casinhl(CMPLXL(cimagl(z), creall(z))); return (CMPLXL(cimagl(w), creall(w))); } long double complex cacosl(long double complex z) { long double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; int sx, sy; int B_is_usable; long double complex w; x = creall(z); y = cimagl(z); sx = signbit(x); sy = signbit(y); ax = fabsl(x); ay = fabsl(y); if (isnan(x) || isnan(y)) { if (isinf(x)) return (CMPLXL(y + y, -INFINITY)); if (isinf(y)) return (CMPLXL(x + x, -y)); if (x == 0) return (CMPLXL(pio2_hi + pio2_lo, y + y)); - return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXL(nan_mix(x, y), nan_mix(x, y))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { w = clog_for_large_values(z); rx = fabsl(cimagl(w)); ry = creall(w) + m_ln2; if (sy == 0) ry = -ry; return (CMPLXL(rx, ry)); } if (x == 1 && y == 0) return (CMPLXL(0, -y)); raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) return (CMPLXL(pio2_hi - (x - pio2_lo), -y)); do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); if (B_is_usable) { if (sx == 0) rx = acosl(B); else rx = acosl(-B); } else { if (sx == 0) rx = atan2l(sqrt_A2mx2, new_x); else rx = atan2l(sqrt_A2mx2, -new_x); } if (sy == 0) ry = -ry; return (CMPLXL(rx, ry)); } long double complex cacoshl(long double complex z) { long double complex w; long double rx, ry; w = cacosl(z); rx = creall(w); ry = cimagl(w); if (isnan(rx) && isnan(ry)) return (CMPLXL(ry, rx)); if (isnan(rx)) return (CMPLXL(fabsl(ry), rx)); if (isnan(ry)) return (CMPLXL(ry, ry)); return (CMPLXL(fabsl(ry), copysignl(rx, cimagl(z)))); } static long double complex clog_for_large_values(long double complex z) { long double x, y; long double ax, ay, t; x = creall(z); y = cimagl(z); ax = fabsl(x); ay = fabsl(y); if (ax < ay) { t = ax; ax = ay; ay = t; } if (ax > HALF_MAX) return (CMPLXL(logl(hypotl(x / m_e, y / m_e)) + 1, atan2l(y, x))); if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) return (CMPLXL(logl(hypotl(x, y)), atan2l(y, x))); return (CMPLXL(logl(ax * ax + ay * ay) / 2, atan2l(y, x))); } static inline long double sum_squares(long double x, long double y) { if (y < SQRT_MIN) return (x * x); return (x * x + y * y); } static inline long double real_part_reciprocal(long double x, long double y) { long double scale; uint16_t hx, hy; int16_t ix, iy; GET_LDBL_EXPSIGN(hx, x); ix = hx & 0x7fff; GET_LDBL_EXPSIGN(hy, y); iy = hy & 0x7fff; #define BIAS (LDBL_MAX_EXP - 1) #define CUTOFF (LDBL_MANT_DIG / 2 + 1) if (ix - iy >= CUTOFF || isinf(x)) return (1 / x); if (iy - ix >= CUTOFF) return (x / y / y); if (ix <= BIAS + LDBL_MAX_EXP / 2 - CUTOFF) return (x / (x * x + y * y)); scale = 1; SET_LDBL_EXPSIGN(scale, 0x7fff - ix); x *= scale; y *= scale; return (x / (x * x + y * y) * scale); } long double complex catanhl(long double complex z) { long double x, y, ax, ay, rx, ry; x = creall(z); y = cimagl(z); ax = fabsl(x); ay = fabsl(y); if (y == 0 && ax <= 1) return (CMPLXL(atanhl(x), y)); if (x == 0) return (CMPLXL(x, atanl(y))); if (isnan(x) || isnan(y)) { if (isinf(x)) return (CMPLXL(copysignl(0, x), y + y)); if (isinf(y)) return (CMPLXL(copysignl(0, x), copysignl(pio2_hi + pio2_lo, y))); - return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXL(nan_mix(x, y), nan_mix(x, y))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) return (CMPLXL(real_part_reciprocal(x, y), copysignl(pio2_hi + pio2_lo, y))); if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { raise_inexact(); return (z); } if (ax == 1 && ay < LDBL_EPSILON) rx = (m_ln2 - logl(ay)) / 2; else rx = log1pl(4 * ax / sum_squares(ax - 1, ay)) / 4; if (ax == 1) ry = atan2l(2, -ay) / 2; else if (ay < LDBL_EPSILON) ry = atan2l(2 * ay, (1 - ax) * (1 + ax)) / 2; else ry = atan2l(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; return (CMPLXL(copysignl(rx, x), copysignl(ry, y))); } long double complex catanl(long double complex z) { long double complex w; w = catanhl(CMPLXL(cimagl(z), creall(z))); return (CMPLXL(cimagl(w), creall(w))); } Index: head/lib/msun/src/e_atan2.c =================================================================== --- head/lib/msun/src/e_atan2.c (revision 336361) +++ head/lib/msun/src/e_atan2.c (revision 336362) @@ -1,129 +1,129 @@ /* @(#)e_atan2.c 1.3 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== * */ #include __FBSDID("$FreeBSD$"); /* __ieee754_atan2(y,x) * Method : * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). * 2. Reduce x to positive by (if x and y are unexceptional): * ARG (x+iy) = arctan(y/x) ... if x > 0, * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, * * Special cases: * * ATAN2((anything), NaN ) is NaN; * ATAN2(NAN , (anything) ) is NaN; * ATAN2(+-0, +(anything but NaN)) is +-0 ; * ATAN2(+-0, -(anything but NaN)) is +-pi ; * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2; * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ; * ATAN2(+-(anything but INF and NaN), -INF) is +-pi; * ATAN2(+-INF,+INF ) is +-pi/4 ; * ATAN2(+-INF,-INF ) is +-3pi/4; * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2; * * Constants: * The hexadecimal values are the intended ones for the following * constants. The decimal values may be used, provided that the * compiler will convert from decimal to binary accurately enough * to produce the hexadecimal values shown. */ #include #include "math.h" #include "math_private.h" static volatile double tiny = 1.0e-300; static const double zero = 0.0, pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */ pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */ pi = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */ static volatile double pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ double __ieee754_atan2(double y, double x) { double z; int32_t k,m,hx,hy,ix,iy; u_int32_t lx,ly; EXTRACT_WORDS(hx,lx,x); ix = hx&0x7fffffff; EXTRACT_WORDS(hy,ly,y); iy = hy&0x7fffffff; if(((ix|((lx|-lx)>>31))>0x7ff00000)|| ((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */ - return x+y; + return nan_mix(x, y); if(hx==0x3ff00000&&lx==0) return atan(y); /* x=1.0 */ m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */ /* when y = 0 */ if((iy|ly)==0) { switch(m) { case 0: case 1: return y; /* atan(+-0,+anything)=+-0 */ case 2: return pi+tiny;/* atan(+0,-anything) = pi */ case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ } } /* when x = 0 */ if((ix|lx)==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; /* when x is INF */ if(ix==0x7ff00000) { if(iy==0x7ff00000) { switch(m) { case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */ case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */ case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/ case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/ } } else { switch(m) { case 0: return zero ; /* atan(+...,+INF) */ case 1: return -zero ; /* atan(-...,+INF) */ case 2: return pi+tiny ; /* atan(+...,-INF) */ case 3: return -pi-tiny ; /* atan(-...,-INF) */ } } } /* when y is INF */ if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; /* compute y/x */ k = (iy-ix)>>20; if(k > 60) { /* |y/x| > 2**60 */ z=pi_o_2+0.5*pi_lo; m&=1; } else if(hx<0&&k<-60) z=0.0; /* 0 > |y|/x > -2**-60 */ else z=atan(fabs(y/x)); /* safe to do y/x */ switch (m) { case 0: return z ; /* atan(+,+) */ case 1: return -z ; /* atan(-,+) */ case 2: return pi-(z-pi_lo);/* atan(+,-) */ default: /* case 3 */ return (z-pi_lo)-pi;/* atan(-,-) */ } } #if LDBL_MANT_DIG == 53 __weak_reference(atan2, atan2l); #endif Index: head/lib/msun/src/e_atan2f.c =================================================================== --- head/lib/msun/src/e_atan2f.c (revision 336361) +++ head/lib/msun/src/e_atan2f.c (revision 336362) @@ -1,96 +1,96 @@ /* e_atan2f.c -- float version of e_atan2.c. * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" static volatile float tiny = 1.0e-30; static const float zero = 0.0, pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */ pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */ pi = 3.1415927410e+00; /* 0x40490fdb */ static volatile float pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */ float __ieee754_atan2f(float y, float x) { float z; int32_t k,m,hx,hy,ix,iy; GET_FLOAT_WORD(hx,x); ix = hx&0x7fffffff; GET_FLOAT_WORD(hy,y); iy = hy&0x7fffffff; if((ix>0x7f800000)|| (iy>0x7f800000)) /* x or y is NaN */ - return x+y; + return nan_mix(x, y); if(hx==0x3f800000) return atanf(y); /* x=1.0 */ m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */ /* when y = 0 */ if(iy==0) { switch(m) { case 0: case 1: return y; /* atan(+-0,+anything)=+-0 */ case 2: return pi+tiny;/* atan(+0,-anything) = pi */ case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ } } /* when x = 0 */ if(ix==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; /* when x is INF */ if(ix==0x7f800000) { if(iy==0x7f800000) { switch(m) { case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */ case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */ case 2: return (float)3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/ case 3: return (float)-3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/ } } else { switch(m) { case 0: return zero ; /* atan(+...,+INF) */ case 1: return -zero ; /* atan(-...,+INF) */ case 2: return pi+tiny ; /* atan(+...,-INF) */ case 3: return -pi-tiny ; /* atan(-...,-INF) */ } } } /* when y is INF */ if(iy==0x7f800000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; /* compute y/x */ k = (iy-ix)>>23; if(k > 26) { /* |y/x| > 2**26 */ z=pi_o_2+(float)0.5*pi_lo; m&=1; } else if(k<-26&&hx<0) z=0.0; /* 0 > |y|/x > -2**-26 */ else z=atanf(fabsf(y/x)); /* safe to do y/x */ switch (m) { case 0: return z ; /* atan(+,+) */ case 1: return -z ; /* atan(-,+) */ case 2: return pi-(z-pi_lo);/* atan(+,-) */ default: /* case 3 */ return (z-pi_lo)-pi;/* atan(-,-) */ } } Index: head/lib/msun/src/e_atan2l.c =================================================================== --- head/lib/msun/src/e_atan2l.c (revision 336361) +++ head/lib/msun/src/e_atan2l.c (revision 336362) @@ -1,120 +1,120 @@ /* @(#)e_atan2.c 1.3 95/01/18 */ /* FreeBSD: head/lib/msun/src/e_atan2.c 176451 2008-02-22 02:30:36Z das */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== * */ #include __FBSDID("$FreeBSD$"); /* * See comments in e_atan2.c. * Converted to long double by David Schultz . */ #include #include "invtrig.h" #include "math.h" #include "math_private.h" static volatile long double tiny = 1.0e-300; static const long double zero = 0.0; #ifdef __i386__ /* XXX Work around the fact that gcc truncates long double constants on i386 */ static volatile double pi1 = 3.14159265358979311600e+00, /* 0x1.921fb54442d18p+1 */ pi2 = 1.22514845490862001043e-16; /* 0x1.1a80000000000p-53 */ #define pi ((long double)pi1 + pi2) #else static const long double pi = 3.14159265358979323846264338327950280e+00L; #endif long double atan2l(long double y, long double x) { union IEEEl2bits ux, uy; long double z; int32_t k,m; int16_t exptx, expsignx, expty, expsigny; uy.e = y; expsigny = uy.xbits.expsign; expty = expsigny & 0x7fff; ux.e = x; expsignx = ux.xbits.expsign; exptx = expsignx & 0x7fff; if ((exptx==BIAS+LDBL_MAX_EXP && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)!=0) || /* x is NaN */ (expty==BIAS+LDBL_MAX_EXP && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* y is NaN */ - return x+y; + return nan_mix(x, y); if (expsignx==BIAS && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0) return atanl(y); /* x=1.0 */ m = ((expsigny>>15)&1)|((expsignx>>14)&2); /* 2*sign(x)+sign(y) */ /* when y = 0 */ if(expty==0 && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)==0) { switch(m) { case 0: case 1: return y; /* atan(+-0,+anything)=+-0 */ case 2: return pi+tiny;/* atan(+0,-anything) = pi */ case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ } } /* when x = 0 */ if(exptx==0 && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl)==0) return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; /* when x is INF */ if(exptx==BIAS+LDBL_MAX_EXP) { if(expty==BIAS+LDBL_MAX_EXP) { switch(m) { case 0: return pio2_hi*0.5+tiny;/* atan(+INF,+INF) */ case 1: return -pio2_hi*0.5-tiny;/* atan(-INF,+INF) */ case 2: return 1.5*pio2_hi+tiny;/*atan(+INF,-INF)*/ case 3: return -1.5*pio2_hi-tiny;/*atan(-INF,-INF)*/ } } else { switch(m) { case 0: return zero ; /* atan(+...,+INF) */ case 1: return -zero ; /* atan(-...,+INF) */ case 2: return pi+tiny ; /* atan(+...,-INF) */ case 3: return -pi-tiny ; /* atan(-...,-INF) */ } } } /* when y is INF */ if(expty==BIAS+LDBL_MAX_EXP) return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny; /* compute y/x */ k = expty-exptx; if(k > LDBL_MANT_DIG+2) { /* |y/x| huge */ z=pio2_hi+pio2_lo; m&=1; } else if(expsignx<0&&k<-LDBL_MANT_DIG-2) z=0.0; /* |y/x| tiny, x<0 */ else z=atanl(fabsl(y/x)); /* safe to do y/x */ switch (m) { case 0: return z ; /* atan(+,+) */ case 1: return -z ; /* atan(-,+) */ case 2: return pi-(z-pi_lo);/* atan(+,-) */ default: /* case 3 */ return (z-pi_lo)-pi;/* atan(-,-) */ } } Index: head/lib/msun/src/e_fmod.c =================================================================== --- head/lib/msun/src/e_fmod.c (revision 336361) +++ head/lib/msun/src/e_fmod.c (revision 336362) @@ -1,138 +1,138 @@ /* @(#)e_fmod.c 1.3 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); /* * __ieee754_fmod(x,y) * Return x mod y in exact arithmetic * Method: shift and subtract */ #include #include "math.h" #include "math_private.h" static const double one = 1.0, Zero[] = {0.0, -0.0,}; double __ieee754_fmod(double x, double y) { int32_t n,hx,hy,hz,ix,iy,sx,i; u_int32_t lx,ly,lz; EXTRACT_WORDS(hx,lx,x); EXTRACT_WORDS(hy,ly,y); sx = hx&0x80000000; /* sign of x */ hx ^=sx; /* |x| */ hy &= 0x7fffffff; /* |y| */ /* purge off exception values */ if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */ - return (x*y)/(x*y); + return nan_mix(x, y)/nan_mix(x, y); if(hx<=hy) { if((hx>31]; /* |x|=|y| return x*0*/ } /* determine ix = ilogb(x) */ if(hx<0x00100000) { /* subnormal x */ if(hx==0) { for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; } else { for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; } } else ix = (hx>>20)-1023; /* determine iy = ilogb(y) */ if(hy<0x00100000) { /* subnormal y */ if(hy==0) { for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; } else { for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; } } else iy = (hy>>20)-1023; /* set up {hx,lx}, {hy,ly} and align y to x */ if(ix >= -1022) hx = 0x00100000|(0x000fffff&hx); else { /* subnormal x, shift x to normal */ n = -1022-ix; if(n<=31) { hx = (hx<>(32-n)); lx <<= n; } else { hx = lx<<(n-32); lx = 0; } } if(iy >= -1022) hy = 0x00100000|(0x000fffff&hy); else { /* subnormal y, shift y to normal */ n = -1022-iy; if(n<=31) { hy = (hy<>(32-n)); ly <<= n; } else { hy = ly<<(n-32); ly = 0; } } /* fix point fmod */ n = ix - iy; while(n--) { hz=hx-hy;lz=lx-ly; if(lx>31); lx = lx+lx;} else { if((hz|lz)==0) /* return sign(x)*0 */ return Zero[(u_int32_t)sx>>31]; hx = hz+hz+(lz>>31); lx = lz+lz; } } hz=hx-hy;lz=lx-ly; if(lx=0) {hx=hz;lx=lz;} /* convert back to floating value and restore the sign */ if((hx|lx)==0) /* return sign(x)*0 */ return Zero[(u_int32_t)sx>>31]; while(hx<0x00100000) { /* normalize x */ hx = hx+hx+(lx>>31); lx = lx+lx; iy -= 1; } if(iy>= -1022) { /* normalize output */ hx = ((hx-0x00100000)|((iy+1023)<<20)); INSERT_WORDS(x,hx|sx,lx); } else { /* subnormal output */ n = -1022 - iy; if(n<=20) { lx = (lx>>n)|((u_int32_t)hx<<(32-n)); hx >>= n; } else if (n<=31) { lx = (hx<<(32-n))|(lx>>n); hx = sx; } else { lx = hx>>(n-32); hx = sx; } INSERT_WORDS(x,hx|sx,lx); x *= one; /* create necessary signal */ } return x; /* exact output */ } #if (LDBL_MANT_DIG == 53) __weak_reference(fmod, fmodl); #endif Index: head/lib/msun/src/e_fmodf.c =================================================================== --- head/lib/msun/src/e_fmodf.c (revision 336361) +++ head/lib/msun/src/e_fmodf.c (revision 336362) @@ -1,104 +1,104 @@ /* e_fmodf.c -- float version of e_fmod.c. * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); /* * __ieee754_fmodf(x,y) * Return x mod y in exact arithmetic * Method: shift and subtract */ #include "math.h" #include "math_private.h" static const float one = 1.0, Zero[] = {0.0, -0.0,}; float __ieee754_fmodf(float x, float y) { int32_t n,hx,hy,hz,ix,iy,sx,i; GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hy,y); sx = hx&0x80000000; /* sign of x */ hx ^=sx; /* |x| */ hy &= 0x7fffffff; /* |y| */ /* purge off exception values */ if(hy==0||(hx>=0x7f800000)|| /* y=0,or x not finite */ (hy>0x7f800000)) /* or y is NaN */ - return (x*y)/(x*y); + return nan_mix(x, y)/nan_mix(x, y); if(hx>31]; /* |x|=|y| return x*0*/ /* determine ix = ilogb(x) */ if(hx<0x00800000) { /* subnormal x */ for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1; } else ix = (hx>>23)-127; /* determine iy = ilogb(y) */ if(hy<0x00800000) { /* subnormal y */ for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1; } else iy = (hy>>23)-127; /* set up {hx,lx}, {hy,ly} and align y to x */ if(ix >= -126) hx = 0x00800000|(0x007fffff&hx); else { /* subnormal x, shift x to normal */ n = -126-ix; hx = hx<= -126) hy = 0x00800000|(0x007fffff&hy); else { /* subnormal y, shift y to normal */ n = -126-iy; hy = hy<>31]; hx = hz+hz; } } hz=hx-hy; if(hz>=0) {hx=hz;} /* convert back to floating value and restore the sign */ if(hx==0) /* return sign(x)*0 */ return Zero[(u_int32_t)sx>>31]; while(hx<0x00800000) { /* normalize x */ hx = hx+hx; iy -= 1; } if(iy>= -126) { /* normalize output */ hx = ((hx-0x00800000)|((iy+127)<<23)); SET_FLOAT_WORD(x,hx|sx); } else { /* subnormal output */ n = -126 - iy; hx >>= n; SET_FLOAT_WORD(x,hx|sx); x *= one; /* create necessary signal */ } return x; /* exact output */ } Index: head/lib/msun/src/e_fmodl.c =================================================================== --- head/lib/msun/src/e_fmodl.c (revision 336361) +++ head/lib/msun/src/e_fmodl.c (revision 336362) @@ -1,149 +1,149 @@ /* @(#)e_fmod.c 1.3 95/01/18 */ /*- * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); #include #include #include "fpmath.h" #include "math.h" #include "math_private.h" #define BIAS (LDBL_MAX_EXP - 1) #if LDBL_MANL_SIZE > 32 typedef uint64_t manl_t; #else typedef uint32_t manl_t; #endif #if LDBL_MANH_SIZE > 32 typedef uint64_t manh_t; #else typedef uint32_t manh_t; #endif /* * These macros add and remove an explicit integer bit in front of the * fractional mantissa, if the architecture doesn't have such a bit by * default already. */ #ifdef LDBL_IMPLICIT_NBIT #define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE)) #define HFRAC_BITS LDBL_MANH_SIZE #else #define SET_NBIT(hx) (hx) #define HFRAC_BITS (LDBL_MANH_SIZE - 1) #endif #define MANL_SHIFT (LDBL_MANL_SIZE - 1) static const long double one = 1.0, Zero[] = {0.0, -0.0,}; /* * fmodl(x,y) * Return x mod y in exact arithmetic * Method: shift and subtract * * Assumptions: * - The low part of the mantissa fits in a manl_t exactly. * - The high part of the mantissa fits in an int64_t with enough room * for an explicit integer bit in front of the fractional bits. */ long double fmodl(long double x, long double y) { union IEEEl2bits ux, uy; int64_t hx,hz; /* We need a carry bit even if LDBL_MANH_SIZE is 32. */ manh_t hy; manl_t lx,ly,lz; int ix,iy,n,sx; ux.e = x; uy.e = y; sx = ux.bits.sign; /* purge off exception values */ if((uy.bits.exp|uy.bits.manh|uy.bits.manl)==0 || /* y=0 */ (ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */ (uy.bits.exp == BIAS + LDBL_MAX_EXP && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */ - return (x*y)/(x*y); + return nan_mix(x, y)/nan_mix(x, y); if(ux.bits.exp<=uy.bits.exp) { if((ux.bits.exp>MANL_SHIFT); lx = lx+lx;} else { if ((hz|lz)==0) /* return sign(x)*0 */ return Zero[sx]; hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz; } } hz=hx-hy;lz=lx-ly; if(lx=0) {hx=hz;lx=lz;} /* convert back to floating value and restore the sign */ if((hx|lx)==0) /* return sign(x)*0 */ return Zero[sx]; while(hx<(1ULL<>MANL_SHIFT); lx = lx+lx; iy -= 1; } ux.bits.manh = hx; /* The mantissa is truncated here if needed. */ ux.bits.manl = lx; if (iy < LDBL_MIN_EXP) { ux.bits.exp = iy + (BIAS + 512); ux.e *= 0x1p-512; } else { ux.bits.exp = iy + BIAS; } x = ux.e * one; /* create necessary signal */ return x; /* exact output */ } Index: head/lib/msun/src/e_hypot.c =================================================================== --- head/lib/msun/src/e_hypot.c (revision 336361) +++ head/lib/msun/src/e_hypot.c (revision 336362) @@ -1,131 +1,131 @@ /* @(#)e_hypot.c 1.3 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); /* __ieee754_hypot(x,y) * * Method : * If (assume round-to-nearest) z=x*x+y*y * has error less than sqrt(2)/2 ulp, than * sqrt(z) has error less than 1 ulp (exercise). * * So, compute sqrt(x*x+y*y) with some care as * follows to get the error below 1 ulp: * * Assume x>y>0; * (if possible, set rounding to round-to-nearest) * 1. if x > 2y use * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y * where x1 = x with lower 32 bits cleared, x2 = x-x1; else * 2. if x <= 2y use * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y)) * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, * y1= y with lower 32 bits chopped, y2 = y-y1. * * NOTE: scaling may be necessary if some argument is too * large or too tiny * * Special cases: * hypot(x,y) is INF if x or y is +INF or -INF; else * hypot(x,y) is NAN if x or y is NAN. * * Accuracy: * hypot(x,y) returns sqrt(x^2+y^2) with error less * than 1 ulps (units in the last place) */ #include #include "math.h" #include "math_private.h" double __ieee754_hypot(double x, double y) { double a,b,t1,t2,y1,y2,w; int32_t j,k,ha,hb; GET_HIGH_WORD(ha,x); ha &= 0x7fffffff; GET_HIGH_WORD(hb,y); hb &= 0x7fffffff; if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} a = fabs(a); b = fabs(b); if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */ k=0; if(ha > 0x5f300000) { /* a>2**500 */ if(ha >= 0x7ff00000) { /* Inf or NaN */ u_int32_t low; /* Use original arg order iff result is NaN; quieten sNaNs. */ - w = fabs(x+0.0)-fabs(y+0.0); + w = fabsl(x+0.0L)-fabs(y+0); GET_LOW_WORD(low,a); if(((ha&0xfffff)|low)==0) w = a; GET_LOW_WORD(low,b); if(((hb^0x7ff00000)|low)==0) w = b; return w; } /* scale a and b by 2**-600 */ ha -= 0x25800000; hb -= 0x25800000; k += 600; SET_HIGH_WORD(a,ha); SET_HIGH_WORD(b,hb); } if(hb < 0x20b00000) { /* b < 2**-500 */ if(hb <= 0x000fffff) { /* subnormal b or 0 */ u_int32_t low; GET_LOW_WORD(low,b); if((hb|low)==0) return a; t1=0; SET_HIGH_WORD(t1,0x7fd00000); /* t1=2^1022 */ b *= t1; a *= t1; k -= 1022; } else { /* scale a and b by 2^600 */ ha += 0x25800000; /* a *= 2^600 */ hb += 0x25800000; /* b *= 2^600 */ k -= 600; SET_HIGH_WORD(a,ha); SET_HIGH_WORD(b,hb); } } /* medium size a and b */ w = a-b; if (w>b) { t1 = 0; SET_HIGH_WORD(t1,ha); t2 = a-t1; w = sqrt(t1*t1-(b*(-b)-t2*(a+t1))); } else { a = a+a; y1 = 0; SET_HIGH_WORD(y1,hb); y2 = b - y1; t1 = 0; SET_HIGH_WORD(t1,ha+0x00100000); t2 = a - t1; w = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b))); } if(k!=0) { u_int32_t high; t1 = 1.0; GET_HIGH_WORD(high,t1); SET_HIGH_WORD(t1,high+(k<<20)); return t1*w; } else return w; } #if LDBL_MANT_DIG == 53 __weak_reference(hypot, hypotl); #endif Index: head/lib/msun/src/e_hypotf.c =================================================================== --- head/lib/msun/src/e_hypotf.c (revision 336361) +++ head/lib/msun/src/e_hypotf.c (revision 336362) @@ -1,83 +1,83 @@ /* e_hypotf.c -- float version of e_hypot.c. * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" float __ieee754_hypotf(float x, float y) { float a,b,t1,t2,y1,y2,w; int32_t j,k,ha,hb; GET_FLOAT_WORD(ha,x); ha &= 0x7fffffff; GET_FLOAT_WORD(hb,y); hb &= 0x7fffffff; if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} a = fabsf(a); b = fabsf(b); if((ha-hb)>0xf000000) {return a+b;} /* x/y > 2**30 */ k=0; if(ha > 0x58800000) { /* a>2**50 */ if(ha >= 0x7f800000) { /* Inf or NaN */ /* Use original arg order iff result is NaN; quieten sNaNs. */ - w = fabsf(x+0.0F)-fabsf(y+0.0F); + w = fabsl(x+0.0L)-fabsf(y+0); if(ha == 0x7f800000) w = a; if(hb == 0x7f800000) w = b; return w; } /* scale a and b by 2**-68 */ ha -= 0x22000000; hb -= 0x22000000; k += 68; SET_FLOAT_WORD(a,ha); SET_FLOAT_WORD(b,hb); } if(hb < 0x26800000) { /* b < 2**-50 */ if(hb <= 0x007fffff) { /* subnormal b or 0 */ if(hb==0) return a; SET_FLOAT_WORD(t1,0x7e800000); /* t1=2^126 */ b *= t1; a *= t1; k -= 126; } else { /* scale a and b by 2^68 */ ha += 0x22000000; /* a *= 2^68 */ hb += 0x22000000; /* b *= 2^68 */ k -= 68; SET_FLOAT_WORD(a,ha); SET_FLOAT_WORD(b,hb); } } /* medium size a and b */ w = a-b; if (w>b) { SET_FLOAT_WORD(t1,ha&0xfffff000); t2 = a-t1; w = __ieee754_sqrtf(t1*t1-(b*(-b)-t2*(a+t1))); } else { a = a+a; SET_FLOAT_WORD(y1,hb&0xfffff000); y2 = b - y1; SET_FLOAT_WORD(t1,(ha+0x00800000)&0xfffff000); t2 = a - t1; w = __ieee754_sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b))); } if(k!=0) { SET_FLOAT_WORD(t1,0x3f800000+(k<<23)); return t1*w; } else return w; } Index: head/lib/msun/src/e_hypotl.c =================================================================== --- head/lib/msun/src/e_hypotl.c (revision 336361) +++ head/lib/msun/src/e_hypotl.c (revision 336362) @@ -1,124 +1,124 @@ /* From: @(#)e_hypot.c 1.3 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); /* long double version of hypot(). See e_hypot.c for most comments. */ #include #include "fpmath.h" #include "math.h" #include "math_private.h" #define GET_LDBL_MAN(h, l, v) do { \ union IEEEl2bits uv; \ \ uv.e = v; \ h = uv.bits.manh; \ l = uv.bits.manl; \ } while (0) #undef GET_HIGH_WORD #define GET_HIGH_WORD(i, v) GET_LDBL_EXPSIGN(i, v) #undef SET_HIGH_WORD #define SET_HIGH_WORD(v, i) SET_LDBL_EXPSIGN(v, i) #define DESW(exp) (exp) /* delta expsign word */ #define ESW(exp) (MAX_EXP - 1 + (exp)) /* expsign word */ #define MANT_DIG LDBL_MANT_DIG #define MAX_EXP LDBL_MAX_EXP #if LDBL_MANL_SIZE > 32 typedef uint64_t man_t; #else typedef uint32_t man_t; #endif long double hypotl(long double x, long double y) { long double a=x,b=y,t1,t2,y1,y2,w; int32_t j,k,ha,hb; GET_HIGH_WORD(ha,x); ha &= 0x7fff; GET_HIGH_WORD(hb,y); hb &= 0x7fff; if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} a = fabsl(a); b = fabsl(b); if((ha-hb)>DESW(MANT_DIG+7)) {return a+b;} /* x/y > 2**(MANT_DIG+7) */ k=0; if(ha > ESW(MAX_EXP/2-12)) { /* a>2**(MAX_EXP/2-12) */ if(ha >= ESW(MAX_EXP)) { /* Inf or NaN */ man_t manh, manl; /* Use original arg order iff result is NaN; quieten sNaNs. */ - w = fabsl(x+0.0)-fabsl(y+0.0); + w = fabsl(x+0.0L)-fabsl(y+0); GET_LDBL_MAN(manh,manl,a); if (manh == LDBL_NBIT && manl == 0) w = a; GET_LDBL_MAN(manh,manl,b); if (hb >= ESW(MAX_EXP) && manh == LDBL_NBIT && manl == 0) w = b; return w; } /* scale a and b by 2**-(MAX_EXP/2+88) */ ha -= DESW(MAX_EXP/2+88); hb -= DESW(MAX_EXP/2+88); k += MAX_EXP/2+88; SET_HIGH_WORD(a,ha); SET_HIGH_WORD(b,hb); } if(hb < ESW(-(MAX_EXP/2-12))) { /* b < 2**-(MAX_EXP/2-12) */ if(hb <= 0) { /* subnormal b or 0 */ man_t manh, manl; GET_LDBL_MAN(manh,manl,b); if((manh|manl)==0) return a; t1=0; SET_HIGH_WORD(t1,ESW(MAX_EXP-2)); /* t1=2^(MAX_EXP-2) */ b *= t1; a *= t1; k -= MAX_EXP-2; } else { /* scale a and b by 2^(MAX_EXP/2+88) */ ha += DESW(MAX_EXP/2+88); hb += DESW(MAX_EXP/2+88); k -= MAX_EXP/2+88; SET_HIGH_WORD(a,ha); SET_HIGH_WORD(b,hb); } } /* medium size a and b */ w = a-b; if (w>b) { t1 = a; union IEEEl2bits uv; uv.e = t1; uv.bits.manl = 0; t1 = uv.e; t2 = a-t1; w = sqrtl(t1*t1-(b*(-b)-t2*(a+t1))); } else { a = a+a; y1 = b; union IEEEl2bits uv; uv.e = y1; uv.bits.manl = 0; y1 = uv.e; y2 = b - y1; t1 = a; uv.e = t1; uv.bits.manl = 0; t1 = uv.e; t2 = a - t1; w = sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b))); } if(k!=0) { u_int32_t high; t1 = 1.0; GET_HIGH_WORD(high,t1); SET_HIGH_WORD(t1,high+DESW(k)); return t1*w; } else return w; } Index: head/lib/msun/src/e_pow.c =================================================================== --- head/lib/msun/src/e_pow.c (revision 336361) +++ head/lib/msun/src/e_pow.c (revision 336362) @@ -1,314 +1,314 @@ /* @(#)e_pow.c 1.5 04/04/22 SMI */ /* * ==================================================== * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. * * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); /* __ieee754_pow(x,y) return x**y * * n * Method: Let x = 2 * (1+f) * 1. Compute and return log2(x) in two pieces: * log2(x) = w1 + w2, * where w1 has 53-24 = 29 bit trailing zeros. * 2. Perform y*log2(x) = n+y' by simulating multi-precision * arithmetic, where |y'|<=0.5. * 3. Return x**y = 2**n*exp(y'*log2) * * Special cases: * 1. (anything) ** 0 is 1 * 2. (anything) ** 1 is itself * 3. (anything) ** NAN is NAN except 1 ** NAN = 1 * 4. NAN ** (anything except 0) is NAN * 5. +-(|x| > 1) ** +INF is +INF * 6. +-(|x| > 1) ** -INF is +0 * 7. +-(|x| < 1) ** +INF is +0 * 8. +-(|x| < 1) ** -INF is +INF * 9. +-1 ** +-INF is 1 * 10. +0 ** (+anything except 0, NAN) is +0 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 * 12. +0 ** (-anything except 0, NAN) is +INF * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF * 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) * 15. +INF ** (+anything except 0,NAN) is +INF * 16. +INF ** (-anything except 0,NAN) is +0 * 17. -INF ** (anything) = -0 ** (-anything) * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) * 19. (-anything except 0 and inf) ** (non-integer) is NAN * * Accuracy: * pow(x,y) returns x**y nearly rounded. In particular * pow(integer,integer) * always returns the correct integer provided it is * representable. * * Constants : * The hexadecimal values are the intended ones for the following * constants. The decimal values may be used, provided that the * compiler will convert from decimal to binary accurately enough * to produce the hexadecimal values shown. */ #include #include "math.h" #include "math_private.h" static const double bp[] = {1.0, 1.5,}, dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ zero = 0.0, half = 0.5, qrtr = 0.25, thrd = 3.3333333333333331e-01, /* 0x3fd55555, 0x55555555 */ one = 1.0, two = 2.0, two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ huge = 1.0e300, tiny = 1.0e-300, /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */ L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */ L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */ L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */ P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */ lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */ lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */ ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */ cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */ cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */ cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/ ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */ ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ double __ieee754_pow(double x, double y) { double z,ax,z_h,z_l,p_h,p_l; double y1,t1,t2,r,s,t,u,v,w; int32_t i,j,k,yisint,n; int32_t hx,hy,ix,iy; u_int32_t lx,ly; EXTRACT_WORDS(hx,lx,x); EXTRACT_WORDS(hy,ly,y); ix = hx&0x7fffffff; iy = hy&0x7fffffff; /* y==zero: x**0 = 1 */ if((iy|ly)==0) return one; /* x==1: 1**y = 1, even if y is NaN */ if (hx==0x3ff00000 && lx == 0) return one; /* y!=zero: result is NaN if either arg is NaN */ if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) || iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) - return (x+0.0)+(y+0.0); + return nan_mix(x, y); /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer * yisint = 1 ... y is an odd int * yisint = 2 ... y is an even int */ yisint = 0; if(hx<0) { if(iy>=0x43400000) yisint = 2; /* even integer y */ else if(iy>=0x3ff00000) { k = (iy>>20)-0x3ff; /* exponent */ if(k>20) { j = ly>>(52-k); if((j<<(52-k))==ly) yisint = 2-(j&1); } else if(ly==0) { j = iy>>(20-k); if((j<<(20-k))==iy) yisint = 2-(j&1); } } } /* special value of y */ if(ly==0) { if (iy==0x7ff00000) { /* y is +-inf */ if(((ix-0x3ff00000)|lx)==0) return one; /* (-1)**+-inf is 1 */ else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */ return (hy>=0)? y: zero; else /* (|x|<1)**-,+inf = inf,0 */ return (hy<0)?-y: zero; } if(iy==0x3ff00000) { /* y is +-1 */ if(hy<0) return one/x; else return x; } if(hy==0x40000000) return x*x; /* y is 2 */ if(hy==0x3fe00000) { /* y is 0.5 */ if(hx>=0) /* x >= +0 */ return sqrt(x); } } ax = fabs(x); /* special value of x */ if(lx==0) { if(ix==0x7ff00000||ix==0||ix==0x3ff00000){ z = ax; /*x is +-0,+-inf,+-1*/ if(hy<0) z = one/z; /* z = (1/|x|) */ if(hx<0) { if(((ix-0x3ff00000)|yisint)==0) { z = (z-z)/(z-z); /* (-1)**non-int is NaN */ } else if(yisint==1) z = -z; /* (x<0)**odd = -(|x|**odd) */ } return z; } } /* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be n = (hx>>31)+1; but ANSI C says a right shift of a signed negative quantity is implementation defined. */ n = ((u_int32_t)hx>>31)-1; /* (x<0)**(non-int) is NaN */ if((n|yisint)==0) return (x-x)/(x-x); s = one; /* s (sign of result -ve**odd) = -1 else = 1 */ if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */ /* |y| is huge */ if(iy>0x41e00000) { /* if |y| > 2**31 */ if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */ if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny; if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny; } /* over/underflow if x is not close to one */ if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny; if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny; /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ t = ax-one; /* t has 20 trailing zeros */ w = (t*t)*(half-t*(thrd-t*qrtr)); u = ivln2_h*t; /* ivln2_h has 21 sig. bits */ v = t*ivln2_l-w*ivln2; t1 = u+v; SET_LOW_WORD(t1,0); t2 = v-(t1-u); } else { double ss,s2,s_h,s_l,t_h,t_l; n = 0; /* take care subnormal number */ if(ix<0x00100000) {ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); } n += ((ix)>>20)-0x3ff; j = ix&0x000fffff; /* determine interval */ ix = j|0x3ff00000; /* normalize ix */ if(j<=0x3988E) k=0; /* |x|>1)|0x20000000)+0x00080000+(k<<18)); t_l = ax - (t_h-bp[k]); s_l = v*((u-s_h*t_h)-s_h*t_l); /* compute log(ax) */ s2 = ss*ss; r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6))))); r += s_l*(s_h+ss); s2 = s_h*s_h; t_h = 3+s2+r; SET_LOW_WORD(t_h,0); t_l = r-((t_h-3)-s2); /* u+v = ss*(1+...) */ u = s_h*t_h; v = s_l*t_h+t_l*ss; /* 2/(3log2)*(ss+...) */ p_h = u+v; SET_LOW_WORD(p_h,0); p_l = v-(p_h-u); z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ z_l = cp_l*p_h+p_l*cp+dp_l[k]; /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */ t = n; t1 = (((z_h+z_l)+dp_h[k])+t); SET_LOW_WORD(t1,0); t2 = z_l-(((t1-t)-dp_h[k])-z_h); } /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ y1 = y; SET_LOW_WORD(y1,0); p_l = (y-y1)*t1+y*t2; p_h = y1*t1; z = p_l+p_h; EXTRACT_WORDS(j,i,z); if (j>=0x40900000) { /* z >= 1024 */ if(((j-0x40900000)|i)!=0) /* if z > 1024 */ return s*huge*huge; /* overflow */ else { if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */ } } else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */ if(((j-0xc090cc00)|i)!=0) /* z < -1075 */ return s*tiny*tiny; /* underflow */ else { if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */ } } /* * compute 2**(p_h+p_l) */ i = j&0x7fffffff; k = (i>>20)-0x3ff; n = 0; if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */ n = j+(0x00100000>>(k+1)); k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */ t = zero; SET_HIGH_WORD(t,n&~(0x000fffff>>k)); n = ((n&0x000fffff)|0x00100000)>>(20-k); if(j<0) n = -n; p_h -= t; } t = p_l+p_h; SET_LOW_WORD(t,0); u = t*lg2_h; v = (p_l-(t-p_h))*lg2+t*lg2_l; z = u+v; w = v-(z-u); t = z*z; t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); r = (z*t1)/(t1-two)-(w+z*w); z = one-(r-z); GET_HIGH_WORD(j,z); j += (n<<20); if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */ else SET_HIGH_WORD(z,j); return s*z; } #if (LDBL_MANT_DIG == 53) __weak_reference(pow, powl); #endif Index: head/lib/msun/src/e_powf.c =================================================================== --- head/lib/msun/src/e_powf.c (revision 336361) +++ head/lib/msun/src/e_powf.c (revision 336362) @@ -1,252 +1,252 @@ /* e_powf.c -- float version of e_pow.c. * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" static const float bp[] = {1.0, 1.5,}, dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */ dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */ zero = 0.0, half = 0.5, qrtr = 0.25, thrd = 3.33333343e-01, /* 0x3eaaaaab */ one = 1.0, two = 2.0, two24 = 16777216.0, /* 0x4b800000 */ huge = 1.0e30, tiny = 1.0e-30, /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ L1 = 6.0000002384e-01, /* 0x3f19999a */ L2 = 4.2857143283e-01, /* 0x3edb6db7 */ L3 = 3.3333334327e-01, /* 0x3eaaaaab */ L4 = 2.7272811532e-01, /* 0x3e8ba305 */ L5 = 2.3066075146e-01, /* 0x3e6c3255 */ L6 = 2.0697501302e-01, /* 0x3e53f142 */ P1 = 1.6666667163e-01, /* 0x3e2aaaab */ P2 = -2.7777778450e-03, /* 0xbb360b61 */ P3 = 6.6137559770e-05, /* 0x388ab355 */ P4 = -1.6533901999e-06, /* 0xb5ddea0e */ P5 = 4.1381369442e-08, /* 0x3331bb4c */ lg2 = 6.9314718246e-01, /* 0x3f317218 */ lg2_h = 6.93145752e-01, /* 0x3f317200 */ lg2_l = 1.42860654e-06, /* 0x35bfbe8c */ ovt = 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */ cp = 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */ cp_h = 9.6191406250e-01, /* 0x3f764000 =12b cp */ cp_l = -1.1736857402e-04, /* 0xb8f623c6 =tail of cp_h */ ivln2 = 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */ ivln2_h = 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/ ivln2_l = 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/ float __ieee754_powf(float x, float y) { float z,ax,z_h,z_l,p_h,p_l; float y1,t1,t2,r,s,sn,t,u,v,w; int32_t i,j,k,yisint,n; int32_t hx,hy,ix,iy,is; GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hy,y); ix = hx&0x7fffffff; iy = hy&0x7fffffff; /* y==zero: x**0 = 1 */ if(iy==0) return one; /* x==1: 1**y = 1, even if y is NaN */ if (hx==0x3f800000) return one; /* y!=zero: result is NaN if either arg is NaN */ if(ix > 0x7f800000 || iy > 0x7f800000) - return (x+0.0F)+(y+0.0F); + return nan_mix(x, y); /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer * yisint = 1 ... y is an odd int * yisint = 2 ... y is an even int */ yisint = 0; if(hx<0) { if(iy>=0x4b800000) yisint = 2; /* even integer y */ else if(iy>=0x3f800000) { k = (iy>>23)-0x7f; /* exponent */ j = iy>>(23-k); if((j<<(23-k))==iy) yisint = 2-(j&1); } } /* special value of y */ if (iy==0x7f800000) { /* y is +-inf */ if (ix==0x3f800000) return one; /* (-1)**+-inf is NaN */ else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */ return (hy>=0)? y: zero; else /* (|x|<1)**-,+inf = inf,0 */ return (hy<0)?-y: zero; } if(iy==0x3f800000) { /* y is +-1 */ if(hy<0) return one/x; else return x; } if(hy==0x40000000) return x*x; /* y is 2 */ if(hy==0x3f000000) { /* y is 0.5 */ if(hx>=0) /* x >= +0 */ return __ieee754_sqrtf(x); } ax = fabsf(x); /* special value of x */ if(ix==0x7f800000||ix==0||ix==0x3f800000){ z = ax; /*x is +-0,+-inf,+-1*/ if(hy<0) z = one/z; /* z = (1/|x|) */ if(hx<0) { if(((ix-0x3f800000)|yisint)==0) { z = (z-z)/(z-z); /* (-1)**non-int is NaN */ } else if(yisint==1) z = -z; /* (x<0)**odd = -(|x|**odd) */ } return z; } n = ((u_int32_t)hx>>31)-1; /* (x<0)**(non-int) is NaN */ if((n|yisint)==0) return (x-x)/(x-x); sn = one; /* s (sign of result -ve**odd) = -1 else = 1 */ if((n|(yisint-1))==0) sn = -one;/* (-ve)**(odd int) */ /* |y| is huge */ if(iy>0x4d000000) { /* if |y| > 2**27 */ /* over/underflow if x is not close to one */ if(ix<0x3f7ffff8) return (hy<0)? sn*huge*huge:sn*tiny*tiny; if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny; /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ t = ax-1; /* t has 20 trailing zeros */ w = (t*t)*(half-t*(thrd-t*qrtr)); u = ivln2_h*t; /* ivln2_h has 16 sig. bits */ v = t*ivln2_l-w*ivln2; t1 = u+v; GET_FLOAT_WORD(is,t1); SET_FLOAT_WORD(t1,is&0xfffff000); t2 = v-(t1-u); } else { float s2,s_h,s_l,t_h,t_l; n = 0; /* take care subnormal number */ if(ix<0x00800000) {ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); } n += ((ix)>>23)-0x7f; j = ix&0x007fffff; /* determine interval */ ix = j|0x3f800000; /* normalize ix */ if(j<=0x1cc471) k=0; /* |x|>1)&0xfffff000)|0x20000000; SET_FLOAT_WORD(t_h,is+0x00400000+(k<<21)); t_l = ax - (t_h-bp[k]); s_l = v*((u-s_h*t_h)-s_h*t_l); /* compute log(ax) */ s2 = s*s; r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6))))); r += s_l*(s_h+s); s2 = s_h*s_h; t_h = 3+s2+r; GET_FLOAT_WORD(is,t_h); SET_FLOAT_WORD(t_h,is&0xfffff000); t_l = r-((t_h-3)-s2); /* u+v = s*(1+...) */ u = s_h*t_h; v = s_l*t_h+t_l*s; /* 2/(3log2)*(s+...) */ p_h = u+v; GET_FLOAT_WORD(is,p_h); SET_FLOAT_WORD(p_h,is&0xfffff000); p_l = v-(p_h-u); z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ z_l = cp_l*p_h+p_l*cp+dp_l[k]; /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ t = n; t1 = (((z_h+z_l)+dp_h[k])+t); GET_FLOAT_WORD(is,t1); SET_FLOAT_WORD(t1,is&0xfffff000); t2 = z_l-(((t1-t)-dp_h[k])-z_h); } /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ GET_FLOAT_WORD(is,y); SET_FLOAT_WORD(y1,is&0xfffff000); p_l = (y-y1)*t1+y*t2; p_h = y1*t1; z = p_l+p_h; GET_FLOAT_WORD(j,z); if (j>0x43000000) /* if z > 128 */ return sn*huge*huge; /* overflow */ else if (j==0x43000000) { /* if z == 128 */ if(p_l+ovt>z-p_h) return sn*huge*huge; /* overflow */ } else if ((j&0x7fffffff)>0x43160000) /* z <= -150 */ return sn*tiny*tiny; /* underflow */ else if (j==0xc3160000){ /* z == -150 */ if(p_l<=z-p_h) return sn*tiny*tiny; /* underflow */ } /* * compute 2**(p_h+p_l) */ i = j&0x7fffffff; k = (i>>23)-0x7f; n = 0; if(i>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */ n = j+(0x00800000>>(k+1)); k = ((n&0x7fffffff)>>23)-0x7f; /* new k for n */ SET_FLOAT_WORD(t,n&~(0x007fffff>>k)); n = ((n&0x007fffff)|0x00800000)>>(23-k); if(j<0) n = -n; p_h -= t; } t = p_l+p_h; GET_FLOAT_WORD(is,t); SET_FLOAT_WORD(t,is&0xffff8000); u = t*lg2_h; v = (p_l-(t-p_h))*lg2+t*lg2_l; z = u+v; w = v-(z-u); t = z*z; t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); r = (z*t1)/(t1-two)-(w+z*w); z = one-(r-z); GET_FLOAT_WORD(j,z); j += (n<<23); if((j>>23)<=0) z = scalbnf(z,n); /* subnormal output */ else SET_FLOAT_WORD(z,j); return sn*z; } Index: head/lib/msun/src/e_remainder.c =================================================================== --- head/lib/msun/src/e_remainder.c (revision 336361) +++ head/lib/msun/src/e_remainder.c (revision 336362) @@ -1,79 +1,79 @@ /* @(#)e_remainder.c 1.3 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); /* __ieee754_remainder(x,p) * Return : * returns x REM p = x - [x/p]*p as if in infinite * precise arithmetic, where [x/p] is the (infinite bit) * integer nearest x/p (in half way case choose the even one). * Method : * Based on fmod() return x-[x/p]chopped*p exactlp. */ #include #include "math.h" #include "math_private.h" static const double zero = 0.0; double __ieee754_remainder(double x, double p) { int32_t hx,hp; u_int32_t sx,lx,lp; double p_half; EXTRACT_WORDS(hx,lx,x); EXTRACT_WORDS(hp,lp,p); sx = hx&0x80000000; hp &= 0x7fffffff; hx &= 0x7fffffff; /* purge off exception values */ - if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */ - if((hx>=0x7ff00000)|| /* x not finite */ + if(((hp|lp)==0)|| /* p = 0 */ + (hx>=0x7ff00000)|| /* x not finite */ ((hp>=0x7ff00000)&& /* p is NaN */ (((hp-0x7ff00000)|lp)!=0))) - return ((long double)x*p)/((long double)x*p); + return nan_mix(x, p)/nan_mix(x, p); if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */ if (((hx-hp)|(lx-lp))==0) return zero*x; x = fabs(x); p = fabs(p); if (hp<0x00200000) { if(x+x>p) { x-=p; if(x+x>=p) x -= p; } } else { p_half = 0.5*p; if(x>p_half) { x-=p; if(x>=p_half) x -= p; } } GET_HIGH_WORD(hx,x); if ((hx&0x7fffffff)==0) hx = 0; SET_HIGH_WORD(x,hx^sx); return x; } #if LDBL_MANT_DIG == 53 __weak_reference(remainder, remainderl); #endif Index: head/lib/msun/src/e_remainderf.c =================================================================== --- head/lib/msun/src/e_remainderf.c (revision 336361) +++ head/lib/msun/src/e_remainderf.c (revision 336362) @@ -1,65 +1,65 @@ /* e_remainderf.c -- float version of e_remainder.c. * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" static const float zero = 0.0; float __ieee754_remainderf(float x, float p) { int32_t hx,hp; u_int32_t sx; float p_half; GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hp,p); sx = hx&0x80000000; hp &= 0x7fffffff; hx &= 0x7fffffff; /* purge off exception values */ - if(hp==0) return (x*p)/(x*p); /* p = 0 */ - if((hx>=0x7f800000)|| /* x not finite */ + if((hp==0)|| /* p = 0 */ + (hx>=0x7f800000)|| /* x not finite */ ((hp>0x7f800000))) /* p is NaN */ - return ((long double)x*p)/((long double)x*p); + return nan_mix(x, p)/nan_mix(x, p); if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */ if ((hx-hp)==0) return zero*x; x = fabsf(x); p = fabsf(p); if (hp<0x01000000) { if(x+x>p) { x-=p; if(x+x>=p) x -= p; } } else { p_half = (float)0.5*p; if(x>p_half) { x-=p; if(x>=p_half) x -= p; } } GET_FLOAT_WORD(hx,x); if ((hx&0x7fffffff)==0) hx = 0; SET_FLOAT_WORD(x,hx^sx); return x; } Index: head/lib/msun/src/math_private.h =================================================================== --- head/lib/msun/src/math_private.h (revision 336361) +++ head/lib/msun/src/math_private.h (revision 336362) @@ -1,834 +1,852 @@ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* * from: @(#)fdlibm.h 5.1 93/09/24 * $FreeBSD$ */ #ifndef _MATH_PRIVATE_H_ #define _MATH_PRIVATE_H_ #include #include /* * The original fdlibm code used statements like: * n0 = ((*(int*)&one)>>29)^1; * index of high word * * ix0 = *(n0+(int*)&x); * high word of x * * ix1 = *((1-n0)+(int*)&x); * low word of x * * to dig two 32 bit words out of the 64 bit IEEE floating point * value. That is non-ANSI, and, moreover, the gcc instruction * scheduler gets it wrong. We instead use the following macros. * Unlike the original code, we determine the endianness at compile * time, not at run time; I don't see much benefit to selecting * endianness at run time. */ /* * A union which permits us to convert between a double and two 32 bit * ints. */ #ifdef __arm__ #if defined(__VFP_FP__) || defined(__ARM_EABI__) #define IEEE_WORD_ORDER BYTE_ORDER #else #define IEEE_WORD_ORDER BIG_ENDIAN #endif #else /* __arm__ */ #define IEEE_WORD_ORDER BYTE_ORDER #endif /* A union which permits us to convert between a long double and four 32 bit ints. */ #if IEEE_WORD_ORDER == BIG_ENDIAN typedef union { long double value; struct { u_int32_t mswhi; u_int32_t mswlo; u_int32_t lswhi; u_int32_t lswlo; } parts32; struct { u_int64_t msw; u_int64_t lsw; } parts64; } ieee_quad_shape_type; #endif #if IEEE_WORD_ORDER == LITTLE_ENDIAN typedef union { long double value; struct { u_int32_t lswlo; u_int32_t lswhi; u_int32_t mswlo; u_int32_t mswhi; } parts32; struct { u_int64_t lsw; u_int64_t msw; } parts64; } ieee_quad_shape_type; #endif #if IEEE_WORD_ORDER == BIG_ENDIAN typedef union { double value; struct { u_int32_t msw; u_int32_t lsw; } parts; struct { u_int64_t w; } xparts; } ieee_double_shape_type; #endif #if IEEE_WORD_ORDER == LITTLE_ENDIAN typedef union { double value; struct { u_int32_t lsw; u_int32_t msw; } parts; struct { u_int64_t w; } xparts; } ieee_double_shape_type; #endif /* Get two 32 bit ints from a double. */ #define EXTRACT_WORDS(ix0,ix1,d) \ do { \ ieee_double_shape_type ew_u; \ ew_u.value = (d); \ (ix0) = ew_u.parts.msw; \ (ix1) = ew_u.parts.lsw; \ } while (0) /* Get a 64-bit int from a double. */ #define EXTRACT_WORD64(ix,d) \ do { \ ieee_double_shape_type ew_u; \ ew_u.value = (d); \ (ix) = ew_u.xparts.w; \ } while (0) /* Get the more significant 32 bit int from a double. */ #define GET_HIGH_WORD(i,d) \ do { \ ieee_double_shape_type gh_u; \ gh_u.value = (d); \ (i) = gh_u.parts.msw; \ } while (0) /* Get the less significant 32 bit int from a double. */ #define GET_LOW_WORD(i,d) \ do { \ ieee_double_shape_type gl_u; \ gl_u.value = (d); \ (i) = gl_u.parts.lsw; \ } while (0) /* Set a double from two 32 bit ints. */ #define INSERT_WORDS(d,ix0,ix1) \ do { \ ieee_double_shape_type iw_u; \ iw_u.parts.msw = (ix0); \ iw_u.parts.lsw = (ix1); \ (d) = iw_u.value; \ } while (0) /* Set a double from a 64-bit int. */ #define INSERT_WORD64(d,ix) \ do { \ ieee_double_shape_type iw_u; \ iw_u.xparts.w = (ix); \ (d) = iw_u.value; \ } while (0) /* Set the more significant 32 bits of a double from an int. */ #define SET_HIGH_WORD(d,v) \ do { \ ieee_double_shape_type sh_u; \ sh_u.value = (d); \ sh_u.parts.msw = (v); \ (d) = sh_u.value; \ } while (0) /* Set the less significant 32 bits of a double from an int. */ #define SET_LOW_WORD(d,v) \ do { \ ieee_double_shape_type sl_u; \ sl_u.value = (d); \ sl_u.parts.lsw = (v); \ (d) = sl_u.value; \ } while (0) /* * A union which permits us to convert between a float and a 32 bit * int. */ typedef union { float value; /* FIXME: Assumes 32 bit int. */ unsigned int word; } ieee_float_shape_type; /* Get a 32 bit int from a float. */ #define GET_FLOAT_WORD(i,d) \ do { \ ieee_float_shape_type gf_u; \ gf_u.value = (d); \ (i) = gf_u.word; \ } while (0) /* Set a float from a 32 bit int. */ #define SET_FLOAT_WORD(d,i) \ do { \ ieee_float_shape_type sf_u; \ sf_u.word = (i); \ (d) = sf_u.value; \ } while (0) /* * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long * double. */ #define EXTRACT_LDBL80_WORDS(ix0,ix1,d) \ do { \ union IEEEl2bits ew_u; \ ew_u.e = (d); \ (ix0) = ew_u.xbits.expsign; \ (ix1) = ew_u.xbits.man; \ } while (0) /* * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit * long double. */ #define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \ do { \ union IEEEl2bits ew_u; \ ew_u.e = (d); \ (ix0) = ew_u.xbits.expsign; \ (ix1) = ew_u.xbits.manh; \ (ix2) = ew_u.xbits.manl; \ } while (0) /* Get expsign as a 16 bit int from a long double. */ #define GET_LDBL_EXPSIGN(i,d) \ do { \ union IEEEl2bits ge_u; \ ge_u.e = (d); \ (i) = ge_u.xbits.expsign; \ } while (0) /* * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int * mantissa. */ #define INSERT_LDBL80_WORDS(d,ix0,ix1) \ do { \ union IEEEl2bits iw_u; \ iw_u.xbits.expsign = (ix0); \ iw_u.xbits.man = (ix1); \ (d) = iw_u.e; \ } while (0) /* * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints * comprising the mantissa. */ #define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2) \ do { \ union IEEEl2bits iw_u; \ iw_u.xbits.expsign = (ix0); \ iw_u.xbits.manh = (ix1); \ iw_u.xbits.manl = (ix2); \ (d) = iw_u.e; \ } while (0) /* Set expsign of a long double from a 16 bit int. */ #define SET_LDBL_EXPSIGN(d,v) \ do { \ union IEEEl2bits se_u; \ se_u.e = (d); \ se_u.xbits.expsign = (v); \ (d) = se_u.e; \ } while (0) #ifdef __i386__ /* Long double constants are broken on i386. */ #define LD80C(m, ex, v) { \ .xbits.man = __CONCAT(m, ULL), \ .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0), \ } #else /* The above works on non-i386 too, but we use this to check v. */ #define LD80C(m, ex, v) { .e = (v), } #endif #ifdef FLT_EVAL_METHOD /* * Attempt to get strict C99 semantics for assignment with non-C99 compilers. */ #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0 #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval)) #else #define STRICT_ASSIGN(type, lval, rval) do { \ volatile type __lval; \ \ if (sizeof(type) >= sizeof(long double)) \ (lval) = (rval); \ else { \ __lval = (rval); \ (lval) = __lval; \ } \ } while (0) #endif #endif /* FLT_EVAL_METHOD */ /* Support switching the mode to FP_PE if necessary. */ #if defined(__i386__) && !defined(NO_FPSETPREC) #define ENTERI() ENTERIT(long double) #define ENTERIT(returntype) \ returntype __retval; \ fp_prec_t __oprec; \ \ if ((__oprec = fpgetprec()) != FP_PE) \ fpsetprec(FP_PE) #define RETURNI(x) do { \ __retval = (x); \ if (__oprec != FP_PE) \ fpsetprec(__oprec); \ RETURNF(__retval); \ } while (0) #define ENTERV() \ fp_prec_t __oprec; \ \ if ((__oprec = fpgetprec()) != FP_PE) \ fpsetprec(FP_PE) #define RETURNV() do { \ if (__oprec != FP_PE) \ fpsetprec(__oprec); \ return; \ } while (0) #else #define ENTERI() #define ENTERIT(x) #define RETURNI(x) RETURNF(x) #define ENTERV() #define RETURNV() return #endif /* Default return statement if hack*_t() is not used. */ #define RETURNF(v) return (v) /* * 2sum gives the same result as 2sumF without requiring |a| >= |b| or * a == 0, but is slower. */ #define _2sum(a, b) do { \ __typeof(a) __s, __w; \ \ __w = (a) + (b); \ __s = __w - (a); \ (b) = ((a) - (__w - __s)) + ((b) - __s); \ (a) = __w; \ } while (0) /* * 2sumF algorithm. * * "Normalize" the terms in the infinite-precision expression a + b for * the sum of 2 floating point values so that b is as small as possible * relative to 'a'. (The resulting 'a' is the value of the expression in * the same precision as 'a' and the resulting b is the rounding error.) * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and * exponent overflow or underflow must not occur. This uses a Theorem of * Dekker (1971). See Knuth (1981) 4.2.2 Theorem C. The name "TwoSum" * is apparently due to Skewchuk (1997). * * For this to always work, assignment of a + b to 'a' must not retain any * extra precision in a + b. This is required by C standards but broken * in many compilers. The brokenness cannot be worked around using * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this * algorithm would be destroyed by non-null strict assignments. (The * compilers are correct to be broken -- the efficiency of all floating * point code calculations would be destroyed similarly if they forced the * conversions.) * * Fortunately, a case that works well can usually be arranged by building * any extra precision into the type of 'a' -- 'a' should have type float_t, * double_t or long double. b's type should be no larger than 'a's type. * Callers should use these types with scopes as large as possible, to * reduce their own extra-precision and efficiciency problems. In * particular, they shouldn't convert back and forth just to call here. */ #ifdef DEBUG #define _2sumF(a, b) do { \ __typeof(a) __w; \ volatile __typeof(a) __ia, __ib, __r, __vw; \ \ __ia = (a); \ __ib = (b); \ assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib)); \ \ __w = (a) + (b); \ (b) = ((a) - __w) + (b); \ (a) = __w; \ \ /* The next 2 assertions are weak if (a) is already long double. */ \ assert((long double)__ia + __ib == (long double)(a) + (b)); \ __vw = __ia + __ib; \ __r = __ia - __vw; \ __r += __ib; \ assert(__vw == (a) && __r == (b)); \ } while (0) #else /* !DEBUG */ #define _2sumF(a, b) do { \ __typeof(a) __w; \ \ __w = (a) + (b); \ (b) = ((a) - __w) + (b); \ (a) = __w; \ } while (0) #endif /* DEBUG */ /* * Set x += c, where x is represented in extra precision as a + b. * x must be sufficiently normalized and sufficiently larger than c, * and the result is then sufficiently normalized. * * The details of ordering are that |a| must be >= |c| (so that (a, c) * can be normalized without extra work to swap 'a' with c). The details of * the normalization are that b must be small relative to the normalized 'a'. * Normalization of (a, c) makes the normalized c tiny relative to the * normalized a, so b remains small relative to 'a' in the result. However, * b need not ever be tiny relative to 'a'. For example, b might be about * 2**20 times smaller than 'a' to give about 20 extra bits of precision. * That is usually enough, and adding c (which by normalization is about * 2**53 times smaller than a) cannot change b significantly. However, * cancellation of 'a' with c in normalization of (a, c) may reduce 'a' * significantly relative to b. The caller must ensure that significant * cancellation doesn't occur, either by having c of the same sign as 'a', * or by having |c| a few percent smaller than |a|. Pre-normalization of * (a, b) may help. * * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2 * exercise 19). We gain considerable efficiency by requiring the terms to * be sufficiently normalized and sufficiently increasing. */ #define _3sumF(a, b, c) do { \ __typeof(a) __tmp; \ \ __tmp = (c); \ _2sumF(__tmp, (a)); \ (b) += (a); \ (a) = __tmp; \ } while (0) /* * Common routine to process the arguments to nan(), nanf(), and nanl(). */ void _scan_nan(uint32_t *__words, int __num_words, const char *__s); +/* + * Mix 1 or 2 NaNs. First add 0 to each arg. This normally just turns + * signaling NaNs into quiet NaNs by setting a quiet bit. We do this + * because we want to never return a signaling NaN, and also because we + * don't want the quiet bit to affect the result. Then mix the converted + * args using addition. The result is typically the arg whose mantissa + * bits (considered as in integer) are largest. + * + * Technical complications: the result in bits might depend on the precision + * and/or on compiler optimizations, especially when different register sets + * are used for different precisions. Try to make the result not depend on + * at least the precision by always doing the main mixing step in long double + * precision. Try to reduce dependencies on optimizations by adding the + * the 0's in different precisions (unless everything is in long double + * precision). + */ +#define nan_mix(x, y) (((x) + 0.0L) + ((y) + 0)) + #ifdef _COMPLEX_H /* * C99 specifies that complex numbers have the same representation as * an array of two elements, where the first element is the real part * and the second element is the imaginary part. */ typedef union { float complex f; float a[2]; } float_complex; typedef union { double complex f; double a[2]; } double_complex; typedef union { long double complex f; long double a[2]; } long_double_complex; #define REALPART(z) ((z).a[0]) #define IMAGPART(z) ((z).a[1]) /* * Inline functions that can be used to construct complex values. * * The C99 standard intends x+I*y to be used for this, but x+I*y is * currently unusable in general since gcc introduces many overflow, * underflow, sign and efficiency bugs by rewriting I*y as * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product. * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted * to -0.0+I*0.0. * * The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL() * to construct complex values. Compilers that conform to the C99 * standard require the following functions to avoid the above issues. */ #ifndef CMPLXF static __inline float complex CMPLXF(float x, float y) { float_complex z; REALPART(z) = x; IMAGPART(z) = y; return (z.f); } #endif #ifndef CMPLX static __inline double complex CMPLX(double x, double y) { double_complex z; REALPART(z) = x; IMAGPART(z) = y; return (z.f); } #endif #ifndef CMPLXL static __inline long double complex CMPLXL(long double x, long double y) { long_double_complex z; REALPART(z) = x; IMAGPART(z) = y; return (z.f); } #endif #endif /* _COMPLEX_H */ #ifdef __GNUCLIKE_ASM /* Asm versions of some functions. */ #ifdef __amd64__ static __inline int irint(double x) { int n; asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x)); return (n); } #define HAVE_EFFICIENT_IRINT #endif #ifdef __i386__ static __inline int irint(double x) { int n; asm("fistl %0" : "=m" (n) : "t" (x)); return (n); } #define HAVE_EFFICIENT_IRINT #endif #if defined(__amd64__) || defined(__i386__) static __inline int irintl(long double x) { int n; asm("fistl %0" : "=m" (n) : "t" (x)); return (n); } #define HAVE_EFFICIENT_IRINTL #endif #endif /* __GNUCLIKE_ASM */ #ifdef DEBUG #if defined(__amd64__) || defined(__i386__) #define breakpoint() asm("int $3") #else #include #define breakpoint() raise(SIGTRAP) #endif #endif /* Write a pari script to test things externally. */ #ifdef DOPRINT #include #ifndef DOPRINT_SWIZZLE #define DOPRINT_SWIZZLE 0 #endif #ifdef DOPRINT_LD80 #define DOPRINT_START(xp) do { \ uint64_t __lx; \ uint16_t __hx; \ \ /* Hack to give more-problematic args. */ \ EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \ __lx ^= DOPRINT_SWIZZLE; \ INSERT_LDBL80_WORDS(*xp, __hx, __lx); \ printf("x = %.21Lg; ", (long double)*xp); \ } while (0) #define DOPRINT_END1(v) \ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) #define DOPRINT_END2(hi, lo) \ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ (long double)(hi), (long double)(lo)) #elif defined(DOPRINT_D64) #define DOPRINT_START(xp) do { \ uint32_t __hx, __lx; \ \ EXTRACT_WORDS(__hx, __lx, *xp); \ __lx ^= DOPRINT_SWIZZLE; \ INSERT_WORDS(*xp, __hx, __lx); \ printf("x = %.21Lg; ", (long double)*xp); \ } while (0) #define DOPRINT_END1(v) \ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) #define DOPRINT_END2(hi, lo) \ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ (long double)(hi), (long double)(lo)) #elif defined(DOPRINT_F32) #define DOPRINT_START(xp) do { \ uint32_t __hx; \ \ GET_FLOAT_WORD(__hx, *xp); \ __hx ^= DOPRINT_SWIZZLE; \ SET_FLOAT_WORD(*xp, __hx); \ printf("x = %.21Lg; ", (long double)*xp); \ } while (0) #define DOPRINT_END1(v) \ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) #define DOPRINT_END2(hi, lo) \ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ (long double)(hi), (long double)(lo)) #else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */ #ifndef DOPRINT_SWIZZLE_HIGH #define DOPRINT_SWIZZLE_HIGH 0 #endif #define DOPRINT_START(xp) do { \ uint64_t __lx, __llx; \ uint16_t __hx; \ \ EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \ __llx ^= DOPRINT_SWIZZLE; \ __lx ^= DOPRINT_SWIZZLE_HIGH; \ INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \ printf("x = %.36Lg; ", (long double)*xp); \ } while (0) #define DOPRINT_END1(v) \ printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v)) #define DOPRINT_END2(hi, lo) \ printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \ (long double)(hi), (long double)(lo)) #endif /* DOPRINT_LD80 */ #else /* !DOPRINT */ #define DOPRINT_START(xp) #define DOPRINT_END1(v) #define DOPRINT_END2(hi, lo) #endif /* DOPRINT */ #define RETURNP(x) do { \ DOPRINT_END1(x); \ RETURNF(x); \ } while (0) #define RETURNPI(x) do { \ DOPRINT_END1(x); \ RETURNI(x); \ } while (0) #define RETURN2P(x, y) do { \ DOPRINT_END2((x), (y)); \ RETURNF((x) + (y)); \ } while (0) #define RETURN2PI(x, y) do { \ DOPRINT_END2((x), (y)); \ RETURNI((x) + (y)); \ } while (0) #ifdef STRUCT_RETURN #define RETURNSP(rp) do { \ if (!(rp)->lo_set) \ RETURNP((rp)->hi); \ RETURN2P((rp)->hi, (rp)->lo); \ } while (0) #define RETURNSPI(rp) do { \ if (!(rp)->lo_set) \ RETURNPI((rp)->hi); \ RETURN2PI((rp)->hi, (rp)->lo); \ } while (0) #endif #define SUM2P(x, y) ({ \ const __typeof (x) __x = (x); \ const __typeof (y) __y = (y); \ \ DOPRINT_END2(__x, __y); \ __x + __y; \ }) /* * ieee style elementary functions * * We rename functions here to improve other sources' diffability * against fdlibm. */ #define __ieee754_sqrt sqrt #define __ieee754_acos acos #define __ieee754_acosh acosh #define __ieee754_log log #define __ieee754_log2 log2 #define __ieee754_atanh atanh #define __ieee754_asin asin #define __ieee754_atan2 atan2 #define __ieee754_exp exp #define __ieee754_cosh cosh #define __ieee754_fmod fmod #define __ieee754_pow pow #define __ieee754_lgamma lgamma #define __ieee754_gamma gamma #define __ieee754_lgamma_r lgamma_r #define __ieee754_gamma_r gamma_r #define __ieee754_log10 log10 #define __ieee754_sinh sinh #define __ieee754_hypot hypot #define __ieee754_j0 j0 #define __ieee754_j1 j1 #define __ieee754_y0 y0 #define __ieee754_y1 y1 #define __ieee754_jn jn #define __ieee754_yn yn #define __ieee754_remainder remainder #define __ieee754_scalb scalb #define __ieee754_sqrtf sqrtf #define __ieee754_acosf acosf #define __ieee754_acoshf acoshf #define __ieee754_logf logf #define __ieee754_atanhf atanhf #define __ieee754_asinf asinf #define __ieee754_atan2f atan2f #define __ieee754_expf expf #define __ieee754_coshf coshf #define __ieee754_fmodf fmodf #define __ieee754_powf powf #define __ieee754_lgammaf lgammaf #define __ieee754_gammaf gammaf #define __ieee754_lgammaf_r lgammaf_r #define __ieee754_gammaf_r gammaf_r #define __ieee754_log10f log10f #define __ieee754_log2f log2f #define __ieee754_sinhf sinhf #define __ieee754_hypotf hypotf #define __ieee754_j0f j0f #define __ieee754_j1f j1f #define __ieee754_y0f y0f #define __ieee754_y1f y1f #define __ieee754_jnf jnf #define __ieee754_ynf ynf #define __ieee754_remainderf remainderf #define __ieee754_scalbf scalbf /* fdlibm kernel function */ int __kernel_rem_pio2(double*,double*,int,int,int); /* double precision kernel functions */ #ifndef INLINE_REM_PIO2 int __ieee754_rem_pio2(double,double*); #endif double __kernel_sin(double,double,int); double __kernel_cos(double,double); double __kernel_tan(double,double,int); double __ldexp_exp(double,int); #ifdef _COMPLEX_H double complex __ldexp_cexp(double complex,int); #endif /* float precision kernel functions */ #ifndef INLINE_REM_PIO2F int __ieee754_rem_pio2f(float,double*); #endif #ifndef INLINE_KERNEL_SINDF float __kernel_sindf(double); #endif #ifndef INLINE_KERNEL_COSDF float __kernel_cosdf(double); #endif #ifndef INLINE_KERNEL_TANDF float __kernel_tandf(double,int); #endif float __ldexp_expf(float,int); #ifdef _COMPLEX_H float complex __ldexp_cexpf(float complex,int); #endif /* long double precision kernel functions */ long double __kernel_sinl(long double, long double, int); long double __kernel_cosl(long double, long double); long double __kernel_tanl(long double, long double, int); long double __p1evll(long double, void *, int); long double __polevll(long double, void *, int); #endif /* !_MATH_PRIVATE_H_ */ Index: head/lib/msun/src/s_ccosh.c =================================================================== --- head/lib/msun/src/s_ccosh.c (revision 336361) +++ head/lib/msun/src/s_ccosh.c (revision 336362) @@ -1,158 +1,159 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Hyperbolic cosine of a complex argument z = x + i y. * * cosh(z) = cosh(x+iy) * = cosh(x) cos(y) + i sinh(x) sin(y). * * Exceptional values are noted in the comments within the source code. * These values and the return value were taken from n1124.pdf. * The sign of the result for some exceptional values is unspecified but * must satisfy both cosh(conj(z)) == conj(cosh(z)) and cosh(-z) == cosh(z). */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" static const double huge = 0x1p1023; double complex ccosh(double complex z) { double x, y, h; int32_t hx, hy, ix, iy, lx, ly; x = creal(z); y = cimag(z); EXTRACT_WORDS(hx, lx, x); EXTRACT_WORDS(hy, ly, y); ix = 0x7fffffff & hx; iy = 0x7fffffff & hy; /* Handle the nearly-non-exceptional cases where x and y are finite. */ if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) return (CMPLX(cosh(x), x * y)); if (ix < 0x40360000) /* |x| < 22: normal case */ return (CMPLX(cosh(x) * cos(y), sinh(x) * sin(y))); /* |x| >= 22, so cosh(x) ~= exp(|x|) */ if (ix < 0x40862e42) { /* x < 710: exp(|x|) won't overflow */ h = exp(fabs(x)) * 0.5; return (CMPLX(h * cos(y), copysign(h, x) * sin(y))); } else if (ix < 0x4096bbaa) { /* x < 1455: scale to avoid overflow */ z = __ldexp_cexp(CMPLX(fabs(x), y), -1); return (CMPLX(creal(z), cimag(z) * copysign(1, x))); } else { /* x >= 1455: the result always overflows */ h = huge * x; return (CMPLX(h * h * cos(y), h * sin(y))); } } /* * cosh(+-0 +- I Inf) = dNaN + I (+-)(+-)0. * The sign of 0 in the result is unspecified. Choice = product * of the signs of the argument. Raise the invalid floating-point * exception. * * cosh(+-0 +- I NaN) = d(NaN) + I (+-)(+-)0. * The sign of 0 in the result is unspecified. Choice = product * of the signs of the argument. */ if ((ix | lx) == 0) /* && iy >= 0x7ff00000 */ return (CMPLX(y - y, x * copysign(0, y))); /* * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0. * * cosh(NaN +- I 0) = d(NaN) + I (+-)(+-)0. * The sign of 0 in the result is unspecified. Choice = product * of the signs of the argument. */ if ((iy | ly) == 0) /* && ix >= 0x7ff00000 */ return (CMPLX(x * x, copysign(0, x) * y)); /* * cosh(x +- I Inf) = dNaN + I dNaN. * Raise the invalid floating-point exception for finite nonzero x. * * cosh(x + I NaN) = d(NaN) + I d(NaN). * Optionally raises the invalid floating-point exception for finite * nonzero x. Choice = don't raise (except for signaling NaNs). */ if (ix < 0x7ff00000) /* && iy >= 0x7ff00000 */ return (CMPLX(y - y, x * (y - y))); /* * cosh(+-Inf + I NaN) = +Inf + I d(NaN). * * cosh(+-Inf +- I Inf) = +Inf + I dNaN. * The sign of Inf in the result is unspecified. Choice = always +. * Raise the invalid floating-point exception. * * cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y) */ if (ix == 0x7ff00000 && lx == 0) { if (iy >= 0x7ff00000) return (CMPLX(INFINITY, x * (y - y))); return (CMPLX(INFINITY * cos(y), x * sin(y))); } /* * cosh(NaN + I NaN) = d(NaN) + I d(NaN). * * cosh(NaN +- I Inf) = d(NaN) + I d(NaN). * Optionally raises the invalid floating-point exception. * Choice = raise. * * cosh(NaN + I y) = d(NaN) + I d(NaN). * Optionally raises the invalid floating-point exception for finite * nonzero y. Choice = don't raise (except for signaling NaNs). */ - return (CMPLX((x * x) * (y - y), (x + x) * (y - y))); + return (CMPLX(((long double)x * x) * (y - y), + ((long double)x + x) * (y - y))); } double complex ccos(double complex z) { /* ccos(z) = ccosh(I * z) */ return (ccosh(CMPLX(-cimag(z), creal(z)))); } Index: head/lib/msun/src/s_ccoshf.c =================================================================== --- head/lib/msun/src/s_ccoshf.c (revision 336361) +++ head/lib/msun/src/s_ccoshf.c (revision 336362) @@ -1,103 +1,104 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Float version of ccosh(). See s_ccosh.c for details. */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" static const float huge = 0x1p127; float complex ccoshf(float complex z) { float x, y, h; int32_t hx, hy, ix, iy; x = crealf(z); y = cimagf(z); GET_FLOAT_WORD(hx, x); GET_FLOAT_WORD(hy, y); ix = 0x7fffffff & hx; iy = 0x7fffffff & hy; if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) return (CMPLXF(coshf(x), x * y)); if (ix < 0x41100000) /* |x| < 9: normal case */ return (CMPLXF(coshf(x) * cosf(y), sinhf(x) * sinf(y))); /* |x| >= 9, so cosh(x) ~= exp(|x|) */ if (ix < 0x42b17218) { /* x < 88.7: expf(|x|) won't overflow */ h = expf(fabsf(x)) * 0.5F; return (CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y))); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1); return (CMPLXF(crealf(z), cimagf(z) * copysignf(1, x))); } else { /* x >= 192.7: the result always overflows */ h = huge * x; return (CMPLXF(h * h * cosf(y), h * sinf(y))); } } if (ix == 0) /* && iy >= 0x7f800000 */ return (CMPLXF(y - y, x * copysignf(0, y))); if (iy == 0) /* && ix >= 0x7f800000 */ return (CMPLXF(x * x, copysignf(0, x) * y)); if (ix < 0x7f800000) /* && iy >= 0x7f800000 */ return (CMPLXF(y - y, x * (y - y))); if (ix == 0x7f800000) { if (iy >= 0x7f800000) return (CMPLXF(INFINITY, x * (y - y))); return (CMPLXF(INFINITY * cosf(y), x * sinf(y))); } - return (CMPLXF((x * x) * (y - y), (x + x) * (y - y))); + return (CMPLXF(((long double)x * x) * (y - y), + ((long double)x + x) * (y - y))); } float complex ccosf(float complex z) { return (ccoshf(CMPLXF(-cimagf(z), crealf(z)))); } Index: head/lib/msun/src/s_csinh.c =================================================================== --- head/lib/msun/src/s_csinh.c (revision 336361) +++ head/lib/msun/src/s_csinh.c (revision 336362) @@ -1,158 +1,159 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Hyperbolic sine of a complex argument z = x + i y. * * sinh(z) = sinh(x+iy) * = sinh(x) cos(y) + i cosh(x) sin(y). * * Exceptional values are noted in the comments within the source code. * These values and the return value were taken from n1124.pdf. * The sign of the result for some exceptional values is unspecified but * must satisfy both sinh(conj(z)) == conj(sinh(z)) and sinh(-z) == -sinh(z). */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" static const double huge = 0x1p1023; double complex csinh(double complex z) { double x, y, h; int32_t hx, hy, ix, iy, lx, ly; x = creal(z); y = cimag(z); EXTRACT_WORDS(hx, lx, x); EXTRACT_WORDS(hy, ly, y); ix = 0x7fffffff & hx; iy = 0x7fffffff & hy; /* Handle the nearly-non-exceptional cases where x and y are finite. */ if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) return (CMPLX(sinh(x), y)); if (ix < 0x40360000) /* |x| < 22: normal case */ return (CMPLX(sinh(x) * cos(y), cosh(x) * sin(y))); /* |x| >= 22, so cosh(x) ~= exp(|x|) */ if (ix < 0x40862e42) { /* x < 710: exp(|x|) won't overflow */ h = exp(fabs(x)) * 0.5; return (CMPLX(copysign(h, x) * cos(y), h * sin(y))); } else if (ix < 0x4096bbaa) { /* x < 1455: scale to avoid overflow */ z = __ldexp_cexp(CMPLX(fabs(x), y), -1); return (CMPLX(creal(z) * copysign(1, x), cimag(z))); } else { /* x >= 1455: the result always overflows */ h = huge * x; return (CMPLX(h * cos(y), h * h * sin(y))); } } /* * sinh(+-0 +- I Inf) = +-0 + I dNaN. * The sign of 0 in the result is unspecified. Choice = same sign * as the argument. Raise the invalid floating-point exception. * * sinh(+-0 +- I NaN) = +-0 + I d(NaN). * The sign of 0 in the result is unspecified. Choice = same sign * as the argument. */ if ((ix | lx) == 0) /* && iy >= 0x7ff00000 */ return (CMPLX(x, y - y)); /* * sinh(+-Inf +- I 0) = +-Inf + I +-0. * * sinh(NaN +- I 0) = d(NaN) + I +-0. */ if ((iy | ly) == 0) /* && ix >= 0x7ff00000 */ return (CMPLX(x + x, y)); /* * sinh(x +- I Inf) = dNaN + I dNaN. * Raise the invalid floating-point exception for finite nonzero x. * * sinh(x + I NaN) = d(NaN) + I d(NaN). * Optionally raises the invalid floating-point exception for finite * nonzero x. Choice = don't raise (except for signaling NaNs). */ if (ix < 0x7ff00000) /* && iy >= 0x7ff00000 */ return (CMPLX(y - y, y - y)); /* * sinh(+-Inf + I NaN) = +-Inf + I d(NaN). * The sign of Inf in the result is unspecified. Choice = same sign * as the argument. * * sinh(+-Inf +- I Inf) = +-Inf + I dNaN. * The sign of Inf in the result is unspecified. Choice = same sign * as the argument. Raise the invalid floating-point exception. * * sinh(+-Inf + I y) = +-Inf cos(y) + I Inf sin(y) */ if (ix == 0x7ff00000 && lx == 0) { if (iy >= 0x7ff00000) return (CMPLX(x, y - y)); return (CMPLX(x * cos(y), INFINITY * sin(y))); } /* * sinh(NaN1 + I NaN2) = d(NaN1, NaN2) + I d(NaN1, NaN2). * * sinh(NaN +- I Inf) = d(NaN, dNaN) + I d(NaN, dNaN). * Optionally raises the invalid floating-point exception. * Choice = raise. * * sinh(NaN + I y) = d(NaN) + I d(NaN). * Optionally raises the invalid floating-point exception for finite * nonzero y. Choice = don't raise (except for signaling NaNs). */ - return (CMPLX((x + x) * (y - y), (x * x) * (y - y))); + return (CMPLX(((long double)x + x) * (y - y), + ((long double)x * x) * (y - y))); } double complex csin(double complex z) { /* csin(z) = -I * csinh(I * z) = I * conj(csinh(I * conj(z))). */ z = csinh(CMPLX(cimag(z), creal(z))); return (CMPLX(cimag(z), creal(z))); } Index: head/lib/msun/src/s_csinhf.c =================================================================== --- head/lib/msun/src/s_csinhf.c (revision 336361) +++ head/lib/msun/src/s_csinhf.c (revision 336362) @@ -1,104 +1,105 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2005 Bruce D. Evans and Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Float version of csinh(). See s_csinh.c for details. */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" static const float huge = 0x1p127; float complex csinhf(float complex z) { float x, y, h; int32_t hx, hy, ix, iy; x = crealf(z); y = cimagf(z); GET_FLOAT_WORD(hx, x); GET_FLOAT_WORD(hy, y); ix = 0x7fffffff & hx; iy = 0x7fffffff & hy; if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) return (CMPLXF(sinhf(x), y)); if (ix < 0x41100000) /* |x| < 9: normal case */ return (CMPLXF(sinhf(x) * cosf(y), coshf(x) * sinf(y))); /* |x| >= 9, so cosh(x) ~= exp(|x|) */ if (ix < 0x42b17218) { /* x < 88.7: expf(|x|) won't overflow */ h = expf(fabsf(x)) * 0.5F; return (CMPLXF(copysignf(h, x) * cosf(y), h * sinf(y))); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1); return (CMPLXF(crealf(z) * copysignf(1, x), cimagf(z))); } else { /* x >= 192.7: the result always overflows */ h = huge * x; return (CMPLXF(h * cosf(y), h * h * sinf(y))); } } if (ix == 0) /* && iy >= 0x7f800000 */ return (CMPLXF(x, y - y)); if (iy == 0) /* && ix >= 0x7f800000 */ return (CMPLXF(x + x, y)); if (ix < 0x7f800000) /* && iy >= 0x7f800000 */ return (CMPLXF(y - y, y - y)); if (ix == 0x7f800000) { if (iy >= 0x7f800000) return (CMPLXF(x, y - y)); return (CMPLXF(x * cosf(y), INFINITY * sinf(y))); } - return (CMPLXF((x + x) * (y - y), (x * x) * (y - y))); + return (CMPLXF(((long double)x + x) * (y - y), + ((long double)x * x) * (y - y))); } float complex csinf(float complex z) { z = csinhf(CMPLXF(cimagf(z), crealf(z))); return (CMPLXF(cimagf(z), crealf(z))); } Index: head/lib/msun/src/s_csqrt.c =================================================================== --- head/lib/msun/src/s_csqrt.c (revision 336361) +++ head/lib/msun/src/s_csqrt.c (revision 336362) @@ -1,114 +1,114 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2007 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ #include __FBSDID("$FreeBSD$"); #include #include #include #include "math_private.h" /* * gcc doesn't implement complex multiplication or division correctly, * so we need to handle infinities specially. We turn on this pragma to * notify conforming c99 compilers that the fast-but-incorrect code that * gcc generates is acceptable, since the special cases have already been * handled. */ #pragma STDC CX_LIMITED_RANGE ON /* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */ #define THRESH 0x1.a827999fcef32p+1022 double complex csqrt(double complex z) { double complex result; double a, b; double t; int scale; a = creal(z); b = cimag(z); /* Handle special cases. */ if (z == 0) return (CMPLX(0, b)); if (isinf(b)) return (CMPLX(INFINITY, b)); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ - return (CMPLX(a, t)); /* return NaN + NaN i */ + return (CMPLX(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */ } if (isinf(a)) { /* * csqrt(inf + NaN i) = inf + NaN i * csqrt(inf + y i) = inf + 0 i * csqrt(-inf + NaN i) = NaN +- inf i * csqrt(-inf + y i) = 0 + inf i */ if (signbit(a)) return (CMPLX(fabs(b - b), copysign(a, b))); else return (CMPLX(a, copysign(b - b, b))); } - /* - * The remaining special case (b is NaN) is handled just fine by - * the normal code path below. - */ + if (isnan(b)) { + t = (a - a) / (a - a); /* raise invalid */ + return (CMPLX(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */ + } /* Scale to avoid overflow. */ if (fabs(a) >= THRESH || fabs(b) >= THRESH) { a *= 0.25; b *= 0.25; scale = 1; } else { scale = 0; } /* Algorithm 312, CACM vol 10, Oct 1967. */ if (a >= 0) { t = sqrt((a + hypot(a, b)) * 0.5); result = CMPLX(t, b / (2 * t)); } else { t = sqrt((-a + hypot(a, b)) * 0.5); result = CMPLX(fabs(b) / (2 * t), copysign(t, b)); } /* Rescale. */ if (scale) return (result * 2); else return (result); } #if LDBL_MANT_DIG == 53 __weak_reference(csqrt, csqrtl); #endif Index: head/lib/msun/src/s_csqrtf.c =================================================================== --- head/lib/msun/src/s_csqrtf.c (revision 336361) +++ head/lib/msun/src/s_csqrtf.c (revision 336362) @@ -1,90 +1,90 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2007 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" /* * gcc doesn't implement complex multiplication or division correctly, * so we need to handle infinities specially. We turn on this pragma to * notify conforming c99 compilers that the fast-but-incorrect code that * gcc generates is acceptable, since the special cases have already been * handled. */ #pragma STDC CX_LIMITED_RANGE ON float complex csqrtf(float complex z) { float a = crealf(z), b = cimagf(z); double t; /* Handle special cases. */ if (z == 0) return (CMPLXF(0, b)); if (isinf(b)) return (CMPLXF(INFINITY, b)); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ - return (CMPLXF(a, t)); /* return NaN + NaN i */ + return (CMPLXF(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */ } if (isinf(a)) { /* * csqrtf(inf + NaN i) = inf + NaN i * csqrtf(inf + y i) = inf + 0 i * csqrtf(-inf + NaN i) = NaN +- inf i * csqrtf(-inf + y i) = 0 + inf i */ if (signbit(a)) return (CMPLXF(fabsf(b - b), copysignf(a, b))); else return (CMPLXF(a, copysignf(b - b, b))); } - /* - * The remaining special case (b is NaN) is handled just fine by - * the normal code path below. - */ + if (isnan(b)) { + t = (a - a) / (a - a); /* raise invalid */ + return (CMPLXF(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */ + } /* * We compute t in double precision to avoid overflow and to * provide correct rounding in nearly all cases. * This is Algorithm 312, CACM vol 10, Oct 1967. */ if (a >= 0) { t = sqrt((a + hypot(a, b)) * 0.5); return (CMPLXF(t, b / (2.0 * t))); } else { t = sqrt((-a + hypot(a, b)) * 0.5); return (CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b))); } } Index: head/lib/msun/src/s_csqrtl.c =================================================================== --- head/lib/msun/src/s_csqrtl.c (revision 336361) +++ head/lib/msun/src/s_csqrtl.c (revision 336362) @@ -1,118 +1,118 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2007-2008 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ #include __FBSDID("$FreeBSD$"); #include #include #include #include "math_private.h" /* * gcc doesn't implement complex multiplication or division correctly, * so we need to handle infinities specially. We turn on this pragma to * notify conforming c99 compilers that the fast-but-incorrect code that * gcc generates is acceptable, since the special cases have already been * handled. */ #pragma STDC CX_LIMITED_RANGE ON /* * We risk spurious overflow for components >= LDBL_MAX / (1 + sqrt(2)). * Rather than determining the fully precise value at which we might * overflow, just use a threshold of approximately LDBL_MAX / 4. */ #if LDBL_MAX_EXP != 0x4000 #error "Unsupported long double format" #else #define THRESH 0x1p16382L #endif long double complex csqrtl(long double complex z) { long double complex result; long double a, b; long double t; int scale; a = creall(z); b = cimagl(z); /* Handle special cases. */ if (z == 0) return (CMPLXL(0, b)); if (isinf(b)) return (CMPLXL(INFINITY, b)); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ - return (CMPLXL(a, t)); /* return NaN + NaN i */ + return (CMPLXL(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */ } if (isinf(a)) { /* * csqrt(inf + NaN i) = inf + NaN i * csqrt(inf + y i) = inf + 0 i * csqrt(-inf + NaN i) = NaN +- inf i * csqrt(-inf + y i) = 0 + inf i */ if (signbit(a)) return (CMPLXL(fabsl(b - b), copysignl(a, b))); else return (CMPLXL(a, copysignl(b - b, b))); } - /* - * The remaining special case (b is NaN) is handled just fine by - * the normal code path below. - */ + if (isnan(b)) { + t = (a - a) / (a - a); /* raise invalid */ + return (CMPLXL(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */ + } /* Scale to avoid overflow. */ if (fabsl(a) >= THRESH || fabsl(b) >= THRESH) { a *= 0.25; b *= 0.25; scale = 1; } else { scale = 0; } /* Algorithm 312, CACM vol 10, Oct 1967. */ if (a >= 0) { t = sqrtl((a + hypotl(a, b)) * 0.5); result = CMPLXL(t, b / (2 * t)); } else { t = sqrtl((-a + hypotl(a, b)) * 0.5); result = CMPLXL(fabsl(b) / (2 * t), copysignl(t, b)); } /* Rescale. */ if (scale) return (result * 2); else return (result); } Index: head/lib/msun/src/s_ctanh.c =================================================================== --- head/lib/msun/src/s_ctanh.c (revision 336361) +++ head/lib/msun/src/s_ctanh.c (revision 336362) @@ -1,147 +1,147 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2011 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Hyperbolic tangent of a complex argument z = x + I y. * * The algorithm is from: * * W. Kahan. Branch Cuts for Complex Elementary Functions or Much * Ado About Nothing's Sign Bit. In The State of the Art in * Numerical Analysis, pp. 165 ff. Iserles and Powell, eds., 1987. * * Method: * * Let t = tan(x) * beta = 1/cos^2(y) * s = sinh(x) * rho = cosh(x) * * We have: * * tanh(z) = sinh(z) / cosh(z) * * sinh(x) cos(y) + I cosh(x) sin(y) * = --------------------------------- * cosh(x) cos(y) + I sinh(x) sin(y) * * cosh(x) sinh(x) / cos^2(y) + I tan(y) * = ------------------------------------- * 1 + sinh^2(x) / cos^2(y) * * beta rho s + I t * = ---------------- * 1 + beta s^2 * * Modifications: * * I omitted the original algorithm's handling of overflow in tan(x) after * verifying with nearpi.c that this can't happen in IEEE single or double * precision. I also handle large x differently. */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" double complex ctanh(double complex z) { double x, y; double t, beta, s, rho, denom; uint32_t hx, ix, lx; x = creal(z); y = cimag(z); EXTRACT_WORDS(hx, lx, x); ix = hx & 0x7fffffff; /* * ctanh(NaN +- I 0) = d(NaN) +- I 0 * * ctanh(NaN + I y) = d(NaN,y) + I d(NaN,y) for y != 0 * * The imaginary part has the sign of x*sin(2*y), but there's no * special effort to get this right. * * ctanh(+-Inf +- I Inf) = +-1 +- I 0 * * ctanh(+-Inf + I y) = +-1 + I 0 sin(2y) for y finite * * The imaginary part of the sign is unspecified. This special * case is only needed to avoid a spurious invalid exception when * y is infinite. */ if (ix >= 0x7ff00000) { if ((ix & 0xfffff) | lx) /* x is NaN */ - return (CMPLX((x + 0) * (y + 0), - y == 0 ? y : (x + 0) * (y + 0))); + return (CMPLX(nan_mix(x, y), + y == 0 ? y : nan_mix(x, y))); SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */ return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); } /* * ctanh(x + I NaN) = d(NaN) + I d(NaN) * ctanh(x +- I Inf) = dNaN + I dNaN */ if (!isfinite(y)) return (CMPLX(y - y, y - y)); /* * ctanh(+-huge +- I y) ~= +-1 +- I 2sin(2y)/exp(2x), using the * approximation sinh^2(huge) ~= exp(2*huge) / 4. * We use a modified formula to avoid spurious overflow. */ if (ix >= 0x40360000) { /* |x| >= 22 */ double exp_mx = exp(-fabs(x)); return (CMPLX(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx)); } /* Kahan's algorithm */ t = tan(y); beta = 1.0 + t * t; /* = 1 / cos^2(y) */ s = sinh(x); rho = sqrt(1 + s * s); /* = cosh(x) */ denom = 1 + beta * s * s; return (CMPLX((beta * rho * s) / denom, t / denom)); } double complex ctan(double complex z) { /* ctan(z) = -I * ctanh(I * z) = I * conj(ctanh(I * conj(z))) */ z = ctanh(CMPLX(cimag(z), creal(z))); return (CMPLX(cimag(z), creal(z))); } Index: head/lib/msun/src/s_ctanhf.c =================================================================== --- head/lib/msun/src/s_ctanhf.c (revision 336361) +++ head/lib/msun/src/s_ctanhf.c (revision 336362) @@ -1,87 +1,87 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2011 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * Hyperbolic tangent of a complex argument z. See s_ctanh.c for details. */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" float complex ctanhf(float complex z) { float x, y; float t, beta, s, rho, denom; uint32_t hx, ix; x = crealf(z); y = cimagf(z); GET_FLOAT_WORD(hx, x); ix = hx & 0x7fffffff; if (ix >= 0x7f800000) { if (ix & 0x7fffff) - return (CMPLXF((x + 0) * (y + 0), - y == 0 ? y : (x + 0) * (y + 0))); + return (CMPLXF(nan_mix(x, y), + y == 0 ? y : nan_mix(x, y))); SET_FLOAT_WORD(x, hx - 0x40000000); return (CMPLXF(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y)))); } if (!isfinite(y)) return (CMPLXF(y - y, y - y)); if (ix >= 0x41300000) { /* |x| >= 11 */ float exp_mx = expf(-fabsf(x)); return (CMPLXF(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx)); } t = tanf(y); beta = 1.0 + t * t; s = sinhf(x); rho = sqrtf(1 + s * s); denom = 1 + beta * s * s; return (CMPLXF((beta * rho * s) / denom, t / denom)); } float complex ctanf(float complex z) { z = ctanhf(CMPLXF(cimagf(z), crealf(z))); return (CMPLXF(cimagf(z), crealf(z))); } Index: head/lib/msun/src/s_remquo.c =================================================================== --- head/lib/msun/src/s_remquo.c (revision 336361) +++ head/lib/msun/src/s_remquo.c (revision 336362) @@ -1,159 +1,159 @@ /* @(#)e_fmod.c 1.3 95/01/18 */ /*- * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); #include #include "math.h" #include "math_private.h" static const double Zero[] = {0.0, -0.0,}; /* * Return the IEEE remainder and set *quo to the last n bits of the * quotient, rounded to the nearest integer. We choose n=31 because * we wind up computing all the integer bits of the quotient anyway as * a side-effect of computing the remainder by the shift and subtract * method. In practice, this is far more bits than are needed to use * remquo in reduction algorithms. */ double remquo(double x, double y, int *quo) { int32_t n,hx,hy,hz,ix,iy,sx,i; u_int32_t lx,ly,lz,q,sxy; EXTRACT_WORDS(hx,lx,x); EXTRACT_WORDS(hy,ly,y); sxy = (hx ^ hy) & 0x80000000; sx = hx&0x80000000; /* sign of x */ hx ^=sx; /* |x| */ hy &= 0x7fffffff; /* |y| */ /* purge off exception values */ if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */ - return (x*y)/(x*y); + return nan_mix(x, y)/nan_mix(x, y); if(hx<=hy) { if((hx>31]; /* |x|=|y| return x*0*/ } } /* determine ix = ilogb(x) */ if(hx<0x00100000) { /* subnormal x */ if(hx==0) { for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; } else { for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; } } else ix = (hx>>20)-1023; /* determine iy = ilogb(y) */ if(hy<0x00100000) { /* subnormal y */ if(hy==0) { for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; } else { for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; } } else iy = (hy>>20)-1023; /* set up {hx,lx}, {hy,ly} and align y to x */ if(ix >= -1022) hx = 0x00100000|(0x000fffff&hx); else { /* subnormal x, shift x to normal */ n = -1022-ix; if(n<=31) { hx = (hx<>(32-n)); lx <<= n; } else { hx = lx<<(n-32); lx = 0; } } if(iy >= -1022) hy = 0x00100000|(0x000fffff&hy); else { /* subnormal y, shift y to normal */ n = -1022-iy; if(n<=31) { hy = (hy<>(32-n)); ly <<= n; } else { hy = ly<<(n-32); ly = 0; } } /* fix point fmod */ n = ix - iy; q = 0; while(n--) { hz=hx-hy;lz=lx-ly; if(lx>31); lx = lx+lx;} else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;} q <<= 1; } hz=hx-hy;lz=lx-ly; if(lx=0) {hx=hz;lx=lz;q++;} /* convert back to floating value and restore the sign */ if((hx|lx)==0) { /* return sign(x)*0 */ q &= 0x7fffffff; *quo = (sxy ? -q : q); return Zero[(u_int32_t)sx>>31]; } while(hx<0x00100000) { /* normalize x */ hx = hx+hx+(lx>>31); lx = lx+lx; iy -= 1; } if(iy>= -1022) { /* normalize output */ hx = ((hx-0x00100000)|((iy+1023)<<20)); } else { /* subnormal output */ n = -1022 - iy; if(n<=20) { lx = (lx>>n)|((u_int32_t)hx<<(32-n)); hx >>= n; } else if (n<=31) { lx = (hx<<(32-n))|(lx>>n); hx = 0; } else { lx = hx>>(n-32); hx = 0; } } fixup: INSERT_WORDS(x,hx,lx); y = fabs(y); if (y < 0x1p-1021) { if (x+x>y || (x+x==y && (q & 1))) { q++; x-=y; } } else if (x>0.5*y || (x==0.5*y && (q & 1))) { q++; x-=y; } GET_HIGH_WORD(hx,x); SET_HIGH_WORD(x,hx^sx); q &= 0x7fffffff; *quo = (sxy ? -q : q); return x; } #if LDBL_MANT_DIG == 53 __weak_reference(remquo, remquol); #endif Index: head/lib/msun/src/s_remquof.c =================================================================== --- head/lib/msun/src/s_remquof.c (revision 336361) +++ head/lib/msun/src/s_remquof.c (revision 336362) @@ -1,122 +1,122 @@ /* @(#)e_fmod.c 1.3 95/01/18 */ /*- * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" static const float Zero[] = {0.0, -0.0,}; /* * Return the IEEE remainder and set *quo to the last n bits of the * quotient, rounded to the nearest integer. We choose n=31 because * we wind up computing all the integer bits of the quotient anyway as * a side-effect of computing the remainder by the shift and subtract * method. In practice, this is far more bits than are needed to use * remquo in reduction algorithms. */ float remquof(float x, float y, int *quo) { int32_t n,hx,hy,hz,ix,iy,sx,i; u_int32_t q,sxy; GET_FLOAT_WORD(hx,x); GET_FLOAT_WORD(hy,y); sxy = (hx ^ hy) & 0x80000000; sx = hx&0x80000000; /* sign of x */ hx ^=sx; /* |x| */ hy &= 0x7fffffff; /* |y| */ /* purge off exception values */ if(hy==0||hx>=0x7f800000||hy>0x7f800000) /* y=0,NaN;or x not finite */ - return (x*y)/(x*y); + return nan_mix(x, y)/nan_mix(x, y); if(hx>31]; /* |x|=|y| return x*0*/ } /* determine ix = ilogb(x) */ if(hx<0x00800000) { /* subnormal x */ for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1; } else ix = (hx>>23)-127; /* determine iy = ilogb(y) */ if(hy<0x00800000) { /* subnormal y */ for (iy = -126,i=(hy<<8); i>0; i<<=1) iy -=1; } else iy = (hy>>23)-127; /* set up {hx,lx}, {hy,ly} and align y to x */ if(ix >= -126) hx = 0x00800000|(0x007fffff&hx); else { /* subnormal x, shift x to normal */ n = -126-ix; hx <<= n; } if(iy >= -126) hy = 0x00800000|(0x007fffff&hy); else { /* subnormal y, shift y to normal */ n = -126-iy; hy <<= n; } /* fix point fmod */ n = ix - iy; q = 0; while(n--) { hz=hx-hy; if(hz<0) hx = hx << 1; else {hx = hz << 1; q++;} q <<= 1; } hz=hx-hy; if(hz>=0) {hx=hz;q++;} /* convert back to floating value and restore the sign */ if(hx==0) { /* return sign(x)*0 */ q &= 0x7fffffff; *quo = (sxy ? -q : q); return Zero[(u_int32_t)sx>>31]; } while(hx<0x00800000) { /* normalize x */ hx <<= 1; iy -= 1; } if(iy>= -126) { /* normalize output */ hx = ((hx-0x00800000)|((iy+127)<<23)); } else { /* subnormal output */ n = -126 - iy; hx >>= n; } fixup: SET_FLOAT_WORD(x,hx); y = fabsf(y); if (y < 0x1p-125f) { if (x+x>y || (x+x==y && (q & 1))) { q++; x-=y; } } else if (x>0.5f*y || (x==0.5f*y && (q & 1))) { q++; x-=y; } GET_FLOAT_WORD(hx,x); SET_FLOAT_WORD(x,hx^sx); q &= 0x7fffffff; *quo = (sxy ? -q : q); return x; } Index: head/lib/msun/src/s_remquol.c =================================================================== --- head/lib/msun/src/s_remquol.c (revision 336361) +++ head/lib/msun/src/s_remquol.c (revision 336362) @@ -1,178 +1,178 @@ /* @(#)e_fmod.c 1.3 95/01/18 */ /*- * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); #include #include #include "fpmath.h" #include "math.h" #include "math_private.h" #define BIAS (LDBL_MAX_EXP - 1) #if LDBL_MANL_SIZE > 32 typedef uint64_t manl_t; #else typedef uint32_t manl_t; #endif #if LDBL_MANH_SIZE > 32 typedef uint64_t manh_t; #else typedef uint32_t manh_t; #endif /* * These macros add and remove an explicit integer bit in front of the * fractional mantissa, if the architecture doesn't have such a bit by * default already. */ #ifdef LDBL_IMPLICIT_NBIT #define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE)) #define HFRAC_BITS LDBL_MANH_SIZE #else #define SET_NBIT(hx) (hx) #define HFRAC_BITS (LDBL_MANH_SIZE - 1) #endif #define MANL_SHIFT (LDBL_MANL_SIZE - 1) static const long double Zero[] = {0.0L, -0.0L}; /* * Return the IEEE remainder and set *quo to the last n bits of the * quotient, rounded to the nearest integer. We choose n=31 because * we wind up computing all the integer bits of the quotient anyway as * a side-effect of computing the remainder by the shift and subtract * method. In practice, this is far more bits than are needed to use * remquo in reduction algorithms. * * Assumptions: * - The low part of the mantissa fits in a manl_t exactly. * - The high part of the mantissa fits in an int64_t with enough room * for an explicit integer bit in front of the fractional bits. */ long double remquol(long double x, long double y, int *quo) { union IEEEl2bits ux, uy; int64_t hx,hz; /* We need a carry bit even if LDBL_MANH_SIZE is 32. */ manh_t hy; manl_t lx,ly,lz; int ix,iy,n,q,sx,sxy; ux.e = x; uy.e = y; sx = ux.bits.sign; sxy = sx ^ uy.bits.sign; ux.bits.sign = 0; /* |x| */ uy.bits.sign = 0; /* |y| */ x = ux.e; /* purge off exception values */ if((uy.bits.exp|uy.bits.manh|uy.bits.manl)==0 || /* y=0 */ (ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */ (uy.bits.exp == BIAS + LDBL_MAX_EXP && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */ - return (x*y)/(x*y); + return nan_mix(x, y)/nan_mix(x, y); if(ux.bits.exp<=uy.bits.exp) { if((ux.bits.exp>MANL_SHIFT); lx = lx+lx;} else {hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz; q++;} q <<= 1; } hz=hx-hy;lz=lx-ly; if(lx=0) {hx=hz;lx=lz;q++;} /* convert back to floating value and restore the sign */ if((hx|lx)==0) { /* return sign(x)*0 */ q &= 0x7fffffff; *quo = (sxy ? -q : q); return Zero[sx]; } while(hx<(1ULL<>MANL_SHIFT); lx = lx+lx; iy -= 1; } ux.bits.manh = hx; /* The integer bit is truncated here if needed. */ ux.bits.manl = lx; if (iy < LDBL_MIN_EXP) { ux.bits.exp = iy + (BIAS + 512); ux.e *= 0x1p-512; } else { ux.bits.exp = iy + BIAS; } ux.bits.sign = 0; x = ux.e; fixup: y = fabsl(y); if (y < LDBL_MIN * 2) { if (x+x>y || (x+x==y && (q & 1))) { q++; x-=y; } } else if (x>0.5*y || (x==0.5*y && (q & 1))) { q++; x-=y; } ux.e = x; ux.bits.sign ^= sx; x = ux.e; q &= 0x7fffffff; *quo = (sxy ? -q : q); return x; }