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
Even though clang comes with a number of internal CUDA wrapper headers,
compiling sample CUDA programs will result in errors similar to:
In file included from <built-in>:1:
In file included from /usr/lib/clang/9.0.0/include/__clang_cuda_runtime_wrapper.h:204:
/usr/home/arr/cuda/var/cuda-repo-10-0-local-10.0.130-410.48/usr/local/cuda-10.0//include/crt/math_functions.hpp:2910:7: error: no matching function for call to '__isnan'
if (__isnan(a)) {
^~~~~~~
/usr/lib/clang/9.0.0/include/__clang_cuda_device_functions.h:460:16: note: candidate function not viable: call to __device__ function from __host__ function
__DEVICE__ int __isnan(double __a) { return __nv_isnand(__a); }
^
CUDA expects __isnan() and __isnanf() declarations to be available,
which are glibc specific extensions, equivalent to the regular isnan()
and isnanf().
To provide these, define __isnan() and __isnanf() as aliases of the
already existing static inline functions __inline_isnan() and
__inline_isnanf() from math.h.
Reported by: arrowd
PR: 241550
MFC after: 1 week
Clang from trunk recently added a warning for when implicit int-to-float
conversions cause a loss of precision. The code in question is designed
to be able to handle that, so add explicit casts to silence this.
Submitted by: James Clarke <jrtc27@jrtc27.com>
Reviewed by: dim
Obtained from: CheriBSD
MFC after: 1 week
Sponsored by: DARPA, AFRL
Differential Revision: https://reviews.freebsd.org/D21913
C/C++) in exp(3), expf(3), expm1(3) and expm1f(3) during intermediate
computations that compute the IEEE-754 bit pattern for |2**k| for
integer |k|.
The implementations of exp(3), expf(3), expm1(3) and expm1f(3) need to
compute IEEE-754 bit patterns for 2**k in certain places. (k is an
integer and 2**k is exactly representable in IEEE-754.)
Currently they do things like 0x3FF0'0000+(k<<20), which is to say they
take the bit pattern representing 1 and then add directly to the
exponent field to get the desired power of two. This is fine when k is
non-negative.
But when k<0 (and certain classes of input trigger this), this
left-shifts a negative number -- an operation with undefined behavior in
C and C++.
The desired semantics can be achieved by instead adding the
possibly-negative k to the IEEE-754 exponent bias to get the desired
exponent field, _then_ shifting that into its proper overall position.
(Note that in case of s_expm1.c and s_expm1f.c, there are SET_HIGH_WORD
and SET_FLOAT_WORD uses further down in each of these files that perform
shift operations involving k, but by these points k's range has been
restricted to 2 < k <= 56, and the shift operations under those
circumstances can't do anything that would be UB.)
Submitted by: Jeff Walden, https://github.com/jswalden
Obtained from: https://github.com/freebsd/freebsd/pull/411
Obtained from: https://github.com/freebsd/freebsd/pull/412
MFC after: 3 days
Replace calls to sinf(x) and cosf(x) with a single call to sincosf().
Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>
Reviewed by: bde
Approved by: grog
MFC after: 3 days
j is int32_t and thus j<<31 is undefined if j==1.
Hinted by: muusl-lib (git 688d3da0f1730daddbc954bbc2d27cc96ceee04c)
Discussed with: freebsd-numerics (kargl)
The long double aliases of double functions are only exposed as aliases if
LDBL_MANT_DIG is 53 (same as DBL_MANT_DIG). Without float.h included these
files were not exposing weak aliases as expected, leading to link failures
if programs use the *l functions. This should fix editors/calligra on
targets with 64-bit long double, which uses erfl and erfcl. Found on
powerpc64.
Reviewed by: kargl@
of NaNs before possible returning a NaN.
The remquo*() and remainder*() functions should now give bitwise identical
results across arches and implementations, and bitwise consistent results
(with lower precisions having truncated mantissas) across precisions. x86
already had consistency across amd64 and i386 and precisions by using the
i387 consistently and normally not using the C versions. Inconsistencies
for C reqmquol() were first detected on sparc64.
Remove double second clearing of the sign bit and extra blank lines.
remainder*(x, y) and remquo*(x, y, quo) were broken for y = 0 by changing
multiplication by y to addition of y. (When y is 0, the result should be
NaN but became 1 for finite x.)
Use a new macro nan_mix_op() to give more control over the mixing, and
expand comments.
Recent re-testing missed finding this bug since I only tested the macro
version on amd64 and i386 and these arches don't use the C versions (they
use either asm versions or builtins).
Reported by: enh via freebsd-numerics
This is a follow-up to r336299.
* lib/msun/Makefile:
. Remove polevll.c
* lib/msun/ld80/e_powl.c:
. Copy contents of polevll.c to here. This is the only consumer of
these functions. Make functions 'static inline'.
. Make reducl a 'static inline' function.
* lib/msun/man/exp.3:
. Remove BUGS section that no longer applies.
* lib/msun/src/math_private.h:
. Remove prototypes of __p1evll() and __polevll()
* lib/msun/src/s_cpow.c:
* lib/msun/src/s_cpowf.c:
* lib/msun/src/s_cpowl.c
. Include math_private.h.
. Use the CMPLX macro from either C99 or math_private.h (depends on
compiler support) instead of the problematic use of complex I.
Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>
PR: 229876
MFC after: 1 week
This was open-coded in range reduction for trig and exp functions. Now
there are 3 static inline functions rnint[fl]() that replace open-coded
expressions, and type-generic irint() and i64rint() macros that hide the
complications for efficiently using non-generic irint() and irintl()
functions and casts.
Special details:
ld128/e_rem_pio2l.h needs to use i64rint() since it needs a 46-bit integer
result. Everything else only needs a (less than) 32-bit integer result so
uses irint().
Float and double cases now use float_t and double_t locally instead of
STRICT_ASSIGN() to avoid bugs in extra precision.
On amd64, inline asm is now only used for irint() on long doubles. The SSE
asm for irint() on amd64 only existed because the ifdef tangles made the
correct method of simply casting to int for this case non-obvious.
cc1: warnings being treated as errors
/usr/src/lib/msun/src/s_cpow.c: In function 'cpow':
/usr/src/lib/msun/src/s_cpow.c:63: warning: implicit declaration of function 'CMPLX'
This is a follow-up to r336299.
* lib/msun/Makefile:
. Remove polevll.c
* lib/msun/ld80/e_powl.c:
. Copy contents of polevll.c to here. This is the only consumer of
these functions. Make functions 'static inline'.
. Make reducl a 'static inline' function.
* lib/msun/man/exp.3:
. Remove BUGS section that no longer applies.
* lib/msun/src/math_private.h:
. Remove prototypes of __p1evll() and __polevll()
* lib/msun/src/s_cpow.c:
* lib/msun/src/s_cpowf.c:
* lib/msun/src/s_cpowl.c
. Use the CMPLX macro from either C99 or math_private.h (depends of
compiler support) instead of the problematic use of complex I.
Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>
PR: 229876
MFC after: 1 week
with 1 huge component and 1 tiny (but nowhere near denormal) component.
Rescale earlier so that a scale factor of 2 can be combined with a non-
scale divisor of 2, so that the division doesn't shift out a bit. In the
usual case where the scale factor is just 1, the division may shift out a
bit, but then the underflow is not spurious and the inaccuracies are harder
to fix.
Remove the STDC CX_LIMITED_RANGE pragma and its verbose comment. We still
don't have any C99 compilers (that support fenv pragmas), and if we did
then there are thousands of other places in libm that would need to use
them more than here.
The other cleanups are smaller.
and csqrtl().
When one component is huge and the other is tiny, scaling down the tiny
component gave spurious underflow.
When both components are denormal, not scaling them up gave inaccuracies
of 34+ ulps on not very carefully selected args. Fixing this reduces the
maximum error to 1.6 ulps on the same set of args (mosly not denormal ones).
The scaling used multiplication of a complex variable by 2, but clang messes
this on amd64 up by losing the sign of -0.0. Calculate the components
separately, as is well known to be needed for operations on more exceptional
values.
independent of the precision in most cases. This is mainly to simplify
checking for errors. r176266 did this for e_pow[f].c using a less
refined expression that often didn't work. r176276 fixes an error in
the log message for r176266. The main refinement is to always expand
to long double precision. See old log messages (especially these 2)
and the comment on the macro for more general details.
Specific details:
- using nan_mix() consistently for the new and old pow*() functions was
the only thing needed to make my consistency test for powl() vs pow()
pass on amd64.
- catrig[fl].c already had all the refinements, but open-coded.
- e_atan2[fl].c, e_fmod[fl].c and s_remquo[fl] only had primitive NaN
mixing.
- e_hypot[fl].c already had a different refined version of r176266. Refine
this further. nan_mix() is not directly usable here since we want to
clear the sign bit.
- e_remainder[f].c already had an earlier version of r176266.
- s_ccosh[f].c,/s_csinh[f].c already had a version equivalent to r176266.
Refine this further. nan_mix() is not directly usable here since the
expression has to handle some non-NaN cases.
- s_csqrt.[fl]: the mixing was special and mostly wrong. Partially fix the
special version.
- s_ctanh[f].c already had a version of r176266.
This corresponds to the latest status (hasn't changed in 9+
years) from openbsd of ld80/ld128 powl, and source cpowf, cpow,
cpowl (the complex power functions for float complex, double
complex, and long double complex) which are required for C99
compliance and were missing from FreeBSD. Also required for
some numerical codes using complex numbered Hamiltonians.
Thanks to jhb for tracking down the issue with making
weak_reference compile on powerpc.
When asked to review, bde said "I don't like it" - but
provided no actionable feedback or superior implementations.
Discussed with: jhb
Submitted by: jmd
Differential Revision: https://reviews.freebsd.org/D15919
Remove unnecessary casts, use integer literal constants instead of
floating point constants where possible, and introduce three const
static variables to hold 0.5, 0.25, and 1/3.
PR: 229420
Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>
MFC after: 1 week
As a component of atan2(y, x), the case of x == 1.0 is farmed out to
atan(y). The current implementation of this comparison is vulnerable
to signed integer underflow (that is, undefined behavior), and it's
performed in a somewhat more complicated way than it need be. Change
it to not be quite so cute, rather directly comparing the high/low
bits of x to the specific IEEE-754 bit pattern that encodes 1.0.
Note that while there are three different e_atan* files in the
relevant directory, only this one needs fixing. e_atan2f.c already
compares against the full bit pattern encoding 1.0f, while
e_atan2l.cuses bitwise-ands/ors/nots and so doesn't require a change.
Closes#130
Submitted by: Jeff Walden (@jswalden github PR #130)
Reviewed by: bde
MFC After: 1 month
Mainly focus on files that use BSD 2-Clause license, however the tool I
was using mis-identified many licenses so this was mostly a manual - error
prone - task.
The Software Package Data Exchange (SPDX) group provides a specification
to make it easier for automated tools to detect and summarize well known
opensource licenses. We are gradually adopting the specification, noting
that the tags are considered only advisory and do not, in any way,
superceed or replace the license texts.
LDBL_MAX is broken on i386:
https://lists.freebsd.org/pipermail/freebsd-numerics/2012-September/000288.html
Gcc has produced +Infinity for LDBL_MAX on i386 and amd64 with -m32
for some time, and newer versions of gcc are now warning that the
"floating constant exceeds range of 'long double'". Avoid this by
referring to proxy values instead.
Reviewed by: bde
Approved by: markj (mentor)
Sponsored by: Dell EMC Isilon
LDBL_MAX is broken on i386:
https://lists.freebsd.org/pipermail/freebsd-numerics/2012-September/000288.html
Gcc has produced +Infinity for LDBL_MAX on i386 and amd64 with -m32
for some time, and newer versions of gcc are now warning that the
"floating constant exceeds range of 'long double'". Avoid this by
referring to half the value of LDBL_MAX instead.
Reviewed by: bde
Approved by: markj (mentor)
Sponsored by: Dell EMC Isilon
The primary benefit of these functions is that argument
reduction is done once instead of twice in independent
calls to sin() and cos().
* lib/msun/Makefile:
. Add s_sincos[fl].c to the build.
. Add sincos.3 documentation.
. Add appropriate MLINKS.
* lib/msun/Symbol.map:
. Expose sincos[fl] symbols in dynamic libm.so.
* lib/msun/man/sincos.3:
. Documentation for sincos[fl].
* lib/msun/src/k_sincos.h:
. Kernel for sincos() function. This merges the individual kernels
for sin() and cos(). The merger offered an opportunity to re-arrange
the individual kernels for better performance.
* lib/msun/src/k_sincosf.h:
. Kernel for sincosf() function. This merges the individual kernels
for sinf() and cosf(). The merger offered an opportunity to re-arrange
the individual kernels for better performance.
* lib/msun/src/k_sincosl.h:
. Kernel for sincosl() function. This merges the individual kernels
for sinl() and cosl(). The merger offered an opportunity to re-arrange
the individual kernels for better performance.
* lib/msun/src/math.h:
. Add prototytpes for sincos[fl]().
* lib/msun/src/math_private.h:
. Add RETURNV macros. This is needed to reset fpsetprec on I386
hardware for a function with type void.
* lib/msun/src/s_sincos.c:
. Implementation of sincos() where sin() and cos() were merged into
one routine and possibly re-arranged for better performance.
* lib/msun/src/s_sincosf.c:
. Implementation of sincosf() where sinf() and cosf() were merged into
one routine and possibly re-arranged for better performance.
* lib/msun/src/s_sincosl.c:
. Implementation of sincosl() where sinl() and cosl() were merged into
one routine and possibly re-arranged for better performance.
PR: 215977, 218300
Submitted by: Steven G. Kargl <sgk@troutmask.apl.washington.edu>
MFC after: 1 month
Differential Revision: https://reviews.freebsd.org/D10765
