diff options
author | sjg <sjg@FreeBSD.org> | 2015-05-27 01:19:58 +0000 |
---|---|---|
committer | sjg <sjg@FreeBSD.org> | 2015-05-27 01:19:58 +0000 |
commit | 65145fa4c81da358fcbc3b650156dab705dfa34e (patch) | |
tree | 55c065b6730aaac2afb6c29933ee6ec5fa4c4249 /lib/msun | |
parent | 60ff4eb0dff94a04d75d0d52a3957aaaf5f8c693 (diff) | |
parent | e6b664c390af88d4a87208bc042ce503da664c3b (diff) | |
download | FreeBSD-src-65145fa4c81da358fcbc3b650156dab705dfa34e.zip FreeBSD-src-65145fa4c81da358fcbc3b650156dab705dfa34e.tar.gz |
Merge sync of head
Diffstat (limited to 'lib/msun')
45 files changed, 649 insertions, 328 deletions
diff --git a/lib/msun/Makefile b/lib/msun/Makefile index afb7067..094fcba 100644 --- a/lib/msun/Makefile +++ b/lib/msun/Makefile @@ -223,6 +223,8 @@ MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3 .include <src.opts.mk> -.include <bsd.arch.inc.mk> +.if ${MK_TESTS} != "no" +SUBDIR+= tests +.endif .include <bsd.lib.mk> diff --git a/lib/msun/Makefile.amd64 b/lib/msun/Makefile.amd64 deleted file mode 100644 index dd0f5b0..0000000 --- a/lib/msun/Makefile.amd64 +++ /dev/null @@ -1,6 +0,0 @@ -# $FreeBSD$ - -.if ${MK_TESTS} != "no" -SUBDIR+= tests -.endif - diff --git a/lib/msun/Makefile.i386 b/lib/msun/Makefile.i386 deleted file mode 100644 index dd0f5b0..0000000 --- a/lib/msun/Makefile.i386 +++ /dev/null @@ -1,6 +0,0 @@ -# $FreeBSD$ - -.if ${MK_TESTS} != "no" -SUBDIR+= tests -.endif - diff --git a/lib/msun/aarch64/Makefile.inc b/lib/msun/aarch64/Makefile.inc new file mode 100644 index 0000000..286a608 --- /dev/null +++ b/lib/msun/aarch64/Makefile.inc @@ -0,0 +1,4 @@ +# $FreeBSD$ + +LDBL_PREC = 113 + diff --git a/lib/msun/aarch64/fenv.c b/lib/msun/aarch64/fenv.c new file mode 100644 index 0000000..4ec362f --- /dev/null +++ b/lib/msun/aarch64/fenv.c @@ -0,0 +1,56 @@ +/*- + * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG> + * Copyright (c) 2013 Andrew Turner <andrew@FreeBSD.ORG> + * 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. + * + * $FreeBSD$ + */ + +#define __fenv_static +#include "fenv.h" + +/* + * Hopefully the system ID byte is immutable, so it's valid to use + * this as a default environment. + */ +const fenv_t __fe_dfl_env = 0; + +#ifdef __GNUC_GNU_INLINE__ +#error "This file must be compiled with C99 'inline' semantics" +#endif + +extern inline int feclearexcept(int __excepts); +extern inline int fegetexceptflag(fexcept_t *__flagp, int __excepts); +extern inline int fesetexceptflag(const fexcept_t *__flagp, int __excepts); +extern inline int feraiseexcept(int __excepts); +extern inline int fetestexcept(int __excepts); +extern inline int fegetround(void); +extern inline int fesetround(int __round); +extern inline int fegetenv(fenv_t *__envp); +extern inline int feholdexcept(fenv_t *__envp); +extern inline int fesetenv(const fenv_t *__envp); +extern inline int feupdateenv(const fenv_t *__envp); +extern inline int feenableexcept(int __mask); +extern inline int fedisableexcept(int __mask); +extern inline int fegetexcept(void); diff --git a/lib/msun/aarch64/fenv.h b/lib/msun/aarch64/fenv.h new file mode 100644 index 0000000..2bd2987 --- /dev/null +++ b/lib/msun/aarch64/fenv.h @@ -0,0 +1,246 @@ +/*- + * Copyright (c) 2004-2005 David Schultz <das@FreeBSD.ORG> + * 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. + * + * $FreeBSD$ + */ + +#ifndef _FENV_H_ +#define _FENV_H_ + +#include <sys/_types.h> + +#ifndef __fenv_static +#define __fenv_static static +#endif + +typedef __uint64_t fenv_t; +typedef __uint64_t fexcept_t; + +/* Exception flags */ +#define FE_INVALID 0x00000001 +#define FE_DIVBYZERO 0x00000002 +#define FE_OVERFLOW 0x00000004 +#define FE_UNDERFLOW 0x00000008 +#define FE_INEXACT 0x00000010 +#define FE_ALL_EXCEPT (FE_DIVBYZERO | FE_INEXACT | \ + FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW) + +/* + * Rounding modes + * + * We can't just use the hardware bit values here, because that would + * make FE_UPWARD and FE_DOWNWARD negative, which is not allowed. + */ +#define FE_TONEAREST 0x0 +#define FE_UPWARD 0x1 +#define FE_DOWNWARD 0x2 +#define FE_TOWARDZERO 0x3 +#define _ROUND_MASK (FE_TONEAREST | FE_DOWNWARD | \ + FE_UPWARD | FE_TOWARDZERO) +#define _ROUND_SHIFT 22 + +__BEGIN_DECLS + +/* Default floating-point environment */ +extern const fenv_t __fe_dfl_env; +#define FE_DFL_ENV (&__fe_dfl_env) + +/* We need to be able to map status flag positions to mask flag positions */ +#define _FPUSW_SHIFT 8 +#define _ENABLE_MASK (FE_ALL_EXCEPT << _FPUSW_SHIFT) + +#define __mrs_fpcr(__r) __asm __volatile("mrs %0, fpcr" : "=r" (__r)) +#define __msr_fpcr(__r) __asm __volatile("msr fpcr, %0" : : "r" (__r)) + +#define __mrs_fpsr(__r) __asm __volatile("mrs %0, fpsr" : "=r" (__r)) +#define __msr_fpsr(__r) __asm __volatile("msr fpsr, %0" : : "r" (__r)) + +__fenv_static __inline int +feclearexcept(int __excepts) +{ + fexcept_t __r; + + __mrs_fpsr(__r); + __r &= ~__excepts; + __msr_fpsr(__r); + return (0); +} + +__fenv_static inline int +fegetexceptflag(fexcept_t *__flagp, int __excepts) +{ + fexcept_t __r; + + __mrs_fpsr(__r); + *__flagp = __r & __excepts; + return (0); +} + +__fenv_static inline int +fesetexceptflag(const fexcept_t *__flagp, int __excepts) +{ + fexcept_t __r; + + __mrs_fpsr(__r); + __r &= ~__excepts; + __r |= *__flagp & __excepts; + __msr_fpsr(__r); + return (0); +} + +__fenv_static inline int +feraiseexcept(int __excepts) +{ + fexcept_t __r; + + __mrs_fpsr(__r); + __r |= __excepts; + __msr_fpsr(__r); + return (0); +} + +__fenv_static inline int +fetestexcept(int __excepts) +{ + fexcept_t __r; + + __mrs_fpsr(__r); + return (__r & __excepts); +} + +__fenv_static inline int +fegetround(void) +{ + fenv_t __r; + + __mrs_fpcr(__r); + return ((__r >> _ROUND_SHIFT) & _ROUND_MASK); +} + +__fenv_static inline int +fesetround(int __round) +{ + fenv_t __r; + + if (__round & ~_ROUND_MASK) + return (-1); + __mrs_fpcr(__r); + __r &= ~(_ROUND_MASK << _ROUND_SHIFT); + __r |= __round << _ROUND_SHIFT; + __msr_fpcr(__r); + return (0); +} + +__fenv_static inline int +fegetenv(fenv_t *__envp) +{ + fenv_t __r; + + __mrs_fpcr(__r); + *__envp = __r & _ENABLE_MASK; + + __mrs_fpsr(__r); + *__envp |= __r & (FE_ALL_EXCEPT | (_ROUND_MASK << _ROUND_SHIFT)); + + return (0); +} + +__fenv_static inline int +feholdexcept(fenv_t *__envp) +{ + fenv_t __r; + + __mrs_fpcr(__r); + *__envp = __r & _ENABLE_MASK; + __r &= ~(_ENABLE_MASK); + __msr_fpcr(__r); + + __mrs_fpsr(__r); + *__envp |= __r & (FE_ALL_EXCEPT | (_ROUND_MASK << _ROUND_SHIFT)); + __r &= ~(_ENABLE_MASK); + __msr_fpsr(__r); + return (0); +} + +__fenv_static inline int +fesetenv(const fenv_t *__envp) +{ + + __msr_fpcr((*__envp) & _ENABLE_MASK); + __msr_fpsr((*__envp) & (FE_ALL_EXCEPT | (_ROUND_MASK << _ROUND_SHIFT))); + return (0); +} + +__fenv_static inline int +feupdateenv(const fenv_t *__envp) +{ + fexcept_t __r; + + __mrs_fpsr(__r); + fesetenv(__envp); + feraiseexcept(__r & FE_ALL_EXCEPT); + return (0); +} + +#if __BSD_VISIBLE + +/* We currently provide no external definitions of the functions below. */ + +static inline int +feenableexcept(int __mask) +{ + fenv_t __old_r, __new_r; + + __mrs_fpcr(__old_r); + __new_r = __old_r | ((__mask & FE_ALL_EXCEPT) << _FPUSW_SHIFT); + __msr_fpcr(__new_r); + return ((__old_r >> _FPUSW_SHIFT) & FE_ALL_EXCEPT); +} + +static inline int +fedisableexcept(int __mask) +{ + fenv_t __old_r, __new_r; + + __mrs_fpcr(__old_r); + __new_r = __old_r & ~((__mask & FE_ALL_EXCEPT) << _FPUSW_SHIFT); + __msr_fpcr(__new_r); + return ((__old_r >> _FPUSW_SHIFT) & FE_ALL_EXCEPT); +} + +static inline int +fegetexcept(void) +{ + fenv_t __r; + + __mrs_fpcr(__r); + return ((__r & _ENABLE_MASK) >> _FPUSW_SHIFT); +} + +#endif /* __BSD_VISIBLE */ + +__END_DECLS + +#endif /* !_FENV_H_ */ diff --git a/lib/msun/ld128/k_expl.h b/lib/msun/ld128/k_expl.h index a5668fd..e0a48fc 100644 --- a/lib/msun/ld128/k_expl.h +++ b/lib/msun/ld128/k_expl.h @@ -322,7 +322,7 @@ __ldexp_cexpl(long double complex z, int expt) scale2 = 1; SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt); - return (cpackl(cos(y) * exp_x * scale1 * scale2, + return (CMPLXL(cos(y) * exp_x * scale1 * scale2, sinl(y) * exp_x * scale1 * scale2)); } #endif /* _COMPLEX_H */ diff --git a/lib/msun/ld80/k_expl.h b/lib/msun/ld80/k_expl.h index ebfb9a8..9b081fa 100644 --- a/lib/msun/ld80/k_expl.h +++ b/lib/msun/ld80/k_expl.h @@ -299,7 +299,7 @@ __ldexp_cexpl(long double complex z, int expt) scale2 = 1; SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt); - return (cpackl(cos(y) * exp_x * scale1 * scale2, + return (CMPLXL(cos(y) * exp_x * scale1 * scale2, sinl(y) * exp_x * scale1 * scale2)); } #endif /* _COMPLEX_H */ diff --git a/lib/msun/man/cexp.3 b/lib/msun/man/cexp.3 index 97e36c1..776e6ce 100644 --- a/lib/msun/man/cexp.3 +++ b/lib/msun/man/cexp.3 @@ -103,7 +103,7 @@ is not finite, the sign of the result is indeterminate. .Sh SEE ALSO .Xr complex 3 , .Xr exp 3 , -.Xr math 3 , +.Xr math 3 .Sh STANDARDS The .Fn cexp diff --git a/lib/msun/man/complex.3 b/lib/msun/man/complex.3 index 34eb03e..d286494 100644 --- a/lib/msun/man/complex.3 +++ b/lib/msun/man/complex.3 @@ -103,9 +103,9 @@ ctan tangent ctanh hyperbolic tangent .El .Sh SEE ALSO -.Xr math 3 , .Xr fenv 3 , .Xr ieee 3 , +.Xr math 3 , .Xr tgmath 3 .Rs .%T "ISO/IEC 9899:TC3" diff --git a/lib/msun/man/csqrt.3 b/lib/msun/man/csqrt.3 index 1a1dfa0..d8066ae 100644 --- a/lib/msun/man/csqrt.3 +++ b/lib/msun/man/csqrt.3 @@ -85,7 +85,7 @@ an \*(Na is generated, an invalid exception will be thrown. .Sh SEE ALSO .Xr cabs 3 , .Xr fenv 3 , -.Xr math 3 , +.Xr math 3 .Sh STANDARDS The .Fn csqrt , diff --git a/lib/msun/man/j0.3 b/lib/msun/man/j0.3 index 91849da..587b72e 100644 --- a/lib/msun/man/j0.3 +++ b/lib/msun/man/j0.3 @@ -28,7 +28,7 @@ .\" from: @(#)j0.3 6.7 (Berkeley) 4/19/91 .\" $FreeBSD$ .\" -.Dd February 18, 2008 +.Dd March 10, 2015 .Dt J0 3 .Os .Sh NAME @@ -77,24 +77,17 @@ The functions .Fn j0 , .Fn j0f , -.Fn j1 +.Fn j1 , and .Fn j1f -compute the -.Em Bessel function of the first kind of the order -0 and the -.Em order -1, respectively, -for the -real value +compute the Bessel function of the first kind of orders +0 and 1 for the real value .Fa x ; the functions .Fn jn and .Fn jnf -compute the -.Em Bessel function of the first kind of the integer -.Em order +compute the Bessel function of the first kind of the integer order .Fa n for the real value .Fa x . @@ -105,13 +98,8 @@ The functions .Fn y1 , and .Fn y1f -compute the linearly independent -.Em Bessel function of the second kind of the order -0 and the -.Em order -1, respectively, -for the -positive +compute the linearly independent Bessel function of the second kind +of orders 0 and 1 for the positive .Em real value .Fa x ; @@ -119,9 +107,7 @@ the functions .Fn yn and .Fn ynf -compute the -.Em Bessel function of the second kind for the integer -.Em order +compute the Bessel function of the second kind for the integer order .Fa n for the positive .Em real @@ -141,11 +127,20 @@ and .Fn ynf . If .Fa x -is negative, these routines will generate an invalid exception and -return \*(Na. +is negative, including -\*(If, these routines will generate an invalid +exception and return \*(Na. +If +.Fa x +is \*(Pm0, these routines +will generate a divide-by-zero exception and return -\*(If. If .Fa x -is 0 or a sufficiently small positive number, these routines +is a sufficiently small positive number, then +.Fn y1 , +.Fn y1f , +.Fn yn , +and +.Fn ynf will generate an overflow exception and return -\*(If. .Sh SEE ALSO .Xr math 3 diff --git a/lib/msun/man/lgamma.3 b/lib/msun/man/lgamma.3 index d244722..c8a22a2 100644 --- a/lib/msun/man/lgamma.3 +++ b/lib/msun/man/lgamma.3 @@ -99,7 +99,7 @@ returns the sign of \(*G(x). and .Fn lgammal_r x signgamp provide the same functionality as -.Fn lgamma x , +.Fn lgamma x , .Fn lgammaf x , and .Fn lgammal x , @@ -124,7 +124,6 @@ are deprecated aliases for and .Fn lgammaf_r , respectively. - .Sh IDIOSYNCRASIES Do not use the expression .Dq Li signgam\(**exp(lgamma(x)) diff --git a/lib/msun/man/nextafter.3 b/lib/msun/man/nextafter.3 index c8c4aa3..01c472d 100644 --- a/lib/msun/man/nextafter.3 +++ b/lib/msun/man/nextafter.3 @@ -78,11 +78,11 @@ routines conform to They implement the Nextafter function recommended by .St -ieee754 , with the extension that -.Fn nextafter +0.0, -0.0 +.Fn nextafter "+0.0" "-0.0" returns .Li -0.0 , and -.Fn nextafter -0.0, +0.0 +.Fn nextafter "-0.0" "+0.0" returns .Li +0.0 . .Sh HISTORY diff --git a/lib/msun/man/sin.3 b/lib/msun/man/sin.3 index c7daf09..4a5352e 100644 --- a/lib/msun/man/sin.3 +++ b/lib/msun/man/sin.3 @@ -70,9 +70,9 @@ functions return the sine value. .Xr asin 3 , .Xr atan 3 , .Xr atan2 3 , -.Xr csin 3 , .Xr cos 3 , .Xr cosh 3 , +.Xr csin 3 , .Xr math 3 , .Xr sinh 3 , .Xr tan 3 , diff --git a/lib/msun/src/catrig.c b/lib/msun/src/catrig.c index 200977c..c0f5f55 100644 --- a/lib/msun/src/catrig.c +++ b/lib/msun/src/catrig.c @@ -286,19 +286,19 @@ casinh(double complex z) if (isnan(x) || isnan(y)) { /* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */ if (isinf(x)) - return (cpack(x, y + y)); + return (CMPLX(x, y + y)); /* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */ if (isinf(y)) - return (cpack(y, x + x)); + return (CMPLX(y, x + x)); /* casinh(NaN + I*0) = NaN + I*0 */ if (y == 0) - return (cpack(x + x, y)); + 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 (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { @@ -307,7 +307,7 @@ casinh(double complex z) w = clog_for_large_values(z) + m_ln2; else w = clog_for_large_values(-z) + m_ln2; - return (cpack(copysign(creal(w), x), copysign(cimag(w), y))); + return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y))); } /* Avoid spuriously raising inexact for z = 0. */ @@ -325,7 +325,7 @@ casinh(double complex z) ry = asin(B); else ry = atan2(new_y, sqrt_A2my2); - return (cpack(copysign(rx, x), copysign(ry, y))); + return (CMPLX(copysign(rx, x), copysign(ry, y))); } /* @@ -335,9 +335,9 @@ casinh(double complex z) double complex casin(double complex z) { - double complex w = casinh(cpack(cimag(z), creal(z))); + double complex w = casinh(CMPLX(cimag(z), creal(z))); - return (cpack(cimag(w), creal(w))); + return (CMPLX(cimag(w), creal(w))); } /* @@ -370,19 +370,19 @@ cacos(double complex z) if (isnan(x) || isnan(y)) { /* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */ if (isinf(x)) - return (cpack(y + y, -INFINITY)); + return (CMPLX(y + y, -INFINITY)); /* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */ if (isinf(y)) - return (cpack(x + x, -y)); + return (CMPLX(x + x, -y)); /* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */ if (x == 0) - return (cpack(pio2_hi + pio2_lo, y + y)); + 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 (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { @@ -392,18 +392,18 @@ cacos(double complex z) ry = creal(w) + m_ln2; if (sy == 0) ry = -ry; - return (cpack(rx, ry)); + return (CMPLX(rx, ry)); } /* Avoid spuriously raising inexact for z = 1. */ if (x == 1 && y == 0) - return (cpack(0, -y)); + return (CMPLX(0, -y)); /* All remaining cases are inexact. */ raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) - return (cpack(pio2_hi - (x - pio2_lo), -y)); + 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) { @@ -419,7 +419,7 @@ cacos(double complex z) } if (sy == 0) ry = -ry; - return (cpack(rx, ry)); + return (CMPLX(rx, ry)); } /* @@ -437,15 +437,15 @@ cacosh(double complex z) ry = cimag(w); /* cacosh(NaN + I*NaN) = NaN + I*NaN */ if (isnan(rx) && isnan(ry)) - return (cpack(ry, rx)); + return (CMPLX(ry, rx)); /* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */ /* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */ if (isnan(rx)) - return (cpack(fabs(ry), rx)); + return (CMPLX(fabs(ry), rx)); /* cacosh(0 + I*NaN) = NaN + I*NaN */ if (isnan(ry)) - return (cpack(ry, ry)); - return (cpack(fabs(ry), copysign(rx, cimag(z)))); + return (CMPLX(ry, ry)); + return (CMPLX(fabs(ry), copysign(rx, cimag(z)))); } /* @@ -475,16 +475,16 @@ clog_for_large_values(double complex z) * this method is still poor since it is uneccessarily slow. */ if (ax > DBL_MAX / 2) - return (cpack(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x))); + 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 (cpack(log(hypot(x, y)), atan2(y, x))); + return (CMPLX(log(hypot(x, y)), atan2(y, x))); - return (cpack(log(ax * ax + ay * ay) / 2, atan2(y, x))); + return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x))); } /* @@ -575,30 +575,30 @@ catanh(double complex z) /* This helps handle many cases. */ if (y == 0 && ax <= 1) - return (cpack(atanh(x), y)); + return (CMPLX(atanh(x), y)); /* To ensure the same accuracy as atan(), and to filter out z = 0. */ if (x == 0) - return (cpack(x, atan(y))); + return (CMPLX(x, atan(y))); if (isnan(x) || isnan(y)) { /* catanh(+-Inf + I*NaN) = +-0 + I*NaN */ if (isinf(x)) - return (cpack(copysign(0, x), y + y)); + return (CMPLX(copysign(0, x), y + y)); /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */ if (isinf(y)) - return (cpack(copysign(0, x), + 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 (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) - return (cpack(real_part_reciprocal(x, y), + return (CMPLX(real_part_reciprocal(x, y), copysign(pio2_hi + pio2_lo, y))); if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { @@ -623,7 +623,7 @@ catanh(double complex z) else ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; - return (cpack(copysign(rx, x), copysign(ry, y))); + return (CMPLX(copysign(rx, x), copysign(ry, y))); } /* @@ -633,7 +633,7 @@ catanh(double complex z) double complex catan(double complex z) { - double complex w = catanh(cpack(cimag(z), creal(z))); + double complex w = catanh(CMPLX(cimag(z), creal(z))); - return (cpack(cimag(w), creal(w))); + return (CMPLX(cimag(w), creal(w))); } diff --git a/lib/msun/src/catrigf.c b/lib/msun/src/catrigf.c index 08ebef7..ca84ce2 100644 --- a/lib/msun/src/catrigf.c +++ b/lib/msun/src/catrigf.c @@ -156,12 +156,12 @@ casinhf(float complex z) if (isnan(x) || isnan(y)) { if (isinf(x)) - return (cpackf(x, y + y)); + return (CMPLXF(x, y + y)); if (isinf(y)) - return (cpackf(y, x + x)); + return (CMPLXF(y, x + x)); if (y == 0) - return (cpackf(x + x, y)); - return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXF(x + x, y)); + return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { @@ -169,7 +169,7 @@ casinhf(float complex z) w = clog_for_large_values(z) + m_ln2; else w = clog_for_large_values(-z) + m_ln2; - return (cpackf(copysignf(crealf(w), x), + return (CMPLXF(copysignf(crealf(w), x), copysignf(cimagf(w), y))); } @@ -186,15 +186,15 @@ casinhf(float complex z) ry = asinf(B); else ry = atan2f(new_y, sqrt_A2my2); - return (cpackf(copysignf(rx, x), copysignf(ry, y))); + return (CMPLXF(copysignf(rx, x), copysignf(ry, y))); } float complex casinf(float complex z) { - float complex w = casinhf(cpackf(cimagf(z), crealf(z))); + float complex w = casinhf(CMPLXF(cimagf(z), crealf(z))); - return (cpackf(cimagf(w), crealf(w))); + return (CMPLXF(cimagf(w), crealf(w))); } float complex @@ -214,12 +214,12 @@ cacosf(float complex z) if (isnan(x) || isnan(y)) { if (isinf(x)) - return (cpackf(y + y, -INFINITY)); + return (CMPLXF(y + y, -INFINITY)); if (isinf(y)) - return (cpackf(x + x, -y)); + return (CMPLXF(x + x, -y)); if (x == 0) - return (cpackf(pio2_hi + pio2_lo, y + y)); - return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXF(pio2_hi + pio2_lo, y + y)); + return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { @@ -228,16 +228,16 @@ cacosf(float complex z) ry = crealf(w) + m_ln2; if (sy == 0) ry = -ry; - return (cpackf(rx, ry)); + return (CMPLXF(rx, ry)); } if (x == 1 && y == 0) - return (cpackf(0, -y)); + return (CMPLXF(0, -y)); raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) - return (cpackf(pio2_hi - (x - pio2_lo), -y)); + 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) { @@ -253,7 +253,7 @@ cacosf(float complex z) } if (sy == 0) ry = -ry; - return (cpackf(rx, ry)); + return (CMPLXF(rx, ry)); } float complex @@ -266,12 +266,12 @@ cacoshf(float complex z) rx = crealf(w); ry = cimagf(w); if (isnan(rx) && isnan(ry)) - return (cpackf(ry, rx)); + return (CMPLXF(ry, rx)); if (isnan(rx)) - return (cpackf(fabsf(ry), rx)); + return (CMPLXF(fabsf(ry), rx)); if (isnan(ry)) - return (cpackf(ry, ry)); - return (cpackf(fabsf(ry), copysignf(rx, cimagf(z)))); + return (CMPLXF(ry, ry)); + return (CMPLXF(fabsf(ry), copysignf(rx, cimagf(z)))); } static float complex @@ -291,13 +291,13 @@ clog_for_large_values(float complex z) } if (ax > FLT_MAX / 2) - return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1, + return (CMPLXF(logf(hypotf(x / m_e, y / m_e)) + 1, atan2f(y, x))); if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) - return (cpackf(logf(hypotf(x, y)), atan2f(y, x))); + return (CMPLXF(logf(hypotf(x, y)), atan2f(y, x))); - return (cpackf(logf(ax * ax + ay * ay) / 2, atan2f(y, x))); + return (CMPLXF(logf(ax * ax + ay * ay) / 2, atan2f(y, x))); } static inline float @@ -346,22 +346,22 @@ catanhf(float complex z) ay = fabsf(y); if (y == 0 && ax <= 1) - return (cpackf(atanhf(x), y)); + return (CMPLXF(atanhf(x), y)); if (x == 0) - return (cpackf(x, atanf(y))); + return (CMPLXF(x, atanf(y))); if (isnan(x) || isnan(y)) { if (isinf(x)) - return (cpackf(copysignf(0, x), y + y)); + return (CMPLXF(copysignf(0, x), y + y)); if (isinf(y)) - return (cpackf(copysignf(0, x), + return (CMPLXF(copysignf(0, x), copysignf(pio2_hi + pio2_lo, y))); - return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) - return (cpackf(real_part_reciprocal(x, y), + return (CMPLXF(real_part_reciprocal(x, y), copysignf(pio2_hi + pio2_lo, y))); if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { @@ -381,13 +381,13 @@ catanhf(float complex z) else ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; - return (cpackf(copysignf(rx, x), copysignf(ry, y))); + return (CMPLXF(copysignf(rx, x), copysignf(ry, y))); } float complex catanf(float complex z) { - float complex w = catanhf(cpackf(cimagf(z), crealf(z))); + float complex w = catanhf(CMPLXF(cimagf(z), crealf(z))); - return (cpackf(cimagf(w), crealf(w))); + return (CMPLXF(cimagf(w), crealf(w))); } diff --git a/lib/msun/src/e_j0.c b/lib/msun/src/e_j0.c index 8320f25..a253ed1 100644 --- a/lib/msun/src/e_j0.c +++ b/lib/msun/src/e_j0.c @@ -62,7 +62,9 @@ __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" -static double pzero(double), qzero(double); +static __inline double pzero(double), qzero(double); + +static const volatile double vone = 1, vzero = 0; static const double huge = 1e300, @@ -115,7 +117,7 @@ __ieee754_j0(double x) if(ix<0x3f200000) { /* |x| < 2**-13 */ if(huge+x>one) { /* raise inexact if x != 0 */ if(ix<0x3e400000) return one; /* |x|<2**-27 */ - else return one - 0.25*x*x; + else return one - x*x/4; } } z = x*x; @@ -150,10 +152,16 @@ __ieee754_y0(double x) EXTRACT_WORDS(hx,lx,x); ix = 0x7fffffff&hx; - /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ - if(ix>=0x7ff00000) return one/(x+x*x); - if((ix|lx)==0) return -one/zero; - if(hx<0) return zero/zero; + /* + * y0(NaN) = NaN. + * y0(Inf) = 0. + * y0(-Inf) = NaN and raise invalid exception. + */ + if(ix>=0x7ff00000) return vone/(x+x*x); + /* y0(+-0) = -inf and raise divide-by-zero exception. */ + if((ix|lx)==0) return -one/vzero; + /* y0(x<0) = NaN and raise invalid exception. */ + if(hx<0) return vzero/vzero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) * where x0 = x-pi/4 @@ -268,7 +276,8 @@ static const double pS2[5] = { 1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */ }; - static double pzero(double x) +static __inline double +pzero(double x) { const double *p,*q; double z,r,s; @@ -278,7 +287,7 @@ static const double pS2[5] = { if(ix>=0x40200000) {p = pR8; q= pS8;} else if(ix>=0x40122E8B){p = pR5; q= pS5;} else if(ix>=0x4006DB6D){p = pR3; q= pS3;} - else if(ix>=0x40000000){p = pR2; q= pS2;} + else {p = pR2; q= pS2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); @@ -363,7 +372,8 @@ static const double qS2[6] = { -5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */ }; - static double qzero(double x) +static __inline double +qzero(double x) { const double *p,*q; double s,r,z; @@ -373,7 +383,7 @@ static const double qS2[6] = { if(ix>=0x40200000) {p = qR8; q= qS8;} else if(ix>=0x40122E8B){p = qR5; q= qS5;} else if(ix>=0x4006DB6D){p = qR3; q= qS3;} - else if(ix>=0x40000000){p = qR2; q= qS2;} + else {p = qR2; q= qS2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); diff --git a/lib/msun/src/e_j0f.c b/lib/msun/src/e_j0f.c index c45faf3..3e2f7e8 100644 --- a/lib/msun/src/e_j0f.c +++ b/lib/msun/src/e_j0f.c @@ -16,10 +16,16 @@ #include <sys/cdefs.h> __FBSDID("$FreeBSD$"); +/* + * See e_j0.c for complete comments. + */ + #include "math.h" #include "math_private.h" -static float pzerof(float), qzerof(float); +static __inline float pzerof(float), qzerof(float); + +static const volatile float vone = 1, vzero = 0; static const float huge = 1e30, @@ -62,17 +68,17 @@ __ieee754_j0f(float x) * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) */ - if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x); + if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(x); /* |x|>2**49 */ else { u = pzerof(x); v = qzerof(x); z = invsqrtpi*(u*cc-v*ss)/sqrtf(x); } return z; } - if(ix<0x39000000) { /* |x| < 2**-13 */ + if(ix<0x3b000000) { /* |x| < 2**-9 */ if(huge+x>one) { /* raise inexact if x != 0 */ - if(ix<0x32000000) return one; /* |x|<2**-27 */ - else return one - (float)0.25*x*x; + if(ix<0x39800000) return one; /* |x|<2**-12 */ + else return one - x*x/4; } } z = x*x; @@ -107,10 +113,9 @@ __ieee754_y0f(float x) GET_FLOAT_WORD(hx,x); ix = 0x7fffffff&hx; - /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ - if(ix>=0x7f800000) return one/(x+x*x); - if(ix==0) return -one/zero; - if(hx<0) return zero/zero; + if(ix>=0x7f800000) return vone/(x+x*x); + if(ix==0) return -one/vzero; + if(hx<0) return vzero/vzero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) * where x0 = x-pi/4 @@ -136,14 +141,14 @@ __ieee754_y0f(float x) if ((s*c)<zero) cc = z/ss; else ss = z/cc; } - if(ix>0x80000000) z = (invsqrtpi*ss)/sqrtf(x); + if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */ else { u = pzerof(x); v = qzerof(x); z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); } return z; } - if(ix<=0x32000000) { /* x < 2**-27 */ + if(ix<=0x39000000) { /* x < 2**-13 */ return(u00 + tpi*__ieee754_logf(x)); } z = x*x; @@ -224,7 +229,8 @@ static const float pS2[5] = { 1.4657617569e+01, /* 0x416a859a */ }; - static float pzerof(float x) +static __inline float +pzerof(float x) { const float *p,*q; float z,r,s; @@ -232,9 +238,9 @@ static const float pS2[5] = { GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; if(ix>=0x41000000) {p = pR8; q= pS8;} - else if(ix>=0x40f71c58){p = pR5; q= pS5;} - else if(ix>=0x4036db68){p = pR3; q= pS3;} - else if(ix>=0x40000000){p = pR2; q= pS2;} + else if(ix>=0x409173eb){p = pR5; q= pS5;} + else if(ix>=0x4036d917){p = pR3; q= pS3;} + else {p = pR2; q= pS2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); @@ -319,7 +325,8 @@ static const float qS2[6] = { -5.3109550476e+00, /* 0xc0a9f358 */ }; - static float qzerof(float x) +static __inline float +qzerof(float x) { const float *p,*q; float s,r,z; @@ -327,9 +334,9 @@ static const float qS2[6] = { GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; if(ix>=0x41000000) {p = qR8; q= qS8;} - else if(ix>=0x40f71c58){p = qR5; q= qS5;} - else if(ix>=0x4036db68){p = qR3; q= qS3;} - else if(ix>=0x40000000){p = qR2; q= qS2;} + else if(ix>=0x409173eb){p = qR5; q= qS5;} + else if(ix>=0x4036d917){p = qR3; q= qS3;} + else {p = qR2; q= qS2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); diff --git a/lib/msun/src/e_j1.c b/lib/msun/src/e_j1.c index 63800ad..74a7c69 100644 --- a/lib/msun/src/e_j1.c +++ b/lib/msun/src/e_j1.c @@ -62,7 +62,9 @@ __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" -static double pone(double), qone(double); +static __inline double pone(double), qone(double); + +static const volatile double vone = 1, vzero = 0; static const double huge = 1e300, @@ -147,10 +149,16 @@ __ieee754_y1(double x) EXTRACT_WORDS(hx,lx,x); ix = 0x7fffffff&hx; - /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if(ix>=0x7ff00000) return one/(x+x*x); - if((ix|lx)==0) return -one/zero; - if(hx<0) return zero/zero; + /* + * y1(NaN) = NaN. + * y1(Inf) = 0. + * y1(-Inf) = NaN and raise invalid exception. + */ + if(ix>=0x7ff00000) return vone/(x+x*x); + /* y1(+-0) = -inf and raise divide-by-zero exception. */ + if((ix|lx)==0) return -one/vzero; + /* y1(x<0) = NaN and raise invalid exception. */ + if(hx<0) return vzero/vzero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ s = sin(x); c = cos(x); @@ -262,7 +270,8 @@ static const double ps2[5] = { 8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */ }; - static double pone(double x) +static __inline double +pone(double x) { const double *p,*q; double z,r,s; @@ -272,7 +281,7 @@ static const double ps2[5] = { if(ix>=0x40200000) {p = pr8; q= ps8;} else if(ix>=0x40122E8B){p = pr5; q= ps5;} else if(ix>=0x4006DB6D){p = pr3; q= ps3;} - else if(ix>=0x40000000){p = pr2; q= ps2;} + else {p = pr2; q= ps2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); @@ -358,7 +367,8 @@ static const double qs2[6] = { -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */ }; - static double qone(double x) +static __inline double +qone(double x) { const double *p,*q; double s,r,z; @@ -368,7 +378,7 @@ static const double qs2[6] = { if(ix>=0x40200000) {p = qr8; q= qs8;} else if(ix>=0x40122E8B){p = qr5; q= qs5;} else if(ix>=0x4006DB6D){p = qr3; q= qs3;} - else if(ix>=0x40000000){p = qr2; q= qs2;} + else {p = qr2; q= qs2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); diff --git a/lib/msun/src/e_j1f.c b/lib/msun/src/e_j1f.c index 88e2d83..ec7f381 100644 --- a/lib/msun/src/e_j1f.c +++ b/lib/msun/src/e_j1f.c @@ -16,10 +16,16 @@ #include <sys/cdefs.h> __FBSDID("$FreeBSD$"); +/* + * See e_j1.c for complete comments. + */ + #include "math.h" #include "math_private.h" -static float ponef(float), qonef(float); +static __inline float ponef(float), qonef(float); + +static const volatile float vone = 1, vzero = 0; static const float huge = 1e30, @@ -63,7 +69,7 @@ __ieee754_j1f(float x) * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) */ - if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(y); + if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(y); /* |x|>2**49 */ else { u = ponef(y); v = qonef(y); z = invsqrtpi*(u*cc-v*ss)/sqrtf(y); @@ -71,7 +77,7 @@ __ieee754_j1f(float x) if(hx<0) return -z; else return z; } - if(ix<0x32000000) { /* |x|<2**-27 */ + if(ix<0x39000000) { /* |x|<2**-13 */ if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */ } z = x*x; @@ -104,10 +110,9 @@ __ieee754_y1f(float x) GET_FLOAT_WORD(hx,x); ix = 0x7fffffff&hx; - /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if(ix>=0x7f800000) return one/(x+x*x); - if(ix==0) return -one/zero; - if(hx<0) return zero/zero; + if(ix>=0x7f800000) return vone/(x+x*x); + if(ix==0) return -one/vzero; + if(hx<0) return vzero/vzero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ s = sinf(x); c = cosf(x); @@ -129,14 +134,14 @@ __ieee754_y1f(float x) * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one. */ - if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x); + if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */ else { u = ponef(x); v = qonef(x); z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); } return z; } - if(ix<=0x24800000) { /* x < 2**-54 */ + if(ix<=0x33000000) { /* x < 2**-25 */ return(-tpi/x); } z = x*x; @@ -219,7 +224,8 @@ static const float ps2[5] = { 8.3646392822e+00, /* 0x4105d590 */ }; - static float ponef(float x) +static __inline float +ponef(float x) { const float *p,*q; float z,r,s; @@ -227,9 +233,9 @@ static const float ps2[5] = { GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; if(ix>=0x41000000) {p = pr8; q= ps8;} - else if(ix>=0x40f71c58){p = pr5; q= ps5;} - else if(ix>=0x4036db68){p = pr3; q= ps3;} - else if(ix>=0x40000000){p = pr2; q= ps2;} + else if(ix>=0x409173eb){p = pr5; q= ps5;} + else if(ix>=0x4036d917){p = pr3; q= ps3;} + else {p = pr2; q= ps2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); @@ -315,17 +321,18 @@ static const float qs2[6] = { -4.9594988823e+00, /* 0xc09eb437 */ }; - static float qonef(float x) +static __inline float +qonef(float x) { const float *p,*q; float s,r,z; int32_t ix; GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; - if(ix>=0x40200000) {p = qr8; q= qs8;} - else if(ix>=0x40f71c58){p = qr5; q= qs5;} - else if(ix>=0x4036db68){p = qr3; q= qs3;} - else if(ix>=0x40000000){p = qr2; q= qs2;} + if(ix>=0x41000000) {p = qr8; q= qs8;} + else if(ix>=0x409173eb){p = qr5; q= qs5;} + else if(ix>=0x4036d917){p = qr3; q= qs3;} + else {p = qr2; q= qs2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); diff --git a/lib/msun/src/e_jn.c b/lib/msun/src/e_jn.c index 8b0bc62..eefd4ff 100644 --- a/lib/msun/src/e_jn.c +++ b/lib/msun/src/e_jn.c @@ -43,6 +43,8 @@ __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" +static const volatile double vone = 1, vzero = 0; + static const double invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */ @@ -220,10 +222,12 @@ __ieee754_yn(int n, double x) EXTRACT_WORDS(hx,lx,x); ix = 0x7fffffff&hx; - /* if Y(n,NaN) is NaN */ + /* yn(n,NaN) = NaN */ if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x; - if((ix|lx)==0) return -one/zero; - if(hx<0) return zero/zero; + /* yn(n,+-0) = -inf and raise divide-by-zero exception. */ + if((ix|lx)==0) return -one/vzero; + /* yn(n,x<0) = NaN and raise invalid exception. */ + if(hx<0) return vzero/vzero; sign = 1; if(n<0){ n = -n; diff --git a/lib/msun/src/e_jnf.c b/lib/msun/src/e_jnf.c index f564aec..965feeb 100644 --- a/lib/msun/src/e_jnf.c +++ b/lib/msun/src/e_jnf.c @@ -16,9 +16,15 @@ #include <sys/cdefs.h> __FBSDID("$FreeBSD$"); +/* + * See e_jn.c for complete comments. + */ + #include "math.h" #include "math_private.h" +static const volatile float vone = 1, vzero = 0; + static const float two = 2.0000000000e+00, /* 0x40000000 */ one = 1.0000000000e+00; /* 0x3F800000 */ @@ -172,10 +178,9 @@ __ieee754_ynf(int n, float x) GET_FLOAT_WORD(hx,x); ix = 0x7fffffff&hx; - /* if Y(n,NaN) is NaN */ if(ix>0x7f800000) return x+x; - if(ix==0) return -one/zero; - if(hx<0) return zero/zero; + if(ix==0) return -one/vzero; + if(hx<0) return vzero/vzero; sign = 1; if(n<0){ n = -n; diff --git a/lib/msun/src/k_exp.c b/lib/msun/src/k_exp.c index f592f69..ecef54c 100644 --- a/lib/msun/src/k_exp.c +++ b/lib/msun/src/k_exp.c @@ -103,6 +103,6 @@ __ldexp_cexp(double complex z, int expt) half_expt = expt - half_expt; INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); - return (cpack(cos(y) * exp_x * scale1 * scale2, + return (CMPLX(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2)); } diff --git a/lib/msun/src/k_expf.c b/lib/msun/src/k_expf.c index 548a008..f8c254d 100644 --- a/lib/msun/src/k_expf.c +++ b/lib/msun/src/k_expf.c @@ -82,6 +82,6 @@ __ldexp_cexpf(float complex z, int expt) half_expt = expt - half_expt; SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23); - return (cpackf(cosf(y) * exp_x * scale1 * scale2, + return (CMPLXF(cosf(y) * exp_x * scale1 * scale2, sinf(y) * exp_x * scale1 * scale2)); } diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h index 8af2c65..afaf201 100644 --- a/lib/msun/src/math_private.h +++ b/lib/msun/src/math_private.h @@ -454,9 +454,15 @@ typedef union { * (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 -cpackf(float x, float y) +CMPLXF(float x, float y) { float_complex z; @@ -464,9 +470,11 @@ cpackf(float x, float y) IMAGPART(z) = y; return (z.f); } +#endif +#ifndef CMPLX static __inline double complex -cpack(double x, double y) +CMPLX(double x, double y) { double_complex z; @@ -474,9 +482,11 @@ cpack(double x, double y) IMAGPART(z) = y; return (z.f); } +#endif +#ifndef CMPLXL static __inline long double complex -cpackl(long double x, long double y) +CMPLXL(long double x, long double y) { long_double_complex z; @@ -484,6 +494,8 @@ cpackl(long double x, long double y) IMAGPART(z) = y; return (z.f); } +#endif + #endif /* _COMPLEX_H */ #ifdef __GNUCLIKE_ASM diff --git a/lib/msun/src/s_ccosh.c b/lib/msun/src/s_ccosh.c index 9ea962b..fbe15fc 100644 --- a/lib/msun/src/s_ccosh.c +++ b/lib/msun/src/s_ccosh.c @@ -62,23 +62,23 @@ ccosh(double complex z) /* Handle the nearly-non-exceptional cases where x and y are finite. */ if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) - return (cpack(cosh(x), x * y)); + return (CMPLX(cosh(x), x * y)); if (ix < 0x40360000) /* small x: normal case */ - return (cpack(cosh(x) * cos(y), sinh(x) * sin(y))); + 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 (cpack(h * cos(y), copysign(h, x) * sin(y))); + return (CMPLX(h * cos(y), copysign(h, x) * sin(y))); } else if (ix < 0x4096bbaa) { /* x < 1455: scale to avoid overflow */ - z = __ldexp_cexp(cpack(fabs(x), y), -1); - return (cpack(creal(z), cimag(z) * copysign(1, x))); + 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 (cpack(h * h * cos(y), h * sin(y))); + return (CMPLX(h * h * cos(y), h * sin(y))); } } @@ -92,7 +92,7 @@ ccosh(double complex z) * the same as d(NaN). */ if ((ix | lx) == 0 && iy >= 0x7ff00000) - return (cpack(y - y, copysign(0, x * (y - y)))); + return (CMPLX(y - y, copysign(0, x * (y - y)))); /* * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0. @@ -102,8 +102,8 @@ ccosh(double complex z) */ if ((iy | ly) == 0 && ix >= 0x7ff00000) { if (((hx & 0xfffff) | lx) == 0) - return (cpack(x * x, copysign(0, x) * y)); - return (cpack(x * x, copysign(0, (x + x) * y))); + return (CMPLX(x * x, copysign(0, x) * y)); + return (CMPLX(x * x, copysign(0, (x + x) * y))); } /* @@ -115,7 +115,7 @@ ccosh(double complex z) * nonzero x. Choice = don't raise (except for signaling NaNs). */ if (ix < 0x7ff00000 && iy >= 0x7ff00000) - return (cpack(y - y, x * (y - y))); + return (CMPLX(y - y, x * (y - y))); /* * cosh(+-Inf + I NaN) = +Inf + I d(NaN). @@ -128,8 +128,8 @@ ccosh(double complex z) */ if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { if (iy >= 0x7ff00000) - return (cpack(x * x, x * (y - y))); - return (cpack((x * x) * cos(y), x * sin(y))); + return (CMPLX(x * x, x * (y - y))); + return (CMPLX((x * x) * cos(y), x * sin(y))); } /* @@ -143,7 +143,7 @@ ccosh(double complex z) * Optionally raises the invalid floating-point exception for finite * nonzero y. Choice = don't raise (except for signaling NaNs). */ - return (cpack((x * x) * (y - y), (x + x) * (y - y))); + return (CMPLX((x * x) * (y - y), (x + x) * (y - y))); } double complex @@ -151,5 +151,5 @@ ccos(double complex z) { /* ccos(z) = ccosh(I * z) */ - return (ccosh(cpack(-cimag(z), creal(z)))); + return (ccosh(CMPLX(-cimag(z), creal(z)))); } diff --git a/lib/msun/src/s_ccoshf.c b/lib/msun/src/s_ccoshf.c index 1de9ad4..fe8cf89 100644 --- a/lib/msun/src/s_ccoshf.c +++ b/lib/msun/src/s_ccoshf.c @@ -55,50 +55,50 @@ ccoshf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) - return (cpackf(coshf(x), x * y)); + return (CMPLXF(coshf(x), x * y)); if (ix < 0x41100000) /* small x: normal case */ - return (cpackf(coshf(x) * cosf(y), sinhf(x) * sinf(y))); + 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 (cpackf(h * cosf(y), copysignf(h, x) * sinf(y))); + return (CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y))); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ - z = __ldexp_cexpf(cpackf(fabsf(x), y), -1); - return (cpackf(crealf(z), cimagf(z) * copysignf(1, x))); + 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 (cpackf(h * h * cosf(y), h * sinf(y))); + return (CMPLXF(h * h * cosf(y), h * sinf(y))); } } if (ix == 0 && iy >= 0x7f800000) - return (cpackf(y - y, copysignf(0, x * (y - y)))); + return (CMPLXF(y - y, copysignf(0, x * (y - y)))); if (iy == 0 && ix >= 0x7f800000) { if ((hx & 0x7fffff) == 0) - return (cpackf(x * x, copysignf(0, x) * y)); - return (cpackf(x * x, copysignf(0, (x + x) * y))); + return (CMPLXF(x * x, copysignf(0, x) * y)); + return (CMPLXF(x * x, copysignf(0, (x + x) * y))); } if (ix < 0x7f800000 && iy >= 0x7f800000) - return (cpackf(y - y, x * (y - y))); + return (CMPLXF(y - y, x * (y - y))); if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) { if (iy >= 0x7f800000) - return (cpackf(x * x, x * (y - y))); - return (cpackf((x * x) * cosf(y), x * sinf(y))); + return (CMPLXF(x * x, x * (y - y))); + return (CMPLXF((x * x) * cosf(y), x * sinf(y))); } - return (cpackf((x * x) * (y - y), (x + x) * (y - y))); + return (CMPLXF((x * x) * (y - y), (x + x) * (y - y))); } float complex ccosf(float complex z) { - return (ccoshf(cpackf(-cimagf(z), crealf(z)))); + return (ccoshf(CMPLXF(-cimagf(z), crealf(z)))); } diff --git a/lib/msun/src/s_cexp.c b/lib/msun/src/s_cexp.c index abe178f..9e11d51 100644 --- a/lib/msun/src/s_cexp.c +++ b/lib/msun/src/s_cexp.c @@ -50,22 +50,22 @@ cexp(double complex z) /* cexp(x + I 0) = exp(x) + I 0 */ if ((hy | ly) == 0) - return (cpack(exp(x), y)); + return (CMPLX(exp(x), y)); EXTRACT_WORDS(hx, lx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ if (((hx & 0x7fffffff) | lx) == 0) - return (cpack(cos(y), sin(y))); + return (CMPLX(cos(y), sin(y))); if (hy >= 0x7ff00000) { if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ - return (cpack(y - y, y - y)); + return (CMPLX(y - y, y - y)); } else if (hx & 0x80000000) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ - return (cpack(0.0, 0.0)); + return (CMPLX(0.0, 0.0)); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ - return (cpack(x, y - y)); + return (CMPLX(x, y - y)); } } @@ -84,6 +84,6 @@ cexp(double complex z) * - x = NaN (spurious inexact exception from y) */ exp_x = exp(x); - return (cpack(exp_x * cos(y), exp_x * sin(y))); + return (CMPLX(exp_x * cos(y), exp_x * sin(y))); } } diff --git a/lib/msun/src/s_cexpf.c b/lib/msun/src/s_cexpf.c index 0e30d08..0138cd1 100644 --- a/lib/msun/src/s_cexpf.c +++ b/lib/msun/src/s_cexpf.c @@ -50,22 +50,22 @@ cexpf(float complex z) /* cexp(x + I 0) = exp(x) + I 0 */ if (hy == 0) - return (cpackf(expf(x), y)); + return (CMPLXF(expf(x), y)); GET_FLOAT_WORD(hx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ if ((hx & 0x7fffffff) == 0) - return (cpackf(cosf(y), sinf(y))); + return (CMPLXF(cosf(y), sinf(y))); if (hy >= 0x7f800000) { if ((hx & 0x7fffffff) != 0x7f800000) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ - return (cpackf(y - y, y - y)); + return (CMPLXF(y - y, y - y)); } else if (hx & 0x80000000) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ - return (cpackf(0.0, 0.0)); + return (CMPLXF(0.0, 0.0)); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ - return (cpackf(x, y - y)); + return (CMPLXF(x, y - y)); } } @@ -84,6 +84,6 @@ cexpf(float complex z) * - x = NaN (spurious inexact exception from y) */ exp_x = expf(x); - return (cpackf(exp_x * cosf(y), exp_x * sinf(y))); + return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y))); } } diff --git a/lib/msun/src/s_conj.c b/lib/msun/src/s_conj.c index 5770c29..d0dd051 100644 --- a/lib/msun/src/s_conj.c +++ b/lib/msun/src/s_conj.c @@ -34,5 +34,5 @@ double complex conj(double complex z) { - return (cpack(creal(z), -cimag(z))); + return (CMPLX(creal(z), -cimag(z))); } diff --git a/lib/msun/src/s_conjf.c b/lib/msun/src/s_conjf.c index b090760..c4cf127 100644 --- a/lib/msun/src/s_conjf.c +++ b/lib/msun/src/s_conjf.c @@ -34,5 +34,5 @@ float complex conjf(float complex z) { - return (cpackf(crealf(z), -cimagf(z))); + return (CMPLXF(crealf(z), -cimagf(z))); } diff --git a/lib/msun/src/s_conjl.c b/lib/msun/src/s_conjl.c index 0e431ef..a4604b6 100644 --- a/lib/msun/src/s_conjl.c +++ b/lib/msun/src/s_conjl.c @@ -34,5 +34,5 @@ long double complex conjl(long double complex z) { - return (cpackl(creall(z), -cimagl(z))); + return (CMPLXL(creall(z), -cimagl(z))); } diff --git a/lib/msun/src/s_cproj.c b/lib/msun/src/s_cproj.c index 8e9404c..2f0db6e 100644 --- a/lib/msun/src/s_cproj.c +++ b/lib/msun/src/s_cproj.c @@ -39,7 +39,7 @@ cproj(double complex z) if (!isinf(creal(z)) && !isinf(cimag(z))) return (z); else - return (cpack(INFINITY, copysign(0.0, cimag(z)))); + return (CMPLX(INFINITY, copysign(0.0, cimag(z)))); } #if LDBL_MANT_DIG == 53 diff --git a/lib/msun/src/s_cprojf.c b/lib/msun/src/s_cprojf.c index 68ea77b..bd27bdc 100644 --- a/lib/msun/src/s_cprojf.c +++ b/lib/msun/src/s_cprojf.c @@ -39,5 +39,5 @@ cprojf(float complex z) if (!isinf(crealf(z)) && !isinf(cimagf(z))) return (z); else - return (cpackf(INFINITY, copysignf(0.0, cimagf(z)))); + return (CMPLXF(INFINITY, copysignf(0.0, cimagf(z)))); } diff --git a/lib/msun/src/s_cprojl.c b/lib/msun/src/s_cprojl.c index 07385bc..5cccc91 100644 --- a/lib/msun/src/s_cprojl.c +++ b/lib/msun/src/s_cprojl.c @@ -39,5 +39,5 @@ cprojl(long double complex z) if (!isinf(creall(z)) && !isinf(cimagl(z))) return (z); else - return (cpackl(INFINITY, copysignl(0.0, cimagl(z)))); + return (CMPLXL(INFINITY, copysignl(0.0, cimagl(z)))); } diff --git a/lib/msun/src/s_csinh.c b/lib/msun/src/s_csinh.c index c192f30..37f46da 100644 --- a/lib/msun/src/s_csinh.c +++ b/lib/msun/src/s_csinh.c @@ -62,23 +62,23 @@ csinh(double complex z) /* Handle the nearly-non-exceptional cases where x and y are finite. */ if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) - return (cpack(sinh(x), y)); + return (CMPLX(sinh(x), y)); if (ix < 0x40360000) /* small x: normal case */ - return (cpack(sinh(x) * cos(y), cosh(x) * sin(y))); + 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 (cpack(copysign(h, x) * cos(y), h * sin(y))); + return (CMPLX(copysign(h, x) * cos(y), h * sin(y))); } else if (ix < 0x4096bbaa) { /* x < 1455: scale to avoid overflow */ - z = __ldexp_cexp(cpack(fabs(x), y), -1); - return (cpack(creal(z) * copysign(1, x), cimag(z))); + 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 (cpack(h * cos(y), h * h * sin(y))); + return (CMPLX(h * cos(y), h * h * sin(y))); } } @@ -92,7 +92,7 @@ csinh(double complex z) * the same as d(NaN). */ if ((ix | lx) == 0 && iy >= 0x7ff00000) - return (cpack(copysign(0, x * (y - y)), y - y)); + return (CMPLX(copysign(0, x * (y - y)), y - y)); /* * sinh(+-Inf +- I 0) = +-Inf + I +-0. @@ -101,8 +101,8 @@ csinh(double complex z) */ if ((iy | ly) == 0 && ix >= 0x7ff00000) { if (((hx & 0xfffff) | lx) == 0) - return (cpack(x, y)); - return (cpack(x, copysign(0, y))); + return (CMPLX(x, y)); + return (CMPLX(x, copysign(0, y))); } /* @@ -114,7 +114,7 @@ csinh(double complex z) * nonzero x. Choice = don't raise (except for signaling NaNs). */ if (ix < 0x7ff00000 && iy >= 0x7ff00000) - return (cpack(y - y, x * (y - y))); + return (CMPLX(y - y, x * (y - y))); /* * sinh(+-Inf + I NaN) = +-Inf + I d(NaN). @@ -129,8 +129,8 @@ csinh(double complex z) */ if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { if (iy >= 0x7ff00000) - return (cpack(x * x, x * (y - y))); - return (cpack(x * cos(y), INFINITY * sin(y))); + return (CMPLX(x * x, x * (y - y))); + return (CMPLX(x * cos(y), INFINITY * sin(y))); } /* @@ -144,7 +144,7 @@ csinh(double complex z) * Optionally raises the invalid floating-point exception for finite * nonzero y. Choice = don't raise (except for signaling NaNs). */ - return (cpack((x * x) * (y - y), (x + x) * (y - y))); + return (CMPLX((x * x) * (y - y), (x + x) * (y - y))); } double complex @@ -152,6 +152,6 @@ csin(double complex z) { /* csin(z) = -I * csinh(I * z) */ - z = csinh(cpack(-cimag(z), creal(z))); - return (cpack(cimag(z), -creal(z))); + z = csinh(CMPLX(-cimag(z), creal(z))); + return (CMPLX(cimag(z), -creal(z))); } diff --git a/lib/msun/src/s_csinhf.c b/lib/msun/src/s_csinhf.c index c523125..6348272 100644 --- a/lib/msun/src/s_csinhf.c +++ b/lib/msun/src/s_csinhf.c @@ -55,51 +55,51 @@ csinhf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) - return (cpackf(sinhf(x), y)); + return (CMPLXF(sinhf(x), y)); if (ix < 0x41100000) /* small x: normal case */ - return (cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y))); + 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 (cpackf(copysignf(h, x) * cosf(y), h * sinf(y))); + return (CMPLXF(copysignf(h, x) * cosf(y), h * sinf(y))); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ - z = __ldexp_cexpf(cpackf(fabsf(x), y), -1); - return (cpackf(crealf(z) * copysignf(1, x), cimagf(z))); + 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 (cpackf(h * cosf(y), h * h * sinf(y))); + return (CMPLXF(h * cosf(y), h * h * sinf(y))); } } if (ix == 0 && iy >= 0x7f800000) - return (cpackf(copysignf(0, x * (y - y)), y - y)); + return (CMPLXF(copysignf(0, x * (y - y)), y - y)); if (iy == 0 && ix >= 0x7f800000) { if ((hx & 0x7fffff) == 0) - return (cpackf(x, y)); - return (cpackf(x, copysignf(0, y))); + return (CMPLXF(x, y)); + return (CMPLXF(x, copysignf(0, y))); } if (ix < 0x7f800000 && iy >= 0x7f800000) - return (cpackf(y - y, x * (y - y))); + return (CMPLXF(y - y, x * (y - y))); if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) { if (iy >= 0x7f800000) - return (cpackf(x * x, x * (y - y))); - return (cpackf(x * cosf(y), INFINITY * sinf(y))); + return (CMPLXF(x * x, x * (y - y))); + return (CMPLXF(x * cosf(y), INFINITY * sinf(y))); } - return (cpackf((x * x) * (y - y), (x + x) * (y - y))); + return (CMPLXF((x * x) * (y - y), (x + x) * (y - y))); } float complex csinf(float complex z) { - z = csinhf(cpackf(-cimagf(z), crealf(z))); - return (cpackf(cimagf(z), -crealf(z))); + z = csinhf(CMPLXF(-cimagf(z), crealf(z))); + return (CMPLXF(cimagf(z), -crealf(z))); } diff --git a/lib/msun/src/s_csqrt.c b/lib/msun/src/s_csqrt.c index 18a7ae3..75e12d8 100644 --- a/lib/msun/src/s_csqrt.c +++ b/lib/msun/src/s_csqrt.c @@ -58,12 +58,12 @@ csqrt(double complex z) /* Handle special cases. */ if (z == 0) - return (cpack(0, b)); + return (CMPLX(0, b)); if (isinf(b)) - return (cpack(INFINITY, b)); + return (CMPLX(INFINITY, b)); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ - return (cpack(a, t)); /* return NaN + NaN i */ + return (CMPLX(a, t)); /* return NaN + NaN i */ } if (isinf(a)) { /* @@ -73,9 +73,9 @@ csqrt(double complex z) * csqrt(-inf + y i) = 0 + inf i */ if (signbit(a)) - return (cpack(fabs(b - b), copysign(a, b))); + return (CMPLX(fabs(b - b), copysign(a, b))); else - return (cpack(a, copysign(b - b, b))); + return (CMPLX(a, copysign(b - b, b))); } /* * The remaining special case (b is NaN) is handled just fine by @@ -94,10 +94,10 @@ csqrt(double complex z) /* Algorithm 312, CACM vol 10, Oct 1967. */ if (a >= 0) { t = sqrt((a + hypot(a, b)) * 0.5); - result = cpack(t, b / (2 * t)); + result = CMPLX(t, b / (2 * t)); } else { t = sqrt((-a + hypot(a, b)) * 0.5); - result = cpack(fabs(b) / (2 * t), copysign(t, b)); + result = CMPLX(fabs(b) / (2 * t), copysign(t, b)); } /* Rescale. */ diff --git a/lib/msun/src/s_csqrtf.c b/lib/msun/src/s_csqrtf.c index da7fe18..7ee513f 100644 --- a/lib/msun/src/s_csqrtf.c +++ b/lib/msun/src/s_csqrtf.c @@ -49,12 +49,12 @@ csqrtf(float complex z) /* Handle special cases. */ if (z == 0) - return (cpackf(0, b)); + return (CMPLXF(0, b)); if (isinf(b)) - return (cpackf(INFINITY, b)); + return (CMPLXF(INFINITY, b)); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ - return (cpackf(a, t)); /* return NaN + NaN i */ + return (CMPLXF(a, t)); /* return NaN + NaN i */ } if (isinf(a)) { /* @@ -64,9 +64,9 @@ csqrtf(float complex z) * csqrtf(-inf + y i) = 0 + inf i */ if (signbit(a)) - return (cpackf(fabsf(b - b), copysignf(a, b))); + return (CMPLXF(fabsf(b - b), copysignf(a, b))); else - return (cpackf(a, copysignf(b - b, b))); + return (CMPLXF(a, copysignf(b - b, b))); } /* * The remaining special case (b is NaN) is handled just fine by @@ -80,9 +80,9 @@ csqrtf(float complex z) */ if (a >= 0) { t = sqrt((a + hypot(a, b)) * 0.5); - return (cpackf(t, b / (2.0 * t))); + return (CMPLXF(t, b / (2.0 * t))); } else { t = sqrt((-a + hypot(a, b)) * 0.5); - return (cpackf(fabsf(b) / (2.0 * t), copysignf(t, b))); + return (CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b))); } } diff --git a/lib/msun/src/s_csqrtl.c b/lib/msun/src/s_csqrtl.c index dd18e1e..ea2ef27 100644 --- a/lib/msun/src/s_csqrtl.c +++ b/lib/msun/src/s_csqrtl.c @@ -58,12 +58,12 @@ csqrtl(long double complex z) /* Handle special cases. */ if (z == 0) - return (cpackl(0, b)); + return (CMPLXL(0, b)); if (isinf(b)) - return (cpackl(INFINITY, b)); + return (CMPLXL(INFINITY, b)); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ - return (cpackl(a, t)); /* return NaN + NaN i */ + return (CMPLXL(a, t)); /* return NaN + NaN i */ } if (isinf(a)) { /* @@ -73,9 +73,9 @@ csqrtl(long double complex z) * csqrt(-inf + y i) = 0 + inf i */ if (signbit(a)) - return (cpackl(fabsl(b - b), copysignl(a, b))); + return (CMPLXL(fabsl(b - b), copysignl(a, b))); else - return (cpackl(a, copysignl(b - b, b))); + return (CMPLXL(a, copysignl(b - b, b))); } /* * The remaining special case (b is NaN) is handled just fine by @@ -94,10 +94,10 @@ csqrtl(long double complex z) /* Algorithm 312, CACM vol 10, Oct 1967. */ if (a >= 0) { t = sqrtl((a + hypotl(a, b)) * 0.5); - result = cpackl(t, b / (2 * t)); + result = CMPLXL(t, b / (2 * t)); } else { t = sqrtl((-a + hypotl(a, b)) * 0.5); - result = cpackl(fabsl(b) / (2 * t), copysignl(t, b)); + result = CMPLXL(fabsl(b) / (2 * t), copysignl(t, b)); } /* Rescale. */ diff --git a/lib/msun/src/s_ctanh.c b/lib/msun/src/s_ctanh.c index d427e28..b15c814 100644 --- a/lib/msun/src/s_ctanh.c +++ b/lib/msun/src/s_ctanh.c @@ -102,9 +102,9 @@ ctanh(double complex z) */ if (ix >= 0x7ff00000) { if ((ix & 0xfffff) | lx) /* x is NaN */ - return (cpack(x, (y == 0 ? y : x * y))); + return (CMPLX(x, (y == 0 ? y : x * y))); SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */ - return (cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); + return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); } /* @@ -112,7 +112,7 @@ ctanh(double complex z) * ctanh(x +- i Inf) = NaN + i NaN */ if (!isfinite(y)) - return (cpack(y - y, y - y)); + return (CMPLX(y - y, y - y)); /* * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the @@ -121,7 +121,7 @@ ctanh(double complex z) */ if (ix >= 0x40360000) { /* x >= 22 */ double exp_mx = exp(-fabs(x)); - return (cpack(copysign(1, x), + return (CMPLX(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx)); } @@ -131,7 +131,7 @@ ctanh(double complex z) s = sinh(x); rho = sqrt(1 + s * s); /* = cosh(x) */ denom = 1 + beta * s * s; - return (cpack((beta * rho * s) / denom, t / denom)); + return (CMPLX((beta * rho * s) / denom, t / denom)); } double complex @@ -139,6 +139,6 @@ ctan(double complex z) { /* ctan(z) = -I * ctanh(I * z) */ - z = ctanh(cpack(-cimag(z), creal(z))); - return (cpack(cimag(z), -creal(z))); + z = ctanh(CMPLX(-cimag(z), creal(z))); + return (CMPLX(cimag(z), -creal(z))); } diff --git a/lib/msun/src/s_ctanhf.c b/lib/msun/src/s_ctanhf.c index 4be28d8..9b0cda1 100644 --- a/lib/msun/src/s_ctanhf.c +++ b/lib/msun/src/s_ctanhf.c @@ -51,18 +51,18 @@ ctanhf(float complex z) if (ix >= 0x7f800000) { if (ix & 0x7fffff) - return (cpackf(x, (y == 0 ? y : x * y))); + return (CMPLXF(x, (y == 0 ? y : x * y))); SET_FLOAT_WORD(x, hx - 0x40000000); - return (cpackf(x, + return (CMPLXF(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y)))); } if (!isfinite(y)) - return (cpackf(y - y, y - y)); + return (CMPLXF(y - y, y - y)); if (ix >= 0x41300000) { /* x >= 11 */ float exp_mx = expf(-fabsf(x)); - return (cpackf(copysignf(1, x), + return (CMPLXF(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx)); } @@ -71,14 +71,14 @@ ctanhf(float complex z) s = sinhf(x); rho = sqrtf(1 + s * s); denom = 1 + beta * s * s; - return (cpackf((beta * rho * s) / denom, t / denom)); + return (CMPLXF((beta * rho * s) / denom, t / denom)); } float complex ctanf(float complex z) { - z = ctanhf(cpackf(-cimagf(z), crealf(z))); - return (cpackf(cimagf(z), -crealf(z))); + z = ctanhf(CMPLXF(-cimagf(z), crealf(z))); + return (CMPLXF(cimagf(z), -crealf(z))); } diff --git a/lib/msun/src/s_scalbln.c b/lib/msun/src/s_scalbln.c index d609d4e..5e4e42e 100644 --- a/lib/msun/src/s_scalbln.c +++ b/lib/msun/src/s_scalbln.c @@ -27,50 +27,28 @@ #include <sys/cdefs.h> __FBSDID("$FreeBSD$"); -#include <limits.h> #include <math.h> +#define NMAX 65536 +#define NMIN -65536 + double -scalbln (double x, long n) +scalbln(double x, long n) { - int in; - in = (int)n; - if (in != n) { - if (n > 0) - in = INT_MAX; - else - in = INT_MIN; - } - return (scalbn(x, in)); + return (scalbn(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n)); } float -scalblnf (float x, long n) +scalblnf(float x, long n) { - int in; - in = (int)n; - if (in != n) { - if (n > 0) - in = INT_MAX; - else - in = INT_MIN; - } - return (scalbnf(x, in)); + return (scalbnf(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n)); } long double -scalblnl (long double x, long n) +scalblnl(long double x, long n) { - int in; - in = (int)n; - if (in != n) { - if (n > 0) - in = INT_MAX; - else - in = INT_MIN; - } - return (scalbnl(x, (int)n)); + return (scalbnl(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n)); } diff --git a/lib/msun/tests/Makefile b/lib/msun/tests/Makefile index 4261e48..0479cfb 100644 --- a/lib/msun/tests/Makefile +++ b/lib/msun/tests/Makefile @@ -6,13 +6,12 @@ TESTSRC= ${SRCTOP}/contrib/netbsd-tests/lib/libm TESTSDIR= ${TESTSBASE}/lib/msun -.if ${MACHINE} == "sparc" || ${MACHINE} == "i386" \ - || ${MACHINE} == "amd64" || ${MACHINE_CPU} == "arm" \ - || ${MACHINE} == "sparc64" +# All architectures on FreeBSD have fenv.h CFLAGS+= -DHAVE_FENV_H -.endif -.if ${MACHINE} == "amd64" || ${MACHINE} == "i386" +# Not sure why this isn't defined for all architectures, since most +# have long double. +.if ${MACHINE_CPUARCH} == "amd64" || ${MACHINE_CPUARCH} == "i386" CFLAGS+= -D__HAVE_LONG_DOUBLE .endif @@ -41,8 +40,7 @@ NETBSD_ATF_TESTS_C+= tanh_test CSTD= c99 -LDADD+= -lm -DPADD+= ${LIBM} +LIBADD+= m #COPTS+= -Wfloat-equal # Copied from lib/msun/Makefile |