This is a fixed version of 888796ade284.
PR: 277783
Reported by: Victor Stinner
Reviewed by: emaste
MFC after: 1 week
(cherry picked from commit 888796ade2842486d3167067e8034254c38aadd3)
(cherry picked from commit e77ad954bb825983b4346b9cc646c9c910b1be24)
(cherry picked from commit 34f746cc7f8a8dd261027a8b392b76e70adc8438)
These do not use __FBSDID but instead use bare char arrays.
Reviewed by: imp, emaste
Differential Revision: https://reviews.freebsd.org/D41957
(cherry picked from commit eba230afba4932f02a1ca44efc797cf7499a5cb0)
This removes the __ieee754_ prefix from a number of the math functions.
msun/src/math_private.h contains the statement that
/*
* ieee style elementary functions
*
* We rename functions here to improve other sources' diffability
* against fdlibm.
*/
#define __ieee754_sqrt sqrt
...
Here, fdlibm refers to https://netlib.org/fdlibm. It is seen from
https://netlib.org/fdlibm/readme that this prefix was used to
differentiate between different standards:
Wrapper functions will twist the result of the ieee754
function to comply to the standard specified by the value
of _LIB_VERSION
if _LIB_VERSION = _IEEE_, return the ieee754 result;
if _LIB_VERSION = _SVID_, return SVID result;
if _LIB_VERSION = _XOPEN_, return XOPEN result;
if _LIB_VERSION = _POSIX_, return POSIX/ANSI result.
(These are macros, see fdlibm.h for their definition.)
AFAICT, FreeBSD has never supported these wrappers. In addition, as C99,
principally the long double, functions were added to libm, this
convention was not maintained. Given that only 148 of 324 files under
lib/msun contain a "Copyright (C) 1993 by Sun Microsystems" statement,
the removal of the __ieee754_ prefix provides consistency across all
source files.
The last time someone compared lib/msun to fdlibm appears to be
commit 3f70824172
Author: David Schultz <das@FreeBSD.org>
Date: Fri Feb 4 18:26:06 2005 +0000
Reduce diffs against vendor source (Sun fdlibm 5.3).
The most recent fdlibm RCS string that appears in a Sun Microsystem
copyrighted file is date "95/01/18". With Oracle Corporation's
acquisition of Sun Microsystems in 2009, it is unlikely that fdlibm will
ever be updated. A search for fdlibm at https://opensource.oracle.com/
yields no hits.
Finally, OpenBSD removed the use of this prefix over 21 years ago. pSee
revision 1.6 of OpenBSD's math_private.h.
Note: this does not drop the __ieee754_ prefix from the trigonometric
argument reduction functions, e.g., __ieee754_rem_pio2. These functions
are internal to the libm and exported through Symbol.map; and thus,
reserved for the implementation.
PR: 272783
MFC after: 1 week
patch to fix half-cycle trigonometric functions
Paul Zimmermann, a MPFR developer, contacted me about large errors in
the half-cycle trigonometric functions. I've have investigated these
issues and developed the attached patch. The float, double, and ld80
(long double) changes have been tested.
Caveat emptor: The ld128 changes have not been compiled. The ld128
changes have not been tested. I do not have access to a system that
uses ld128 floating point.
Here is an itemized list of changes:
* lib/msun/src/math_private.h:
. Add fast floor macros to compute the integer part of |x| for
0 <= |x| 01xp(N-1), where N is the precision of the type of x.
These macros are used in the half-cycle trigonometric functions
(e.g., sinpi(x)).
. The FFLOOR80 macros is used with the Intel 80-bit extended double
functions. This macors corrects an off-by-one error, which led to
enormous error for |x| > 0x1p32.
* lib/msun/src/s_cospif.c:
* lib/msun/src/s_cospi.c:
* lib/msun/ld80/s_cospil.c:
. Update Copyright years.
. Use FFLOOR*() macro to get integer part of |x|.
. Correct handle the range 0x1p(N-1) <= |x| < 0x1pN. Here, one needs
to determine if the integral value of |x| is even or odd to choose
+1 or -1. If |x| >= 0x1pN, always return +1.
* lib/msun/src/s_sinpif.c:
* lib/msun/src/s_sinpi.c:
* lib/msun/ld80/s_sinpil.c:
. Update Copyright years.
. Use FFLOOR*() macro to get integer part of |x|.
* lib/msun/src/s_tanpif.c:
* lib/msun/src/s_tanpi.c:
* lib/msun/ld80/s_tanpil.c:
. Update Copyright years.
. For +-0.5, return +-inf. Previously, tanpi[fl]() returned an NaN.
. Use FFLOOR*() to get integer part of |x|. Need to determine if the
integer part is even or odd. This is used to set +-0 for |x|
integral
and +-inf for (n+1/2).
. For 0x1p(N-1) <= |x| < 0x1pN need to determine if x is an even or
odd
integer to select +0 or -0. For |x| >= 0x1pN, it is always an even
integer, select 0.
. Note, tanpi[fl](x) is an odd function, so one needs to consider
tanpi[fl](-|x|) = - tanpi[fl](|x|).
* lib/msun/ld128/s_cospil.c:
* lib/msun/ld128/s_sinpil.c:
* lib/msun/ld128/s_tanpil.c:
. Update Copyright years.
. These routines use an FFLOOR128 macros, which likely should be
replaced by a bit twiddling algorithm.
. The same considerations above are applied to 0x1p112 <= |x| <
0x1p113,
and |x| >= 0x1p113 cases.
. Note, even and odd determination used fmodl(x,2.), which is likely
slow.
PR: 272742
MFC after: 1 week
David Das (das@) committed Bruce Evan's (bde's) WIP code for
expl() and logl() in git revision 25a4d6bfda. That code
included instrumentation that allowed bde to generate pari
scripts used in testing/debugging. This patch removes that
instrumentation as it is unlikely that others will ever use it.
* math/libm/msun/src/math_private.h:
. Remove bde's macros for the generation of pari scripts.
* math/libm/msun/ld128/s_expl.c:
* math/libm/msun/ld128/s_logl.c:
* math/libm/msun/ld80/s_expl.c:
* math/libm/msun/ld80/s_logl.c:
. Remove bde's DOPRINT_START macro.
. Change RETURNP to RETURNF.
. Change RETURN2P to RETURNF. Adjust arguments as needed.
. Change RETURNPI to RETURNI.
. Change RETURN2PI to RETURNI. Adjust arguments as needed.
PR: 272765
MFC after: 1 week
The current versions of lib/msun/src/s_cospi.c, s_sinpi.c and s_tanpi.c
all exhibit the same defect. After checking for various numeric ranges,
they check to see whether the input argument is a NaN or an Infinity.
However, the code uses a value of 0x7f80000 instead of the correct value
of 0x7ff00000.
If you review s_cospif.c, s_sinpif.c, and s_tanpif.c, you will see that
the equivalent statements in these functions are accurate and have
appropriate source comments.
The impact of these defects is to flag some valid input values as
invalid and raise a pole error (divide by zero).
Reported by: Paul Green <Paul.Green@stratus.com>
PR: 272539
MFC after: 1 week
The SPDX folks have obsoleted the BSD-2-Clause-FreeBSD identifier. Catch
up to that fact and revert to their recommended match of BSD-2-Clause.
Discussed with: pfg
MFC After: 3 days
Sponsored by: Netflix
- s/modfied/modified/
- s/minimun/minimum/
While here, fix some mandoc warnings:
- whitespace at end of input line
- unusual Xr punctuation
- missing comma before name
Obtained from: NetBSD
MFC after: 5 days
Damian McGuckin <damianm at esi dot com dot au> noted that the accuracy
claims in the code for cbrt(3) and cbrtl(3) were incorrect. Fix the
comments to more accurately describe the accuracies.
PR: 265603
MFC after: 3 days
Since https://github.com/llvm/llvm-project/commit/ca75ac5f04f2, clang 15
has a new warning about _Generic selection expressions, such as used in
math.h:
lib/libc/gdtoa/_ldtoa.c:82:10: error: due to lvalue conversion of the controlling expression, association of type 'volatile float' will never be selected because it is qualified [-Werror,-Wunreachable-code-generic-assoc]
switch (fpclassify(u.e)) {
^
lib/msun/src/math.h:109:2: note: expanded from macro 'fpclassify'
__fp_type_select(x, __fpclassifyf, __fpclassifyd, __fpclassifyl)
^
lib/msun/src/math.h:85:14: note: expanded from macro '__fp_type_select'
volatile float: f(x), \
^
This is because the controlling expression always undergoes lvalue
conversion first, dropping any cv-qualifiers. The 'const', 'volatile',
and 'volatile const' associations will therefore never be used.
MFC after: 1 week
Reviewed by: theraven
Differential Revision: https://reviews.freebsd.org/D35815
All supported compilers (modern versions of GCC and clang) support
this.
Many places didn't have an #else so would just silently do the wrong
thing. Ancient versions of icc (the original motivation for this) are
no longer a compiler FreeBSD supports.
PR: 263102 (exp-run)
Reviewed by: brooks, imp
Differential Revision: https://reviews.freebsd.org/D34797
. Disconnect imprecise.c from the build. This file can be deleted.
. Add b_tgammal.c to the build for ld80 and ld128 targets. The ld128
is a 'git mv' of imprecise.c to ld128/b_tgammal.c.
* lib/msun/ld80/b_expl.c:
. New file. Implement __exp__D for ld80 targets. This is based on
bsdsrc/b_exp.c.
* lib/msun/ld80/b_logl.c:
. New file. Implement __log__D for ld80 targets. This is based on
bsdsrc/b_log.c.
* lib/msun/ld80/b_tgammal.c b/lib/msun/ld80/b_tgammal.c
. New file. Implement tgammal(x) for ld80 targets.
Submitted by: Steve Kargl
Differential Revision: https://reviews.freebsd.org/D33444
Reviewed by: pfg
This look like a copy and paste leftover.
Reported by: enh@google.com (via freebsd-numerics@)
Reviewed by: Steve Kargl <sgk@troutmask.apl.washington.edu>
MFC after: 3 days
Some of the functions in msun can be implemented using a compiler
builtin function to generate a small number of instructions. Implement
this support in fma, fmax, fmin, and sqrt on arm64.
Care must be taken as the builtin can be implemented as a function
call on some architectures that lack direct support. In these cases
we need to use the original code path.
As we don't set errno on failure build with -fno-math-errno so the
toolchain doesn't convert a builtin into a function call when it
detects a failure, e.g. gcc will add a call to sqrt when the input
is negative leading to an infinite loop.
Sponsored by: The FreeBSD Foundation
Differential Revision: https://reviews.freebsd.org/D32801
In ac76bc1145, I added a few volatiles to work around ctrig_test
failures with {inf,inf}. This is not necessary anymore now, since in
3b00222f15 we added -fp-exception-behavior=maytrap for clang >= 10 in
libm's Makefile. (The flag tells clang to use stricter floating point
semantics, which libm depends on.)
PR: 244732, 254911
Fixes: ac76bc1145
MFC after: 3 days
The change implements cexpl() for both ld80 and ld128 architectures.
Testing was done on x86_64 and aarch64 systems.
Along the way sincos[fl]() use an optimization that reduces the argument
to being done one rather than twice. This optimization actually pointed
to a bug in the ld128 version of sincosl(), which is now fixed. In
addition, the minmax polynomial coefficients for sincosl() have been
updated.
A concise log of the file-by-file changes follows.
* include/complex.h:
. Add a prototype for cexpl().
* lib/msun/Makefile:
. Add s_cexpl.c to the build.
. Setup a link for cexpl.3 to cexp.3.
* lib/msun/Symbol.map:
. Expose cexpl symbol in libm shared library.
* lib/msun/ld128/s_cexpl.c:
* Implementation of cexpl() for 128-bit long double architectures.
Tested on an aarch64 system.
* lib/msun/ld80/s_cexpl.c:
* Implementation of cexpl() for Intel 80-bit long double.
* lib/msun/man/cexp.3:
. Document cexpl().
* lib/msun/man/complex.3:
. Add a BUGS section about cpow[fl].
* lib/msun/src/s_cexp.c:
. Include float.h for weak references on 53-bit long double targets.
. Use sincos() to reduce argument reduction cost.
* lib/msun/src/s_cexpf.c:
. Use sincosf() to reduce argument reduction cost.
* lib/msun/src/k_sincosl.h:
. Catch up with the new minmax polynomial coefficients for the kernel for
the 128-bit cosl() implementation.
. BUG FIX: *cs was used where *sn should have been. This means that sinl()
was no computed correctly when iy != 0.
* lib/msun/src/s_cosl.c:
. Include fpmath.h to get access to IEEEl2bits.
. Replace M_PI_4 with pio4, a 64-bit or 113-bit approximation for pi / 4.
PR: 216862
MFC after: 1 week
The patch fixes the omission of '#include <float.h>', which is needed for
the weak reference on systems with LDBL_MANT_DIG == DBL_MANT_DIG.
PR: 218514
MFC after: 1 week
Both IEEE-754 2008 and ISO/IEC TS 18661-4 define the half-cycle
trignometric functions cospi, sinpi, and tanpi. The attached
patch implements cospi[fl], sinpi[fl], and tanpi[fl]. Limited
testing on the cospi and sinpi reveal a max ULP less than 0.89;
while tanpi is more problematic with a max ULP less than 2.01
in the interval [0,0.5]. The algorithms used in these functions
are documented in {ks}_cospi.c, {ks}_sinpi.c, and s_tanpi.c.
Note. I no longer have access to a system with ld128 and
adequate support to compile and test the ld128 implementations
of these functions. Given the almost complete lack of input from
others on improvements to libm, I doubt that anyone cares. If
someone does care, the ld128 files contain a number of FIXME comments,
and in particular, while the polynomial coefficients are given
I did not update the polynomial algorithms to properly use the
coefficients.
PR: 218514
MFC after: 2 weeks
These files were copied from MUSL. Add the standard copyright notice and
SPDX-License-Identifier: MIT consistent with our new draft license
policy. It reads word for word the same as the MIT license on the SPDX
web site. Add a pointer to the MUSL COPYIRGHT file which contains a list
of all authors of MUSL.
Sponsored by: Netflix
Noticed by: Steve Kargl
Summary:
From Steve Kargl:
Paul Zimmermann has identified a bug in Openlibm's powf(),
which is identical to FreeBSD's libm. Both derived from
fdlibm. https://github.com/JuliaMath/openlibm/issues/212.
Consider
% cat h.c
int
main(void)
{
float x, y, z;
x = 0x1.ffffecp-1F;
y = -0x1.000002p+27F;
z = 0x1.557a86p115F;
printf("%e %e %e <-- should be %e\n", x, y, powf(x,y), z);
return 0;
}
% cc -o h -fno-builtin h.c -lm && ./h
9.999994e-01 -1.342177e+08 inf <-- should be 5.540807e+34
Reviewers: manu
Subscribers: imp, andrew, emaste
Differential Revision: https://reviews.freebsd.org/D31865
After df3b437c1e, older gcc's such as
4.2.1 (still used on earlier branches for e.g. mips and powerpc) and
6.3.0 (still used for some cross-builds) started throwing bogus errors
like:
In file included from /workspace/src/lib/msun/src/s_llround.c:11:0:
/workspace/src/lib/msun/src/s_lround.c:54:31: error: initializer element is not constant
static const type dtype_min = type_min - 0.5;
^~~~~~~~
/workspace/src/lib/msun/src/s_lround.c:55:31: error: initializer element is not constant
static const type dtype_max = type_max + 0.5;
^~~~~~~~
Since 'type_min' and 'type_max' are constants declared just above these
lines this error is nonsensical, but older gcc's are not smart enough.
Work around the error by reusing the (type)DTYPE_MIN and (type)DTYPE_MAX
macros, so I can MFC this right away, unbreaking a few stable builds.
MFC after: immediately
It turned out that the (type)DTYPE_MAX conversions at the top of
s_lround.c are now emitted as cvtsi2sd instructions, at least on SSE
capable CPUs. This caused the FE_INEXACT flag to always be set, at least
for the double and float variants. Under clang 11, the whole INRANGE()
comparisons were still optimized away, but this has "improved" in clang
12, due to stricter adherence to the -ffp-exception-behavior=maytrap
compiler flag.
To avoid run-time integer to float conversions, use static constants
instead, so they are computed at compile time, and the INRANGE()
statements are optimized away again, if applicable.
While here, use an integer instead of a floating type to store the test
results in lround_test.c, as this is more appropriate, and we can also
drop the volatile hack.
Reported by: arichardson
MFC after: 3 days
When compiling parts of math.h with clang using a C standard before C11,
and using -pedantic, it will result in warnings similar to:
bug254714.c:5:11: warning: '_Generic' is a C11 extension [-Wc11-extensions]
return !isfinite(1.0);
^
/usr/include/math.h:111:21: note: expanded from macro 'isfinite'
^
/usr/include/math.h:82:39: note: expanded from macro '__fp_type_select'
^
This is because the block that enables use of _Generic is conditional
not only on C11, but also on whether the compiler advertises support for
C generic selections via __has_extension(c_generic_selections).
To work around the warning without having to pessimize the code, use the
__extension__ keyword, which is supported by both clang and gcc. While
here, remove the check for __clang__, as _Generic has been supported for
a long time by gcc too now.
Reported by: yuri
PR: 254714
MFC after: 1 week
This caused LDBL_MANT_DIG to not be defined and therefore the scalbnl
alias was not being emitted for double==long double platforms.
Fixes: 760b2ffc ("Update scalbn* functions to the musl versions")
Reported by: Jenkins
The only diff compared to musl is a minor change to scalbnl() to replace
musl's union ldshape with union IEEEl2bits.
This fixes the scalbn tests on non-x86 (since x86 has an assembly version
that is used instead).
Musl commit messages:
commit 8c44a060243f04283ca68dad199aab90336141db
Author: Szabolcs Nagy <nsz@port70.net>
Date: Mon Apr 3 02:38:13 2017 +0200
fix scalbn when result is in the subnormal range
in nearest rounding mode scalbn could introduce double rounding error
when an intermediate value and the final result were both in the
subnormal range e.g.
scalbn(0x1.7ffffffffffffp-1, -1073)
returned 0x1p-1073 instead of 0x1p-1074, because the intermediate
computation got rounded to 0x1.8p-1023.
with the fix an intermediate value can only be in the subnormal range
if the final result is 0 which is correct even after double rounding.
(there still can be two roundings so signals may be raised twice, but
that's only observable with trapping exceptions which is not supported.)
commit 2eaed464e2080d8321d3903b71086a1ecfc4ee4a
Author: Szabolcs Nagy <nsz@port70.net>
Date: Wed Sep 4 15:52:54 2013 +0000
math: use float_t and double_t in scalbnf and scalbn
remove STRICT_ASSIGN (c99 semantics is assumed) and use the conventional
union to prepare the scaling factor (so libm.h is no longer needed)
commit 1b77b9072f374bd26eb0574b83a0d5f18d75ec60
Author: Szabolcs Nagy <nsz@port70.net>
Date: Thu Aug 15 10:07:46 2013 +0000
math: minor scalbn*.c simplification
commit c4359e01303da2755fe7e8033826b132eb3659b1
Author: Szabolcs Nagy <nsz@port70.net>
Date: Tue Nov 13 10:55:35 2012 +0100
math: excess precision fix modf, modff, scalbn, scalbnf
old code was correct only if the result was stored (without the
excess precision) or musl was compiled with -ffloat-store.
now we use STRICT_ASSIGN to work around the issue.
(see note 160 in c11 section 6.8.6.4)
commit 666271c105e4137bdfa195e217799d74143370d4
Author: Szabolcs Nagy <nsz@port70.net>
Date: Tue Nov 13 10:30:40 2012 +0100
math: fix scalbn and scalbnf on overflow/underflow
old code was correct only if the result was stored (without the
excess precision) or musl was compiled with -ffloat-store.
(see note 160 in n1570.pdf section 6.8.6.4)
commit 8051e08e10d2b739fcfcbc6bc7466e8d77fa49f1
Author: nsz <nsz@port70.net>
Date: Mon Mar 19 10:54:07 2012 +0100
simplify scalbn*.c implementations
The old scalbn.c was wrong and slow, the new one is just slow.
(scalbn(0x1p+1023,-2097) should give 0x1p-1074, but the old code gave 0)
Reviewed By: dim
Differential Revision: https://reviews.freebsd.org/D28872
This applies musl commit b02eed9c4841913d690a2d0029737d72615384fe by
Szabolcs Nagy and updates the tests accordingly. This also allows
removing an XFAIL from the test.
musl commit message:
complex: fix ctanh(+-0+i*nan) and ctanh(+-0+-i*inf)
These cases were incorrect in C11 as described by
http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1886.htm
PR: 217528
Reviewed By: dim
MFC after: 1 week
Differential Revision: https://reviews.freebsd.org/D28578
This adjusts the factor used to scale the subnormal numbers, so it
becomes the right value after adjusting its exponent. Thanks to Steve
Kargl for finding the most elegant fix.
Also enable the hypot tests, and add a test case for this bug.
PR: 253313
MFC after: 1 week
This sprinkles a few strategic volatiles in an attempt to defeat clang's
optimization interfering with the expected floating-point exception
flags.
Reported by: lwhsu
PR: 244732
MFC after: 3 days
The intel compiler support has badly decayed over the years. Stop
pretending that we support it. Note, I've stopped short of requiring
gcc builtin support with this commit since other compilers may be used
to build non-base software and we need to support those so more
investigation is needed before simplifying further.
by Steve Kargl:
- Use sincos[f] instead of a call to cos[f] and a call to sin[f].
- While here, alphabetize declaration.
Submitted by: sgk at troutmask.apl.washington.edu (Steve Kargl)
The "for" loop on big endian was inverting all the bits instead of
just the words
Issue reported by TestSuite (msun lib nan_test case)
Submitted by: Renato Riolino <renato.riolino@eldorado.org.br>
Submitted by: Fernando Valle <fernando.valle@eldorado.org.br>
Reviewed by: pfg, alfredo
Approved by: jhibbits (mentor)
Sponsored by: Eldorado Research Institute (eldorado.org.br)
Differential Revision: https://reviews.freebsd.org/D23926
negative numbers (invoking undefined behavior)
Summary:
Various paths through hypot(x, y) will multiply x and y by a power of
two, perform the calculation in a range where IEEE-754 provides greater
precision, then undo the multiplication to determine the true result.
Undoing that multiplication is implemented as t1*w, where t1=2**k.
2**k is often computed by taking the high word of 1.0, then adding k<<20
(for doubles or long doubles) or k<<23 (for floats) to it, then
overwriting that high word. But when k is negative this left-shifts a
negative value -- and that's undefined behavior in many editions of C
and C++.
This patch should fix all hypot implementations to compute 2**k without
triggering this particular bit of undefined behavior.
Test Plan: I've only very lightly tested out the hypot(double, double)
change, in SpiderMonkey's JavaScript engine, for consistency with prior
behavior. The other functions' changes have more or less only been
eyeballed. Careful examination appreciated! Do note, however, that an
error in any of these changes would most likely produce a value that is
incorrect by a factor of two, so any mistake would most likely be
glaring if invoked.
Submitted by: Jeff Walden <jwalden@mit.edu>
Obtained from: https://github.com/freebsd/freebsd/pull/414
Reviewed by: dim, lwhsu
MFC after: 1 week
Differential Revision: https://reviews.freebsd.org/D22354
Alter bsd.compat.mk to set MACHINE and MACHINE_ARCH when included
directly so MD paths in Makefiles work. In the process centralize
setting them in LIBCOMPATWMAKEENV.
Alter .PATH and CFLAGS settings in work when the Makefile is included.
While here only support LIB32 on supported platforms rather than always
enabling it and requiring users of MK_LIB32 to filter based
TARGET/MACHINE_ARCH.
The net effect of this change is to make Makefile.libcompat only build
compatability libraries.
Changes relative to r354449:
Correct detection of the compiler type when bsd.compat.mk is used
outside Makefile.libcompat. Previously it always matched the clang
case.
Set LDFLAGS including the linker emulation for mips where -m32 seems to
be insufficent.
Reviewed by: imp, kib (origional version in r354449)
Obtained from: CheriBSD (conceptually)
Sponsored by: DARPA, AFRL
Differential Revision: https://reviews.freebsd.org/D22251
