diff --git a/include/complex.h b/include/complex.h index 892bc55e5145..c31c15d9da4b 100644 --- a/include/complex.h +++ b/include/complex.h @@ -1,138 +1,140 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2001-2011 The FreeBSD Project. * 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 _COMPLEX_H #define _COMPLEX_H #include #ifdef __GNUC__ #if __STDC_VERSION__ < 199901 #define _Complex __complex__ #endif #define _Complex_I ((float _Complex)1.0i) #endif #ifdef __generic _Static_assert(__generic(_Complex_I, float _Complex, 1, 0), "_Complex_I must be of type float _Complex"); #endif #define complex _Complex #define I _Complex_I #if __ISO_C_VISIBLE >= 2011 #ifdef __clang__ #define CMPLX(x, y) ((double complex){ x, y }) #define CMPLXF(x, y) ((float complex){ x, y }) #define CMPLXL(x, y) ((long double complex){ x, y }) #elif __GNUC_PREREQ__(4, 7) #define CMPLX(x, y) __builtin_complex((double)(x), (double)(y)) #define CMPLXF(x, y) __builtin_complex((float)(x), (float)(y)) #define CMPLXL(x, y) __builtin_complex((long double)(x), (long double)(y)) #endif #endif /* __ISO_C_VISIBLE >= 2011 */ __BEGIN_DECLS double cabs(double complex); float cabsf(float complex); long double cabsl(long double complex); double complex cacos(double complex); float complex cacosf(float complex); double complex cacosh(double complex); float complex cacoshf(float complex); long double complex cacoshl(long double complex); long double complex cacosl(long double complex); double carg(double complex); float cargf(float complex); long double cargl(long double complex); double complex casin(double complex); float complex casinf(float complex); double complex casinh(double complex); float complex casinhf(float complex); long double complex casinhl(long double complex); long double complex casinl(long double complex); double complex catan(double complex); float complex catanf(float complex); double complex catanh(double complex); float complex catanhf(float complex); long double complex catanhl(long double complex); long double complex catanl(long double complex); double complex ccos(double complex); float complex ccosf(float complex); double complex ccosh(double complex); float complex ccoshf(float complex); double complex cexp(double complex); float complex cexpf(float complex); +long double complex + cexpl(long double complex); double cimag(double complex) __pure2; float cimagf(float complex) __pure2; long double cimagl(long double complex) __pure2; double complex clog(double complex); float complex clogf(float complex); long double complex clogl(long double complex); 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 cprojl(long double complex) __pure2; double creal(double complex) __pure2; float crealf(float complex) __pure2; long double creall(long double complex) __pure2; double complex csin(double complex); float complex csinf(float complex); double complex csinh(double complex); float complex csinhf(float complex); double complex csqrt(double complex); float complex csqrtf(float complex); long double complex csqrtl(long double complex); double complex ctan(double complex); float complex ctanf(float complex); double complex ctanh(double complex); float complex ctanhf(float complex); __END_DECLS #endif /* _COMPLEX_H */ diff --git a/lib/msun/Makefile b/lib/msun/Makefile index 1d94e371e61f..d7c0e2f88358 100644 --- a/lib/msun/Makefile +++ b/lib/msun/Makefile @@ -1,273 +1,273 @@ # @(#)Makefile 5.1beta 93/09/24 # $FreeBSD$ # # ==================================================== # 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. # ==================================================== # # PACKAGE= clibs ARCH_SUBDIR= ${MACHINE_CPUARCH:S/i386/i387/} .include "${ARCH_SUBDIR}/Makefile.inc" .PATH: ${.CURDIR}/${ARCH_SUBDIR} .if ${MACHINE_CPUARCH} == "i386" || ${MACHINE_CPUARCH} == "amd64" .PATH: ${.CURDIR}/x86 CFLAGS+= -I${.CURDIR}/x86 .endif # long double format .if ${LDBL_PREC} == 64 .PATH: ${.CURDIR}/ld80 CFLAGS+= -I${.CURDIR}/ld80 .elif ${LDBL_PREC} == 113 .PATH: ${.CURDIR}/ld128 CFLAGS+= -I${.CURDIR}/ld128 .endif CFLAGS+= -I${.CURDIR}/${ARCH_SUBDIR} .include .if ${COMPILER_TYPE} == "clang" && ${COMPILER_VERSION} >= 100000 && \ (${MACHINE_CPUARCH} == "amd64" || ${MACHINE_CPUARCH} == "i386") # When using clang with x86_64 CPUs that support AVX, some floating point # transformations may raise exceptions that would not have been raised by the # original code. To avoid this, use the -fp-exception-behavior=maytrap flag, # introduced in clang 10.0.0. # See also: https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=254911 CFLAGS+= -ffp-exception-behavior=maytrap .endif .PATH: ${.CURDIR}/bsdsrc .PATH: ${.CURDIR}/src .PATH: ${.CURDIR}/man LIB= m SHLIBDIR?= /lib SHLIB_MAJOR= 5 WARNS?= 1 IGNORE_PRAGMA= COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \ e_acos.c e_acosf.c e_acosh.c e_acoshf.c e_asin.c e_asinf.c \ e_atan2.c e_atan2f.c e_atanh.c e_atanhf.c e_cosh.c e_coshf.c e_exp.c \ e_expf.c e_fmod.c e_fmodf.c e_gamma.c e_gamma_r.c e_gammaf.c \ e_gammaf_r.c e_hypot.c e_hypotf.c e_j0.c e_j0f.c e_j1.c e_j1f.c \ e_jn.c e_jnf.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \ e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \ e_pow.c e_powf.c e_rem_pio2.c \ e_rem_pio2f.c e_remainder.c e_remainderf.c e_scalb.c e_scalbf.c \ e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \ imprecise.c \ k_cos.c k_cosf.c k_exp.c k_expf.c k_rem_pio2.c k_sin.c k_sinf.c \ k_tan.c k_tanf.c \ s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c s_cargl.c \ s_cbrt.c s_cbrtf.c s_ceil.c s_ceilf.c s_clog.c s_clogf.c \ s_copysign.c s_copysignf.c s_cos.c s_cosf.c \ s_csqrt.c s_csqrtf.c s_erf.c s_erff.c \ s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c \ s_finite.c s_finitef.c \ s_floor.c s_floorf.c s_fma.c s_fmaf.c \ s_fmax.c s_fmaxf.c s_fmin.c \ s_fminf.c s_frexp.c s_frexpf.c s_ilogb.c s_ilogbf.c \ s_ilogbl.c s_isfinite.c s_isnan.c s_isnormal.c \ s_llrint.c s_llrintf.c s_llround.c s_llroundf.c s_llroundl.c \ s_log1p.c s_log1pf.c s_logb.c s_logbf.c s_lrint.c s_lrintf.c \ s_lround.c s_lroundf.c s_lroundl.c s_modff.c \ s_nan.c s_nearbyint.c s_nextafter.c s_nextafterf.c \ s_nexttowardf.c s_remquo.c s_remquof.c \ s_rint.c s_rintf.c s_round.c s_roundf.c \ s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \ s_signgam.c s_significand.c s_significandf.c s_sin.c \ s_sincos.c s_sincosf.c s_sinf.c \ s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \ w_cabs.c w_cabsf.c w_drem.c w_dremf.c # IEEE-754 2008 and ISO/IEC TS 18661-4 half-cycle trignometric functions COMMON_SRCS+= s_cospi.c s_cospif.c \ s_sinpi.c s_sinpif.c \ s_tanpi.c s_tanpif.c # Location of fpmath.h and _fpmath.h .if exists(${LIBCSRCDIR}/${MACHINE_ARCH}) LIBC_ARCH=${MACHINE_ARCH} .else LIBC_ARCH=${MACHINE_CPUARCH} .endif CFLAGS+= -I${.CURDIR}/src -I${LIBCSRCDIR}/include \ -I${LIBCSRCDIR}/${LIBC_ARCH} SYM_MAPS+= ${.CURDIR}/Symbol.map VERSION_DEF= ${LIBCSRCDIR}/Versions.def SYMBOL_MAPS= ${SYM_MAPS} # C99 long double functions COMMON_SRCS+= s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c .if ${LDBL_PREC} != 53 # If long double != double use these; otherwise, we alias the double versions. 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_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 \ + s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cexpl.c \ s_clogl.c s_cosl.c s_cospil.c s_cprojl.c \ s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \ s_fmaxl.c s_fminl.c s_frexpl.c s_logbl.c s_logl.c s_nanl.c \ s_nextafterl.c s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c \ s_scalbnl.c s_sinl.c s_sincosl.c s_sinpil.c \ s_tanhl.c s_tanl.c s_tanpil.c s_truncl.c w_cabsl.c # Work around this warning from gcc: # lib/msun/ld80/e_powl.c:275:1: error: floating constant exceeds range of # 'long double' [-Werror=overflow] # if( y >= LDBL_MAX ) # See also: https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=130067 .if ${COMPILER_TYPE} == "gcc" CFLAGS.e_powl.c+= -Wno-error=overflow .endif .endif # C99 complex functions 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 # FreeBSD's C library supplies these functions: #COMMON_SRCS+= s_fabs.c s_frexp.c s_isnan.c s_ldexp.c s_modf.c # Exclude the generic versions of what we provide in the MD area. .if defined(ARCH_SRCS) .for i in ${ARCH_SRCS} COMMON_SRCS:= ${COMMON_SRCS:N${i:R}.c} .endfor .endif SRCS= ${COMMON_SRCS} ${ARCH_SRCS} 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 cospi.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 \ fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \ lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \ nextafter.3 remainder.3 rint.3 \ round.3 scalbn.3 signbit.3 sin.3 sincos.3 \ sinh.3 sinpi.3 sqrt.3 tan.3 tanh.3 tanpi.3 trunc.3 \ complex.3 MLINKS+=acos.3 acosf.3 acos.3 acosl.3 MLINKS+=acosh.3 acoshf.3 acosh.3 acoshl.3 MLINKS+=asin.3 asinf.3 asin.3 asinl.3 MLINKS+=asinh.3 asinhf.3 asinh.3 asinhl.3 MLINKS+=atan.3 atanf.3 atan.3 atanl.3 MLINKS+=atanh.3 atanhf.3 atanh.3 atanhl.3 MLINKS+=atan2.3 atan2f.3 atan2.3 atan2l.3 \ atan2.3 carg.3 atan2.3 cargf.3 atan2.3 cargl.3 MLINKS+=cacos.3 cacosf.3 cacos.3 cacosl.3 \ cacos.3 cacosh.3 cacos.3 cacoshf.3 cacos.3 cacoshl.3 \ cacos.3 casin.3 cacos.3 casinf.3 cacos.3 casinl.3 \ cacos.3 casinh.3 cacos.3 casinhf.3 cacos.3 casinhl.3 \ cacos.3 catan.3 cacos.3 catanf.3 cacos.3 catanl.3 \ cacos.3 catanh.3 cacos.3 catanhf.3 cacos.3 catanhl.3 MLINKS+=ccos.3 ccosf.3 ccos.3 csin.3 ccos.3 csinf.3 ccos.3 ctan.3 ccos.3 ctanf.3 MLINKS+=ccosh.3 ccoshf.3 ccosh.3 csinh.3 ccosh.3 csinhf.3 \ ccosh.3 ctanh.3 ccosh.3 ctanhf.3 MLINKS+=ceil.3 ceilf.3 ceil.3 ceill.3 -MLINKS+=cexp.3 cexpf.3 +MLINKS+=cexp.3 cexpf.3 cexp.3 cexpl.3 MLINKS+=cimag.3 cimagf.3 cimag.3 cimagl.3 \ cimag.3 conj.3 cimag.3 conjf.3 cimag.3 conjl.3 \ cimag.3 cproj.3 cimag.3 cprojf.3 cimag.3 cprojl.3 \ cimag.3 creal.3 cimag.3 crealf.3 cimag.3 creall.3 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+=cospi.3 cospif.3 cospi.3 cospil.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 \ exp.3 powl.3 exp.3 exp2.3 exp.3 exp2f.3 exp.3 exp2l.3 exp.3 expf.3 \ exp.3 expl.3 MLINKS+=fabs.3 fabsf.3 fabs.3 fabsl.3 MLINKS+=fdim.3 fdimf.3 fdim.3 fdiml.3 MLINKS+=feclearexcept.3 fegetexceptflag.3 feclearexcept.3 feraiseexcept.3 \ feclearexcept.3 fesetexceptflag.3 feclearexcept.3 fetestexcept.3 MLINKS+=feenableexcept.3 fedisableexcept.3 feenableexcept.3 fegetexcept.3 MLINKS+=fegetenv.3 feholdexcept.3 fegetenv.3 fesetenv.3 \ fegetenv.3 feupdateenv.3 MLINKS+=fegetround.3 fesetround.3 MLINKS+=floor.3 floorf.3 floor.3 floorl.3 MLINKS+=fma.3 fmaf.3 fma.3 fmal.3 MLINKS+=fmax.3 fmaxf.3 fmax.3 fmaxl.3 \ fmax.3 fmin.3 fmax.3 fminf.3 fmax.3 fminl.3 MLINKS+=fmod.3 fmodf.3 fmod.3 fmodl.3 MLINKS+=hypot.3 cabs.3 hypot.3 cabsf.3 hypot.3 cabsl.3 \ hypot.3 hypotf.3 hypot.3 hypotl.3 MLINKS+=ieee_test.3 scalb.3 ieee_test.3 scalbf.3 MLINKS+=ieee_test.3 significand.3 ieee_test.3 significandf.3 MLINKS+=ilogb.3 ilogbf.3 ilogb.3 ilogbl.3 \ ilogb.3 logb.3 ilogb.3 logbf.3 ilogb.3 logbl.3 MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3 MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3 MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 \ lgamma.3 lgammaf.3 lgamma.3 lgammal.3 \ lgamma.3 tgamma.3 lgamma.3 tgammaf.3 MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log10l.3 \ log.3 log1p.3 log.3 log1pf.3 log.3 log1pl.3 \ log.3 logf.3 log.3 logl.3 \ log.3 log2.3 log.3 log2f.3 log.3 log2l.3 MLINKS+=lrint.3 llrint.3 lrint.3 llrintf.3 lrint.3 llrintl.3 \ lrint.3 lrintf.3 lrint.3 lrintl.3 MLINKS+=lround.3 llround.3 lround.3 llroundf.3 lround.3 llroundl.3 \ lround.3 lroundf.3 lround.3 lroundl.3 MLINKS+=nan.3 nanf.3 nan.3 nanl.3 MLINKS+=nextafter.3 nextafterf.3 nextafter.3 nextafterl.3 MLINKS+=nextafter.3 nexttoward.3 nextafter.3 nexttowardf.3 MLINKS+=nextafter.3 nexttowardl.3 MLINKS+=remainder.3 remainderf.3 remainder.3 remainderl.3 \ remainder.3 remquo.3 remainder.3 remquof.3 remainder.3 remquol.3 MLINKS+=rint.3 rintf.3 rint.3 rintl.3 \ rint.3 nearbyint.3 rint.3 nearbyintf.3 rint.3 nearbyintl.3 MLINKS+=round.3 roundf.3 round.3 roundl.3 MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3 MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3 MLINKS+=sin.3 sinf.3 sin.3 sinl.3 MLINKS+=sincos.3 sincosf.3 sin.3 sincosl.3 MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3 MLINKS+=sinpi.3 sinpif.3 sinpi.3 sinpil.3 MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \ sqrt.3 sqrtl.3 MLINKS+=tan.3 tanf.3 tan.3 tanl.3 MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3 MLINKS+=tanpi.3 tanpif.3 tanpi.3 tanpil.3 MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3 .include HAS_TESTS= SUBDIR.${MK_TESTS}+= tests .include .if ${COMPILER_TYPE} == "clang" && ${COMPILER_VERSION} >= 120000 # Silence '#pragma FENV_ACCESS' is not supported on this target - ignored CWARNFLAGS+= -Wno-error=ignored-pragmas .endif .include diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map index 7229e7ef31fd..8650d56eb9d5 100644 --- a/lib/msun/Symbol.map +++ b/lib/msun/Symbol.map @@ -1,319 +1,320 @@ /* * $FreeBSD$ */ /* 7.0-CURRENT */ FBSD_1.0 { __fe_dfl_env; tgamma; acos; acosf; acosh; acoshf; asin; asinf; atan2; atan2f; atanh; atanhf; cosh; coshf; exp; expf; fmod; fmodf; gamma; gamma_r; gammaf; gammaf_r; hypot; hypotf; j0; y0; j0f; y0f; j1; y1; j1f; y1f; jn; yn; jnf; ynf; lgamma; lgamma_r; lgammaf; lgammaf_r; log; log10; log10f; logf; pow; powf; remainder; remainderf; scalb; scalbf; sinh; sinhf; sqrt; sqrtf; asinh; asinhf; atan; atanf; cbrt; cbrtf; ceil; ceilf; ceill; cimag; cimagf; cimagl; conj; conjf; conjl; copysign; copysignf; copysignl; cos; cosf; creal; crealf; creall; erf; erfc; erff; erfcf; exp2; exp2f; expm1; expm1f; fabs; fabsf; fabsl; fdim; fdimf; fdiml; finite; finitef; floor; floorf; floorl; fma; fmaf; fmal; fmax; fmaxf; fmaxl; fmin; fminf; fminl; frexp; frexpf; frexpl; ilogb; ilogbf; ilogbl; __isfinite; __isfinitef; __isfinitel; isnanf; __isnanl; __isnormal; __isnormalf; __isnormall; llrint; llrintf; llround; llroundf; llroundl; log1p; log1pf; logb; logbf; lrint; lrintf; lround; lroundf; lroundl; modff; modfl; nearbyint; nearbyintf; nextafter; nexttoward; nexttowardl; nextafterl; nextafterf; nexttowardf; remquo; remquof; rint; rintf; round; roundf; roundl; scalbln; scalblnf; scalblnl; scalbn; scalbnl; scalbnf; ldexpf; ldexpl; __signbit; __signbitf; __signbitl; signgam; significand; significandf; sin; sinf; tan; tanf; tanh; tanhf; trunc; truncf; truncl; cabs; cabsf; drem; dremf; }; /* First added in 8.0-CURRENT */ FBSD_1.1 { carg; cargf; csqrt; csqrtf; logbl; nan; nanf; nanl; llrintl; lrintl; nearbyintl; rintl; exp2l; sinl; cosl; tanl; tgammaf; sqrtl; hypotl; cabsl; csqrtl; remquol; remainderl; fmodl; acosl; asinl; atan2l; atanl; cargl; cproj; cprojf; cprojl; }; /* First added in 9.0-CURRENT */ FBSD_1.2 { __isnanf; cbrtl; cexp; cexpf; log2; log2f; }; /* First added in 10.0-CURRENT */ FBSD_1.3 { feclearexcept; fegetexceptflag; fetestexcept; fegetround; fesetround; fesetenv; acoshl; asinhl; atanhl; cacos; cacosf; cacosh; cacoshf; casin; casinf; casinh; casinhf; catan; catanf; catanh; catanhf; csin; csinf; csinh; csinhf; ccos; ccosf; ccosh; ccoshf; coshl; ctan; ctanf; ctanh; ctanhf; erfcl; erfl; expl; expm1l; lgammal; log10l; log1pl; log2l; logl; powl; sinhl; tanhl; /* Implemented as weak aliases for imprecise versions */ tgammal; }; /* First added in 11.0-CURRENT */ FBSD_1.4 { lgammal_r; }; /* First added in 12.0-CURRENT */ FBSD_1.5 { cacoshl; cacosl; casinhl; casinl; catanl; catanhl; clog; clogf; clogl; cpow; cpowf; cpowl; sincos; sincosf; sincosl; }; /* First added in 14.0-CURRENT */ FBSD_1.7 { + cexpl; cospi; cospif; cospil; sinpi; sinpif; sinpil; tanpi; tanpif; tanpil; }; diff --git a/lib/msun/src/s_cexpf.c b/lib/msun/ld128/s_cexpl.c similarity index 67% copy from lib/msun/src/s_cexpf.c copy to lib/msun/ld128/s_cexpl.c index b815c99af89f..24f5e3792115 100644 --- a/lib/msun/src/s_cexpf.c +++ b/lib/msun/ld128/s_cexpl.c @@ -1,91 +1,94 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * - * Copyright (c) 2011 David Schultz + * Copyright (c) 2019 Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ #include __FBSDID("$FreeBSD$"); #include +#include #include +#include "fpmath.h" #include "math_private.h" +#include "k_expl.h" -static const uint32_t -exp_ovfl = 0x42b17218, /* MAX_EXP * ln2 ~= 88.722839355 */ -cexp_ovfl = 0x43400074; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ +/* XXX cexpl() should be converted to use bits likeo src/s_cexp.c. */ -float complex -cexpf(float complex z) -{ - float x, y, exp_x; - uint32_t hx, hy; +static const long double +cexp_ovfl = 2.27892930024498818830197576893019292e+04L, +exp_ovfl = 1.13565234062941439494919310779707649e+04L; - x = crealf(z); - y = cimagf(z); +long double complex +cexpl(long double complex z) +{ + long double c, exp_x, s, x, y; - GET_FLOAT_WORD(hy, y); - hy &= 0x7fffffff; + x = creall(z); + y = cimagl(z); /* cexp(x + I 0) = exp(x) + I 0 */ - if (hy == 0) - return (CMPLXF(expf(x), y)); - GET_FLOAT_WORD(hx, x); + if (y == 0) + return (CMPLXL(expl(x), y)); /* cexp(0 + I y) = cos(y) + I sin(y) */ - if ((hx & 0x7fffffff) == 0) - return (CMPLXF(cosf(y), sinf(y))); + if (x == 0) { + sincosl(y, &s, &c); + return (CMPLXL(c, s)); + } - if (hy >= 0x7f800000) { - if ((hx & 0x7fffffff) != 0x7f800000) { + if (!isfinite(y)) { + if (isfinite(x) || isnan(x)) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ - return (CMPLXF(y - y, y - y)); - } else if (hx & 0x80000000) { + return (CMPLXL(y - y, y - y)); + } else if (isinf(x) && copysignl(1.L, x) < 0) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ - return (CMPLXF(0.0, 0.0)); + return (CMPLXL(0.0, 0.0)); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ - return (CMPLXF(x, y - y)); + return (CMPLXL(x, y - y)); } } - if (hx >= exp_ovfl && hx <= cexp_ovfl) { + if (x > exp_ovfl && x < cexp_ovfl) { /* - * x is between 88.7 and 192, so we must scale to avoid - * overflow in expf(x). + * x is between exp_ovfl and cexp_ovfl, so we must scale to + * avoid overflow in exp(x). */ - return (__ldexp_cexpf(z, 0)); + return (__ldexp_cexpl(z, 0)); } else { /* * Cases covered here: * - x < exp_ovfl and exp(x) won't overflow (common case) * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0 * - x = +-Inf (generated by exp()) * - x = NaN (spurious inexact exception from y) */ - exp_x = expf(x); - return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y))); + exp_x = expl(x); + sincosl(y, &s, &c); + return (CMPLXL(exp_x * c, exp_x * s)); } } diff --git a/lib/msun/src/s_cexp.c b/lib/msun/ld80/s_cexpl.c similarity index 63% copy from lib/msun/src/s_cexp.c copy to lib/msun/ld80/s_cexpl.c index 2ef8ba1972ca..5cc35f6d2514 100644 --- a/lib/msun/src/s_cexp.c +++ b/lib/msun/ld80/s_cexpl.c @@ -1,91 +1,107 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2011 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, 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. + * + * src/s_cexp.c converted to long double complex by Steven G. Kargl */ #include __FBSDID("$FreeBSD$"); #include -#include +#include +#ifdef __i386__ +#include +#endif +#include "fpmath.h" +#include "math.h" #include "math_private.h" +#include "k_expl.h" -static const uint32_t -exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */ -cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ - -double complex -cexp(double complex z) +long double complex +cexpl (long double complex z) { - double x, y, exp_x; - uint32_t hx, hy, lx, ly; + long double c, exp_x, s, x, y; + uint64_t lx, ly; + uint16_t hx, hy; - x = creal(z); - y = cimag(z); + ENTERI(); - EXTRACT_WORDS(hy, ly, y); - hy &= 0x7fffffff; + x = creall(z); + y = cimagl(z); + + EXTRACT_LDBL80_WORDS(hy, ly, y); + hy &= 0x7fff; /* cexp(x + I 0) = exp(x) + I 0 */ if ((hy | ly) == 0) - return (CMPLX(exp(x), y)); - EXTRACT_WORDS(hx, lx, x); + RETURNI(CMPLXL(expl(x), y)); + EXTRACT_LDBL80_WORDS(hx, lx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ - if (((hx & 0x7fffffff) | lx) == 0) - return (CMPLX(cos(y), sin(y))); + if (((hx & 0x7fff) | lx) == 0) { + sincosl(y, &s, &c); + RETURNI(CMPLXL(c, s)); + } - if (hy >= 0x7ff00000) { - if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { + if (hy >= 0x7fff) { + if ((hx & 0x7fff) < 0x7fff || ((hx & 0x7fff) == 0x7fff && + (lx & 0x7fffffffffffffffULL) != 0)) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ - return (CMPLX(y - y, y - y)); - } else if (hx & 0x80000000) { + RETURNI(CMPLXL(y - y, y - y)); + } else if (hx & 0x8000) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ - return (CMPLX(0.0, 0.0)); + RETURNI(CMPLXL(0.0, 0.0)); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ - return (CMPLX(x, y - y)); + RETURNI(CMPLXL(x, y - y)); } } - if (hx >= exp_ovfl && hx <= cexp_ovfl) { + /* + * exp_ovfl = 11356.5234062941439497 + * cexp_ovfl = 22755.3287906024445633 + */ + if ((hx == 0x400c && lx > 0xb17217f7d1cf79acULL) || + (hx == 0x400d && lx < 0xb1c6a8573de9768cULL)) { /* - * x is between 709.7 and 1454.3, so we must scale to avoid - * overflow in exp(x). + * x is between exp_ovfl and cexp_ovfl, so we must scale to + * avoid overflow in exp(x). */ - return (__ldexp_cexp(z, 0)); + RETURNI(__ldexp_cexpl(z, 0)); } else { /* * Cases covered here: * - x < exp_ovfl and exp(x) won't overflow (common case) * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0 * - x = +-Inf (generated by exp()) * - x = NaN (spurious inexact exception from y) */ - exp_x = exp(x); - return (CMPLX(exp_x * cos(y), exp_x * sin(y))); + exp_x = expl(x); + sincosl(y, &s, &c); + RETURNI(CMPLXL(exp_x * c, exp_x * s)); } } diff --git a/lib/msun/man/cexp.3 b/lib/msun/man/cexp.3 index 776e6cee823e..27ab3e9c2098 100644 --- a/lib/msun/man/cexp.3 +++ b/lib/msun/man/cexp.3 @@ -1,113 +1,118 @@ .\" Copyright (c) 2011 David Schultz .\" All rights reserved. .\" .\" Redistribution and use in source and binary forms, with or without .\" modification, are permitted provided that the following conditions .\" are met: .\" 1. Redistributions of source code must retain the above copyright .\" notice, 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$ .\" -.Dd March 6, 2011 +.Dd November 3, 2021 .Dt CEXP 3 .Os .Sh NAME .Nm cexp , -.Nm cexpf +.Nm cexpf , +.Nm cexpl .Nd complex exponential functions .Sh LIBRARY .Lb libm .Sh SYNOPSIS .In complex.h .Ft double complex .Fn cexp "double complex z" .Ft float complex .Fn cexpf "float complex z" +.Ft long double complex +.Fn cexpl "long double complex z" .Sh DESCRIPTION The -.Fn cexp +.Fn cexp , +.Fn cexpf , and -.Fn cexpf +.Fn cexpl functions compute the complex exponential of .Fa z , also known as .Em cis Ns ( Ns .Fa z Ns ) . .Sh RETURN VALUES For real numbers .Fa x and .Fa y , .Fn cexp behaves according to Euler's formula: .Bd -ragged -offset indent .Fn cexp "x + I*y" = .Po Sy e Ns ** Ns .Fa x * .Em cos Ns Po Ns .Fa y Ns Pc Pc + Po Ns .Sy I * .Sy e Ns ** Ns .Fa x * .Em sin Ns Po Ns .Fa y Ns Pc Pc .Ed .Pp Generally speaking, infinities, zeroes and \*(Nas are handled as would be expected from this identity given the usual rules of floating-point arithmetic. However, care is taken to avoid generating \*(Nas when they are not deserved. For example, mathematically we expect that .Fo cimag .Fn cexp "x + I*0" Fc = 0 regardless of the value of .Fa x , and .Fn cexp preserves this identity even if .Fa x is \*(If or \*(Na. Likewise, .Fn cexp "-\*(If + I*y" = 0 and .Fo creal .Fn cexp "\*(If + I*y" Fc = \*(If for any .Fa y (even though the latter property is only mathematically true for representable .Fa y . ) If .Fa y is not finite, the sign of the result is indeterminate. .Sh SEE ALSO .Xr complex 3 , .Xr exp 3 , .Xr math 3 .Sh STANDARDS The -.Fn cexp +.Fn cexp , +.Fn cexpf , and -.Fn cexpf +.Fn cexpl functions conform to .St -isoC-99 . diff --git a/lib/msun/man/complex.3 b/lib/msun/man/complex.3 index f1acfbe6da74..8cc0b7f97c52 100644 --- a/lib/msun/man/complex.3 +++ b/lib/msun/man/complex.3 @@ -1,123 +1,129 @@ .\" Copyright (c) 2011 Murray Stokely .\" 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 ``AS IS'' AND .\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE .\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE .\" ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE .\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL .\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS .\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) .\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT .\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY .\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF .\" SUCH DAMAGE. .\" .\" $FreeBSD$ .\" -.Dd June 19, 2018 +.Dd November 3, 2021 .Dt COMPLEX 3 .Os .Sh NAME .Nm complex .Nd "complex arithmetic" .Sh LIBRARY .Lb libm .Sh SYNOPSIS .In complex.h .Sh DESCRIPTION These functions support complex arithmetic in the C math library. .Sh "LIST OF FUNCTIONS" Each of the following .Vt "double complex" functions has a .Vt "float complex" counterpart with an .Ql f appended to the name and a .Vt "long double complex" counterpart with an .Ql l appended. As an example, the .Vt "float complex" and .Vt "long double complex" counterparts of .Ft double .Fn cabs "double complex z" are .Ft float .Fn cabsf "float complex z" and .Ft "long double" .Fn cabsl "long double complex z" , respectively. .de Cl .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) csqrt complex square root .El .Ss Exponential Function .Cl cexp exponential base e .El .Ss Natural logarithm Function .Cl clog natural logarithm .El .\" Section 7.3.9 of ISO C99 standard .Ss Manipulation Functions .Cl carg compute the argument (i.e., phase angle) cimag compute the imaginary part conj compute the complex conjugate cproj compute projection onto Riemann sphere creal compute the real part .El .\" Section 7.3.5-6 of ISO C99 standard .Ss Trigonometric and Hyperbolic Functions .Cl cacos arc cosine cacosh arc hyperbolic cosine casin arc sine casinh arc hyperbolic sine catan arc tangent catanh arc hyperbolic tangent ccos cosine ccosh hyperbolic cosine cpow power function csin sine csinh hyperbolic sine ctan tangent ctanh hyperbolic tangent .El .Sh SEE ALSO .Xr fenv 3 , .Xr ieee 3 , .Xr math 3 , .Xr tgmath 3 .Rs .%T "ISO/IEC 9899:TC3" .%U http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1256.pdf .Re .Sh STANDARDS The .In complex.h functions described here conform to .St -isoC-99 . +.Sh BUGS +The power functions, +.Fn cpowf, cpow , +and +.Fn cpowl , +are implemented, but the code was neither reviewed nor tested. diff --git a/lib/msun/src/k_sincosl.h b/lib/msun/src/k_sincosl.h index 4d4dc69f7820..6425f14a1ea0 100644 --- a/lib/msun/src/k_sincosl.h +++ b/lib/msun/src/k_sincosl.h @@ -1,134 +1,137 @@ /*- * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== * * k_sinl.c and k_cosl.c merged by Steven G. Kargl */ #include __FBSDID("$FreeBSD$"); #if LDBL_MANT_DIG == 64 /* ld80 version of k_sincosl.c. */ #if defined(__amd64__) || defined(__i386__) /* Long double constants are slow on these arches, and broken on i386. */ static const volatile double C1hi = 0.041666666666666664, /* 0x15555555555555.0p-57 */ C1lo = 2.2598839032744733e-18, /* 0x14d80000000000.0p-111 */ S1hi = -0.16666666666666666, /* -0x15555555555555.0p-55 */ S1lo = -9.2563760475949941e-18; /* -0x15580000000000.0p-109 */ #define S1 ((long double)S1hi + S1lo) #define C1 ((long double)C1hi + C1lo) #else static const long double C1 = 0.0416666666666666666136L, /* 0xaaaaaaaaaaaaaa9b.0p-68 */ S1 = -0.166666666666666666671L; /* -0xaaaaaaaaaaaaaaab.0p-66 */ #endif static const double C2 = -0.0013888888888888874, /* -0x16c16c16c16c10.0p-62 */ C3 = 0.000024801587301571716, /* 0x1a01a01a018e22.0p-68 */ C4 = -0.00000027557319215507120, /* -0x127e4fb7602f22.0p-74 */ C5 = 0.0000000020876754400407278, /* 0x11eed8caaeccf1.0p-81 */ C6 = -1.1470297442401303e-11, /* -0x19393412bd1529.0p-89 */ C7 = 4.7383039476436467e-14, /* 0x1aac9d9af5c43e.0p-97 */ S2 = 0.0083333333333333332, /* 0x11111111111111.0p-59 */ S3 = -0.00019841269841269427, /* -0x1a01a01a019f81.0p-65 */ S4 = 0.0000027557319223597490, /* 0x171de3a55560f7.0p-71 */ S5 = -0.000000025052108218074604, /* -0x1ae64564f16cad.0p-78 */ S6 = 1.6059006598854211e-10, /* 0x161242b90243b5.0p-85 */ S7 = -7.6429779983024564e-13, /* -0x1ae42ebd1b2e00.0p-93 */ S8 = 2.6174587166648325e-15; /* 0x179372ea0b3f64.0p-101 */ static inline void __kernel_sincosl(long double x, long double y, int iy, long double *sn, long double *cs) { long double hz, r, v, w, z; z = x * x; v = z * x; /* * XXX Replace Horner scheme with an algorithm suitable for CPUs * with more complex pipelines. */ r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * S8))))); if (iy == 0) *sn = x + v * (S1 + z * r); else *sn = x - ((z * (y / 2 - v * r) - y) - v * S1); hz = z / 2; w = 1 - hz; r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 + z * C7)))))); *cs = w + (((1 - w) - hz) + (z * r - x * y)); } #elif LDBL_MANT_DIG == 113 /* ld128 version of k_sincosl.c. */ static const long double -C1 = 0.04166666666666666666666666666666658424671L, -C2 = -0.001388888888888888888888888888863490893732L, -C3 = 0.00002480158730158730158730158600795304914210L, -C4 = -0.2755731922398589065255474947078934284324e-6L, -C5 = 0.2087675698786809897659225313136400793948e-8L, -C6 = -0.1147074559772972315817149986812031204775e-10L, -C7 = 0.4779477332386808976875457937252120293400e-13L, S1 = -0.16666666666666666666666666666666666606732416116558L, S2 = 0.0083333333333333333333333333333331135404851288270047L, S3 = -0.00019841269841269841269841269839935785325638310428717L, S4 = 0.27557319223985890652557316053039946268333231205686e-5L, S5 = -0.25052108385441718775048214826384312253862930064745e-7L, S6 = 0.16059043836821614596571832194524392581082444805729e-9L, S7 = -0.76471637318198151807063387954939213287488216303768e-12L, S8 = 0.28114572543451292625024967174638477283187397621303e-14L; static const double -C8 = -0.1561920696721507929516718307820958119868e-15, -C9 = 0.4110317413744594971475941557607804508039e-18, -C10 = -0.8896592467191938803288521958313920156409e-21, -C11 = 0.1601061435794535138244346256065192782581e-23, S9 = -0.82206352458348947812512122163446202498005154296863e-17, S10 = 0.19572940011906109418080609928334380560135358385256e-19, S11 = -0.38680813379701966970673724299207480965452616911420e-22, S12 = 0.64038150078671872796678569586315881020659912139412e-25; +static const long double +C1 = 4.16666666666666666666666666666666667e-02L, +C2 = -1.38888888888888888888888888888888834e-03L, +C3 = 2.48015873015873015873015873015446795e-05L, +C4 = -2.75573192239858906525573190949988493e-07L, +C5 = 2.08767569878680989792098886701451072e-09L, +C6 = -1.14707455977297247136657111139971865e-11L, +C7 = 4.77947733238738518870113294139830239e-14L, +C8 = -1.56192069685858079920640872925306403e-16L, +C9 = 4.11031762320473354032038893429515732e-19L, +C10= -8.89679121027589608738005163931958096e-22L, +C11= 1.61171797801314301767074036661901531e-24L, +C12= -2.46748624357670948912574279501044295e-27L; + static inline void __kernel_sincosl(long double x, long double y, int iy, long double *sn, long double *cs) { long double hz, r, v, w, z; z = x * x; v = z * x; /* * XXX Replace Horner scheme with an algorithm suitable for CPUs * with more complex pipelines. */ r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 + z * (S9 + z * (S10 + z * (S11 + z * S12))))))))); if (iy == 0) *sn = x + v * (S1 + z * r); else - *cs = x - ((z * (y / 2 - v * r) - y) - v * S1); + *sn = x - ((z * (y / 2 - v * r) - y) - v * S1); hz = z / 2; w = 1 - hz; r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 + - z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11)))))))))); + z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * (C11+z*C12))))))))))); *cs = w + (((1 - w) - hz) + (z * r - x * y)); } #else #error "Unsupported long double format" #endif diff --git a/lib/msun/src/s_cexp.c b/lib/msun/src/s_cexp.c index 2ef8ba1972ca..a1f853eca4cc 100644 --- a/lib/msun/src/s_cexp.c +++ b/lib/msun/src/s_cexp.c @@ -1,91 +1,99 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2011 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ #include __FBSDID("$FreeBSD$"); #include +#include #include #include "math_private.h" static const uint32_t exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */ cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ double complex cexp(double complex z) { - double x, y, exp_x; + double c, exp_x, s, x, y; uint32_t hx, hy, lx, ly; x = creal(z); y = cimag(z); EXTRACT_WORDS(hy, ly, y); hy &= 0x7fffffff; /* cexp(x + I 0) = exp(x) + I 0 */ if ((hy | ly) == 0) 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 (CMPLX(cos(y), sin(y))); + if (((hx & 0x7fffffff) | lx) == 0) { + sincos(y, &s, &c); + return (CMPLX(c, s)); + } if (hy >= 0x7ff00000) { if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ return (CMPLX(y - y, y - y)); } else if (hx & 0x80000000) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ return (CMPLX(0.0, 0.0)); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ return (CMPLX(x, y - y)); } } if (hx >= exp_ovfl && hx <= cexp_ovfl) { /* * x is between 709.7 and 1454.3, so we must scale to avoid * overflow in exp(x). */ return (__ldexp_cexp(z, 0)); } else { /* * Cases covered here: * - x < exp_ovfl and exp(x) won't overflow (common case) * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0 * - x = +-Inf (generated by exp()) * - x = NaN (spurious inexact exception from y) */ exp_x = exp(x); - return (CMPLX(exp_x * cos(y), exp_x * sin(y))); + sincos(y, &s, &c); + return (CMPLX(exp_x * c, exp_x * s)); } } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(cexp, cexpl); +#endif diff --git a/lib/msun/src/s_cexpf.c b/lib/msun/src/s_cexpf.c index b815c99af89f..d905b74ff4bd 100644 --- a/lib/msun/src/s_cexpf.c +++ b/lib/msun/src/s_cexpf.c @@ -1,91 +1,94 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2011 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" static const uint32_t exp_ovfl = 0x42b17218, /* MAX_EXP * ln2 ~= 88.722839355 */ cexp_ovfl = 0x43400074; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ float complex cexpf(float complex z) { - float x, y, exp_x; + float c, exp_x, s, x, y; uint32_t hx, hy; x = crealf(z); y = cimagf(z); GET_FLOAT_WORD(hy, y); hy &= 0x7fffffff; /* cexp(x + I 0) = exp(x) + I 0 */ if (hy == 0) return (CMPLXF(expf(x), y)); GET_FLOAT_WORD(hx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ - if ((hx & 0x7fffffff) == 0) - return (CMPLXF(cosf(y), sinf(y))); + if ((hx & 0x7fffffff) == 0) { + sincosf(y, &s, &c); + return (CMPLXF(c, s)); + } if (hy >= 0x7f800000) { if ((hx & 0x7fffffff) != 0x7f800000) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ return (CMPLXF(y - y, y - y)); } else if (hx & 0x80000000) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ return (CMPLXF(0.0, 0.0)); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ return (CMPLXF(x, y - y)); } } if (hx >= exp_ovfl && hx <= cexp_ovfl) { /* * x is between 88.7 and 192, so we must scale to avoid * overflow in expf(x). */ return (__ldexp_cexpf(z, 0)); } else { /* * Cases covered here: * - x < exp_ovfl and exp(x) won't overflow (common case) * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0 * - x = +-Inf (generated by exp()) * - x = NaN (spurious inexact exception from y) */ exp_x = expf(x); - return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y))); + sincosf(y, &s, &c); + return (CMPLXF(exp_x * c, exp_x * s)); } } diff --git a/lib/msun/src/s_cosl.c b/lib/msun/src/s_cosl.c index 46a2e8620a22..3d066483f227 100644 --- a/lib/msun/src/s_cosl.c +++ b/lib/msun/src/s_cosl.c @@ -1,97 +1,102 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2007 Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice unmodified, this list of conditions, and the following * disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ #include __FBSDID("$FreeBSD$"); /* * Limited testing on pseudorandom numbers drawn within [-2e8:4e8] shows * an accuracy of <= 0.7412 ULP. */ #include #ifdef __i386__ #include #endif +#include "fpmath.h" #include "math.h" #include "math_private.h" #if LDBL_MANT_DIG == 64 #include "../ld80/e_rem_pio2l.h" +static const union IEEEl2bits +pio4u = LD80C(0xc90fdaa22168c235, -00001, 7.85398163397448309628e-01L); +#define pio4 (pio4u.e) #elif LDBL_MANT_DIG == 113 #include "../ld128/e_rem_pio2l.h" +long double pio4 = 7.85398163397448309615660845819875721e-1L; #else #error "Unsupported long double format" #endif long double cosl(long double x) { union IEEEl2bits z; int e0; long double y[2]; long double hi, lo; z.e = x; z.bits.sign = 0; /* If x = +-0 or x is a subnormal number, then cos(x) = 1 */ if (z.bits.exp == 0) return (1.0); /* If x = NaN or Inf, then cos(x) = NaN. */ if (z.bits.exp == 32767) return ((x - x) / (x - x)); ENTERI(); /* Optimize the case where x is already within range. */ - if (z.e < M_PI_4) + if (z.e < pio4) RETURNI(__kernel_cosl(z.e, 0)); e0 = __ieee754_rem_pio2l(x, y); hi = y[0]; lo = y[1]; switch (e0 & 3) { case 0: hi = __kernel_cosl(hi, lo); break; case 1: hi = - __kernel_sinl(hi, lo, 1); break; case 2: hi = - __kernel_cosl(hi, lo); break; case 3: hi = __kernel_sinl(hi, lo, 1); break; } RETURNI(hi); }