Date: Fri, 27 Jul 2018 17:39:36 +0000 (UTC) From: Dimitry Andric <dim@FreeBSD.org> To: src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-stable@freebsd.org, svn-src-stable-11@freebsd.org Subject: svn commit: r336767 - in stable/11: include lib/msun lib/msun/ld128 lib/msun/ld80 lib/msun/man lib/msun/src Message-ID: <201807271739.w6RHdaIZ087956@repo.freebsd.org>
next in thread | raw e-mail | index | archive | help
Author: dim Date: Fri Jul 27 17:39:36 2018 New Revision: 336767 URL: https://svnweb.freebsd.org/changeset/base/336767 Log: MFC r327400 (by eadler): cacos(3): correct spelling of 'I' In some cases we had 'i' instead of 'I'. PR: 195517 Submitted by: stephen MFC r329259 (by eadler): msun: signed overflow in atan2 As a component of atan2(y, x), the case of x == 1.0 is farmed out to atan(y). The current implementation of this comparison is vulnerable to signed integer underflow (that is, undefined behavior), and it's performed in a somewhat more complicated way than it need be. Change it to not be quite so cute, rather directly comparing the high/low bits of x to the specific IEEE-754 bit pattern that encodes 1.0. Note that while there are three different e_atan* files in the relevant directory, only this one needs fixing. e_atan2f.c already compares against the full bit pattern encoding 1.0f, while e_atan2l.cuses bitwise-ands/ors/nots and so doesn't require a change. Closes #130 Submitted by: Jeff Walden (@jswalden github PR #130) Reviewed by: bde MFC r334721 (by cem): clog.3, complex.3: Fix typos and igor style issues PR: 228783 Reported by: Karsten <freebsd-bugzilla AT kkoenig.net> MFC r336299 (by mmacy): msun: add ld80/ld128 powl, cpow, cpowf, cpowl from openbsd This corresponds to the latest status (hasn't changed in 9+ years) from openbsd of ld80/ld128 powl, and source cpowf, cpow, cpowl (the complex power functions for float complex, double complex, and long double complex) which are required for C99 compliance and were missing from FreeBSD. Also required for some numerical codes using complex numbered Hamiltonians. Thanks to jhb for tracking down the issue with making weak_reference compile on powerpc. When asked to review, bde said "I don't like it" - but provided no actionable feedback or superior implementations. Discussed with: jhb Submitted by: jmd Differential Revision: https://reviews.freebsd.org/D15919 MFC r336563: Recommit r336497: Fix powl, cpow, cpowf, and cpowl imports from OpenBSD This is a follow-up to r336299. * lib/msun/Makefile: . Remove polevll.c * lib/msun/ld80/e_powl.c: . Copy contents of polevll.c to here. This is the only consumer of these functions. Make functions 'static inline'. . Make reducl a 'static inline' function. * lib/msun/man/exp.3: . Remove BUGS section that no longer applies. * lib/msun/src/math_private.h: . Remove prototypes of __p1evll() and __polevll() * lib/msun/src/s_cpow.c: * lib/msun/src/s_cpowf.c: * lib/msun/src/s_cpowl.c . Include math_private.h. . Use the CMPLX macro from either C99 or math_private.h (depends on compiler support) instead of the problematic use of complex I. Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu> PR: 229876 Added: stable/11/lib/msun/ld128/e_powl.c - copied unchanged from r336299, head/lib/msun/ld128/e_powl.c stable/11/lib/msun/ld80/e_powl.c - copied, changed from r336299, head/lib/msun/ld80/e_powl.c stable/11/lib/msun/man/cpow.3 - copied unchanged from r336299, head/lib/msun/man/cpow.3 stable/11/lib/msun/src/s_cpow.c - copied, changed from r336299, head/lib/msun/src/s_cpow.c stable/11/lib/msun/src/s_cpowf.c - copied, changed from r336299, head/lib/msun/src/s_cpowf.c stable/11/lib/msun/src/s_cpowl.c - copied, changed from r336299, head/lib/msun/src/s_cpowl.c Modified: stable/11/include/complex.h stable/11/lib/msun/Makefile stable/11/lib/msun/Symbol.map stable/11/lib/msun/man/cacos.3 stable/11/lib/msun/man/clog.3 stable/11/lib/msun/man/complex.3 stable/11/lib/msun/man/exp.3 stable/11/lib/msun/src/e_atan2.c stable/11/lib/msun/src/e_pow.c stable/11/lib/msun/src/imprecise.c stable/11/lib/msun/src/math_private.h Directory Properties: stable/11/ (props changed) Modified: stable/11/include/complex.h ============================================================================== --- stable/11/include/complex.h Fri Jul 27 16:13:06 2018 (r336766) +++ stable/11/include/complex.h Fri Jul 27 17:39:36 2018 (r336767) @@ -95,6 +95,10 @@ double complex conj(double complex) __pure2; float complex conjf(float complex) __pure2; long double complex conjl(long double complex) __pure2; +float complex cpowf(float complex, float complex); +double complex cpow(double complex, double complex); +long double complex + cpowl(long double complex, long double complex); float complex cprojf(float complex) __pure2; double complex cproj(double complex) __pure2; long double complex Modified: stable/11/lib/msun/Makefile ============================================================================== --- stable/11/lib/msun/Makefile Fri Jul 27 16:13:06 2018 (r336766) +++ stable/11/lib/msun/Makefile Fri Jul 27 17:39:36 2018 (r336767) @@ -98,7 +98,7 @@ COMMON_SRCS+= s_copysignl.c s_fabsl.c s_llrintl.c s_lr COMMON_SRCS+= catrigl.c \ e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \ e_coshl.c e_fmodl.c e_hypotl.c \ - e_lgammal.c e_lgammal_r.c \ + e_lgammal.c e_lgammal_r.c e_powl.c \ e_remainderl.c e_sinhl.c e_sqrtl.c \ invtrig.c k_cosl.c k_sinl.c k_tanl.c \ s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c \ @@ -115,6 +115,7 @@ COMMON_SRCS+= catrig.c catrigf.c \ s_ccosh.c s_ccoshf.c s_cexp.c s_cexpf.c \ s_cimag.c s_cimagf.c s_cimagl.c \ s_conj.c s_conjf.c s_conjl.c \ + s_cpow.c s_cpowf.c s_cpowl.c \ s_cproj.c s_cprojf.c s_creal.c s_crealf.c s_creall.c \ s_csinh.c s_csinhf.c s_ctanh.c s_ctanhf.c @@ -134,7 +135,7 @@ INCS+= fenv.h math.h MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \ ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \ - cimag.3 clog.3 copysign.3 cos.3 cosh.3 csqrt.3 erf.3 \ + cimag.3 clog.3 copysign.3 cos.3 cosh.3 cpow.3 csqrt.3 erf.3 \ exp.3 fabs.3 fdim.3 \ feclearexcept.3 feenableexcept.3 fegetenv.3 \ fegetround.3 fenv.3 floor.3 \ @@ -172,6 +173,7 @@ MLINKS+=clog.3 clogf.3 clog.3 clogl.3 MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3 MLINKS+=cos.3 cosf.3 cos.3 cosl.3 MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3 +MLINKS+=cpow.3 cpowf.3 cpow.3 cpowl.3 MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3 MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 erf.3 erfl.3 erf.3 erfcl.3 MLINKS+=exp.3 expm1.3 exp.3 expm1f.3 exp.3 expm1l.3 exp.3 pow.3 exp.3 powf.3 \ Modified: stable/11/lib/msun/Symbol.map ============================================================================== --- stable/11/lib/msun/Symbol.map Fri Jul 27 16:13:06 2018 (r336766) +++ stable/11/lib/msun/Symbol.map Fri Jul 27 17:39:36 2018 (r336767) @@ -274,10 +274,10 @@ FBSD_1.3 { log1pl; log2l; logl; + powl; sinhl; tanhl; /* Implemented as weak aliases for imprecise versions */ - powl; tgammal; }; @@ -297,6 +297,9 @@ FBSD_1.5 { clog; clogf; clogl; + cpow; + cpowf; + cpowl; sincos; sincosf; sincosl; Copied: stable/11/lib/msun/ld128/e_powl.c (from r336299, head/lib/msun/ld128/e_powl.c) ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ stable/11/lib/msun/ld128/e_powl.c Fri Jul 27 17:39:36 2018 (r336767, copy of r336299, head/lib/msun/ld128/e_powl.c) @@ -0,0 +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 <steve@moshier.net> + * + * 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 <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#include <float.h> +#include <math.h> + +#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; + + /* 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|<sqrt(3/2) */ + else if (j < 0xbb67) + k = 1; /* |x|<sqrt(3) */ + else + { + k = 0; + n += 1; + ix -= 0x00010000; + } + + o.value = ax; + o.parts32.mswhi = ix; + ax = o.value; + + /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ + u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */ + v = one / (ax + bp[k]); + s = u * v; + s_h = s; + + o.value = s_h; + o.parts32.lswlo = 0; + o.parts32.lswhi &= 0xf8000000; + s_h = o.value; + /* t_h=ax+bp[k] High */ + t_h = ax + bp[k]; + o.value = t_h; + o.parts32.lswlo = 0; + o.parts32.lswhi &= 0xf8000000; + t_h = o.value; + 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; + u = LN[0] + s2 * (LN[1] + s2 * (LN[2] + s2 * (LN[3] + s2 * LN[4]))); + v = LD[0] + s2 * (LD[1] + s2 * (LD[2] + s2 * (LD[3] + s2 * (LD[4] + s2)))); + r = s2 * s2 * u / v; + r += s_l * (s_h + s); + s2 = s_h * s_h; + t_h = 3.0 + s2 + r; + o.value = t_h; + o.parts32.lswlo = 0; + o.parts32.lswhi &= 0xf8000000; + t_h = o.value; + t_l = r - ((t_h - 3.0) - 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; + o.value = p_h; + o.parts32.lswlo = 0; + o.parts32.lswhi &= 0xf8000000; + p_h = o.value; + 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 = (long double) n; + t1 = (((z_h + z_l) + dp_h[k]) + t); + o.value = t1; + o.parts32.lswlo = 0; + o.parts32.lswhi &= 0xf8000000; + t1 = o.value; + t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); + + /* s (sign of result -ve**odd) = -1 else = 1 */ + s = one; + if (((((u_int32_t) hx >> 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; +} Copied and modified: stable/11/lib/msun/ld80/e_powl.c (from r336299, head/lib/msun/ld80/e_powl.c) ============================================================================== --- head/lib/msun/ld80/e_powl.c Sun Jul 15 00:23:10 2018 (r336299, copy source) +++ stable/11/lib/msun/ld80/e_powl.c Fri Jul 27 17:39:36 2018 (r336767) @@ -14,6 +14,52 @@ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#include <math.h> + +#include "math_private.h" + +/* + * Polynomial evaluator: + * P[0] x^n + P[1] x^(n-1) + ... + P[n] + */ +static inline long double +__polevll(long double x, long double *PP, int n) +{ + long double y; + long double *P; + + P = PP; + y = *P++; + do { + y = y * x + *P++; + } while (--n); + + return (y); +} + +/* + * Polynomial evaluator: + * x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n] + */ +static inline long double +__p1evll(long double x, long double *PP, int n) +{ + long double y; + long double *P; + + P = PP; + n -= 1; + y = x + *P++; + do { + y = y * x + *P++; + } while (--n); + + return (y); +} + /* powl.c * * Power function, long double precision @@ -467,7 +513,7 @@ return( z ); /* Find a multiple of 1/NXT that is within 1/NXT of x. */ -static long double +static inline long double reducl(long double x) { long double t; Modified: stable/11/lib/msun/man/cacos.3 ============================================================================== --- stable/11/lib/msun/man/cacos.3 Fri Jul 27 16:13:06 2018 (r336766) +++ stable/11/lib/msun/man/cacos.3 Fri Jul 27 17:39:36 2018 (r336767) @@ -136,9 +136,9 @@ corresponding ranges for the return values, adopted by .It Sy Function Ta Sy Branch Cut(s) Ta Sy Range .It cacos Ta (-\*(If, -1) \*(Un (1, \*(If) Ta [0, \*(Pi] .It casin Ta (-\*(If, -1) \*(Un (1, \*(If) Ta [-\*(Pi/2, \*(Pi/2] -.It catan Ta (-\*(If*I, -i) \*(Un (I, \*(If*I) Ta [-\*(Pi/2, \*(Pi/2] +.It catan Ta (-\*(If*I, -I) \*(Un (I, \*(If*I) Ta [-\*(Pi/2, \*(Pi/2] .It cacosh Ta (-\*(If, 1) Ta [-\*(Pi*I, \*(Pi*I] -.It casinh Ta (-\*(If*I, -i) \*(Un (I, \*(If*I) Ta [-\*(Pi/2*I, \*(Pi/2*I] +.It casinh Ta (-\*(If*I, -I) \*(Un (I, \*(If*I) Ta [-\*(Pi/2*I, \*(Pi/2*I] .It catanh Ta (-\*(If, -1) \*(Un (1, \*(If) Ta [-\*(Pi/2*I, \*(Pi/2*I] .El .Sh SEE ALSO Modified: stable/11/lib/msun/man/clog.3 ============================================================================== --- stable/11/lib/msun/man/clog.3 Fri Jul 27 16:13:06 2018 (r336766) +++ stable/11/lib/msun/man/clog.3 Fri Jul 27 17:39:36 2018 (r336767) @@ -24,7 +24,7 @@ .\" .\" $FreeBSD$ .\" -.Dd February 13, 2017 +.Dd June 6, 2018 .Dt CLOG 3 .Os .Sh NAME @@ -32,7 +32,7 @@ .Nm clogf , and .Nm clogl -.Nd complex natural logrithm functions +.Nd complex natural logarithm functions .Sh LIBRARY .Lb libm .Sh SYNOPSIS @@ -47,9 +47,9 @@ and The .Fn clog , .Fn clogf , -and +and .Fn clogl -functions compute the complex natural logrithm of +functions compute the complex natural logarithm of .Fa z . with a branch cut along the negative real axis . .Sh RETURN VALUES @@ -58,14 +58,13 @@ The function returns the complex natural logarithm value, in the range of a strip mathematically unbounded along the real axis and in the interval [-I* \*(Pi , +I* \*(Pi ] along the imaginary axis. -The function satisfies the relationship: +The function satisfies the relationship: .Fo clog .Fn conj "z" Fc = .Fo conj .Fn clog "z" Fc . .Pp - .\" Table is formatted for an 80-column xterm. .Bl -column ".Sy +\*(If + I*\*(Na" ".Sy Return value" ".Sy Divide-by-zero exception" .It Sy Argument Ta Sy Return value Ta Sy Comment @@ -86,9 +85,8 @@ The function satisfies the relationship: .It Ta Ta floating-point exception .It Ta Ta for finite y .It \*(Na + I*\*(If Ta +\*(If + I*\*(Na -.It \*(Na + I*\*(Na Ta \*(Na + I*\*(Na +.It \*(Na + I*\*(Na Ta \*(Na + I*\*(Na .El - .Sh SEE ALSO .Xr complex 3 , .Xr log 3 , Modified: stable/11/lib/msun/man/complex.3 ============================================================================== --- stable/11/lib/msun/man/complex.3 Fri Jul 27 16:13:06 2018 (r336766) +++ stable/11/lib/msun/man/complex.3 Fri Jul 27 17:39:36 2018 (r336767) @@ -24,7 +24,7 @@ .\" .\" $FreeBSD$ .\" -.Dd May 13, 2018 +.Dd June 19, 2018 .Dt COMPLEX 3 .Os .Sh NAME @@ -63,28 +63,28 @@ and .Fn cabsl "long double complex z" , respectively. .de Cl -.Bl -column "csqrt" "complex absolute value (i.e. norm, modulus, magnitude)" +.Bl -column "csqrt" "complex absolute value (i.e., norm, modulus, magnitude)" .Em "Name Description" .. .\" Section 7.3.5 - 7.3.7 of ISO C99 standard unimplemented, see BUGS .\" Section 7.3.8 of ISO C99 standard .Ss Absolute-value Functions .Cl -cabs complex absolute value (i.e. norm, modulus, magnitude) +cabs complex absolute value (i.e., norm, modulus, magnitude) csqrt complex square root .El .Ss Exponential Function .Cl cexp exponential base e .El -.Ss Natural logrithm Function +.Ss Natural logarithm Function .Cl -clog natural logrithm +clog natural logarithm .El .\" Section 7.3.9 of ISO C99 standard .Ss Manipulation Functions .Cl -carg compute the argument (i.e. phase angle) +carg compute the argument (i.e., phase angle) cimag compute the imaginary part conj compute the complex conjugate cproj compute projection onto Riemann sphere @@ -101,6 +101,7 @@ catan arc tangent catanh arc hyperbolic tangent ccos cosine ccosh hyperbolic cosine +cpow power function csin sine csinh hyperbolic sine ctan tangent @@ -120,7 +121,3 @@ The .In complex.h functions described here conform to .St -isoC-99 . -.Sh BUGS -The power functions -.Fn cpow -are not implemented. Copied: stable/11/lib/msun/man/cpow.3 (from r336299, head/lib/msun/man/cpow.3) ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ stable/11/lib/msun/man/cpow.3 Fri Jul 27 17:39:36 2018 (r336767, copy of r336299, head/lib/msun/man/cpow.3) @@ -0,0 +1,63 @@ +.\" Copyright (c) 2011 Martynas Venckus <martynas@openbsd.org> +.\" +.\" 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. +.\" $FreeBSD$ +.\" +.Dd $Mdocdate: June 5 2013 $ +.Dt CPOW 3 +.Os +.Sh NAME +.Nm cpow , +.Nm cpowf , +.Nm cpowl +.Nd complex power functions +.Sh SYNOPSIS +.In complex.h +.Ft double complex +.Fn cpow "double complex x" "double complex z" +.Ft float complex +.Fn cpowf "float complex x" "float complex z" +.Ft long double complex +.Fn cpowl "long double complex x" "long double complex z" +.Sh DESCRIPTION +The +.Fn cpow , +.Fn cpowf +and +.Fn cpowl +functions compute the complex number +.Fa x +raised to the complex power +.Fa z , +with a branch cut along the negative real axis for the first argument. +.Sh RETURN VALUES +The +.Fn cpow , +.Fn cpowf +and +.Fn cpowl +functions return the complex number +.Fa x +raised to the complex power +.Fa z . +.Sh SEE ALSO +.Xr cexp 3 , +.Xr clog 3 +.Sh STANDARDS +The +.Fn cpow , +.Fn cpowf +and +.Fn cpowl +functions conform to +.St -isoC-99 . Modified: stable/11/lib/msun/man/exp.3 ============================================================================== --- stable/11/lib/msun/man/exp.3 Fri Jul 27 16:13:06 2018 (r336766) +++ stable/11/lib/msun/man/exp.3 Fri Jul 27 17:39:36 2018 (r336767) @@ -180,16 +180,9 @@ If 0**0 = 1, then then \*(Na**0 = 1 too because x**0 = 1 for all finite and infinite x, i.e., independently of x. .El -.Sh BUGS -To conform with newer C/C++ standards, a stub implementation for -.Nm powl -was committed to the math library, where -.Nm powl -is mapped to -.Nm pow . -Thus, the numerical accuracy is at most that of the 53-bit double -precision implementation. .Sh SEE ALSO +.Xr clog 3 +.Xr cpow 3 .Xr fenv 3 , .Xr ldexp 3 , .Xr log 3 , Modified: stable/11/lib/msun/src/e_atan2.c ============================================================================== --- stable/11/lib/msun/src/e_atan2.c Fri Jul 27 16:13:06 2018 (r336766) +++ stable/11/lib/msun/src/e_atan2.c Fri Jul 27 17:39:36 2018 (r336767) @@ -71,7 +71,7 @@ __ieee754_atan2(double y, double x) if(((ix|((lx|-lx)>>31))>0x7ff00000)|| ((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */ return x+y; - if((hx-0x3ff00000|lx)==0) return atan(y); /* x=1.0 */ + 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 */ Modified: stable/11/lib/msun/src/e_pow.c ============================================================================== --- stable/11/lib/msun/src/e_pow.c Fri Jul 27 16:13:06 2018 (r336766) +++ stable/11/lib/msun/src/e_pow.c Fri Jul 27 17:39:36 2018 (r336767) @@ -57,6 +57,7 @@ __FBSDID("$FreeBSD$"); * to produce the hexadecimal values shown. */ +#include <float.h> #include "math.h" #include "math_private.h" @@ -307,3 +308,7 @@ __ieee754_pow(double x, double y) else SET_HIGH_WORD(z,j); return s*z; } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(pow, powl); +#endif Modified: stable/11/lib/msun/src/imprecise.c ============================================================================== --- stable/11/lib/msun/src/imprecise.c Fri Jul 27 16:13:06 2018 (r336766) +++ stable/11/lib/msun/src/imprecise.c Fri Jul 27 17:39:36 2018 (r336767) @@ -48,14 +48,6 @@ __weak_reference(imprecise_## x, x);\ WARN_IMPRECISE(x) -long double -imprecise_powl(long double x, long double y) -{ - - return pow(x, y); -} -DECLARE_WEAK(powl); - #define DECLARE_IMPRECISE(f) \ long double imprecise_ ## f ## l(long double v) { return f(v); }\ DECLARE_WEAK(f ## l) Modified: stable/11/lib/msun/src/math_private.h ============================================================================== --- stable/11/lib/msun/src/math_private.h Fri Jul 27 16:13:06 2018 (r336766) +++ stable/11/lib/msun/src/math_private.h Fri Jul 27 17:39:36 2018 (r336767) @@ -48,6 +48,47 @@ #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 Copied and modified: stable/11/lib/msun/src/s_cpow.c (from r336299, head/lib/msun/src/s_cpow.c) ============================================================================== --- head/lib/msun/src/s_cpow.c Sun Jul 15 00:23:10 2018 (r336299, copy source) +++ stable/11/lib/msun/src/s_cpow.c Fri Jul 27 17:39:36 2018 (r336767) @@ -49,6 +49,7 @@ __FBSDID("$FreeBSD$"); #include <complex.h> #include <float.h> #include <math.h> +#include "math_private.h" double complex cpow(double complex a, double complex z) @@ -60,7 +61,7 @@ cpow(double complex a, double complex z) y = cimag (z); absa = cabs (a); if (absa == 0.0) { - return (0.0 + 0.0 * I); + return (CMPLX(0.0, 0.0)); } arga = carg (a); r = pow (absa, x); @@ -69,6 +70,6 @@ cpow(double complex a, double complex z) r = r * exp (-y * arga); theta = theta + y * log (absa); } - w = r * cos (theta) + (r * sin (theta)) * I; + w = CMPLX(r * cos (theta), r * sin (theta)); return (w); } Copied and modified: stable/11/lib/msun/src/s_cpowf.c (from r336299, head/lib/msun/src/s_cpowf.c) ============================================================================== --- head/lib/msun/src/s_cpowf.c Sun Jul 15 00:23:10 2018 (r336299, copy source) +++ stable/11/lib/msun/src/s_cpowf.c Fri Jul 27 17:39:36 2018 (r336767) @@ -48,6 +48,7 @@ __FBSDID("$FreeBSD$"); #include <complex.h> #include <math.h> +#include "math_private.h" float complex cpowf(float complex a, float complex z) @@ -59,7 +60,7 @@ cpowf(float complex a, float complex z) y = cimagf(z); absa = cabsf (a); if (absa == 0.0f) { - return (0.0f + 0.0f * I); + return (CMPLXF(0.0f, 0.0f)); } arga = cargf (a); r = powf (absa, x); @@ -68,6 +69,6 @@ cpowf(float complex a, float complex z) r = r * expf (-y * arga); theta = theta + y * logf (absa); } - w = r * cosf (theta) + (r * sinf (theta)) * I; + w = CMPLXF(r * cosf (theta), r * sinf (theta)); return (w); } Copied and modified: stable/11/lib/msun/src/s_cpowl.c (from r336299, head/lib/msun/src/s_cpowl.c) ============================================================================== --- head/lib/msun/src/s_cpowl.c Sun Jul 15 00:23:10 2018 (r336299, copy source) +++ stable/11/lib/msun/src/s_cpowl.c Fri Jul 27 17:39:36 2018 (r336767) @@ -48,6 +48,7 @@ __FBSDID("$FreeBSD$"); #include <complex.h> #include <math.h> +#include "math_private.h" long double complex cpowl(long double complex a, long double complex z) @@ -59,7 +60,7 @@ cpowl(long double complex a, long double complex z) y = cimagl(z); absa = cabsl(a); if (absa == 0.0L) { - return (0.0L + 0.0L * I); + return (CMPLXL(0.0L, 0.0L)); } arga = cargl(a); r = powl(absa, x); @@ -68,6 +69,6 @@ cpowl(long double complex a, long double complex z) r = r * expl(-y * arga); theta = theta + y * logl(absa); } - w = r * cosl(theta) + (r * sinl(theta)) * I; + w = CMPLXL(r * cosl(theta), r * sinl(theta)); return (w); }
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201807271739.w6RHdaIZ087956>