lib/msun: Cleanup after $FreeBSD$ removal

Remove no longer needed explicit inclusion of sys/cdefs.h.

PR: 276669
MFC after: 1 week

Remove $FreeBSD$: one-line .c pattern

Remove /^[\s*]*__FBSDID\("\$FreeBSD\$"\);?\s*\n/

Merge ^/head r352587 through r352763.

Do not left-shift a negative number (inducing undefined behavior in
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|

Do not left-shift a negative number (inducing undefined behavior in
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

show more ...

Merge from head

IFC @ r254014

IFC @253398

Merge fresh head.

Fix some regressions caused by the switch from gcc to clang. The fixes
are workarounds for various symptoms of the problem described in clang
bugs 3929, 8100, 8241, 10409, and 12958.

The regression

Fix some regressions caused by the switch from gcc to clang. The fixes
are workarounds for various symptoms of the problem described in clang
bugs 3929, 8100, 8241, 10409, and 12958.

The regression tests did their job: they failed, someone brought it
up on the mailing lists, and then the issue got ignored for 6 months.
Oops. There may still be some regressions for functions we don't have
test coverage for yet.

show more ...

IFC @230489 (pending review).

IFC @ r226824

Use STRICT_ASSIGN() to ensure that the compiler doesn't screw things
up by storing x in a wider type than it's supposed to.

Submitted by: bde

Remove some unnecessary initializations.

Obtained from: DragonFlyBSD

Fix some rather obscene code that has ambiguous if...if...else...
constructs in it.

As usual, use a minimax polynomial that is specialized for float
precision. The new polynomial has degree 4 instead of 10, and a maximum
error of 2**-30.04 ulps instead of 2**-33.15. This doesn't a

As usual, use a minimax polynomial that is specialized for float
precision. The new polynomial has degree 4 instead of 10, and a maximum
error of 2**-30.04 ulps instead of 2**-33.15. This doesn't affect the
final error significantly; the maximum error was and is about 0.5015
ulps on i386 -O1, and the number of cases with an error of > 0.5 ulps
is increased from 13851 to 14407.

Note that the error is only this close to 0.5 ulps due to excessive
extra precision caused by compiler bugs on i386. The extra precision
could be obtained intentionally, and is useful for keeping the error
of the hyperbolic float functions below 1 ulp, since these functions
are implemented using expm1f. My recent change for scaling by 2**k
had the unintentional side effect of retaining extra precision for
longer, so callers of expm1f see errors of more like 0.0015 ulps than
0.5015 ulps, and for the hyperbolic functions this reduces the maximum
error from nearly about 2 ulps to about 0.75 ulps.

This is about 10% faster on i386 (A64). expm1* is still very slow,
but now the float version is actually significantly faster. The
algorithm is very sophisticated but not very good except on machines
with fast division.

show more ...

Use a better method of scaling by 2**k. Instead of adding to the
exponent bits of the reduced result, construct 2**k (hopefully in
parallel with the construction of the reduced result) and multiply

Use a better method of scaling by 2**k. Instead of adding to the
exponent bits of the reduced result, construct 2**k (hopefully in
parallel with the construction of the reduced result) and multiply by
it. This tends to be much faster if the construction of 2**k is
actually in parallel, and might be faster even with no parallelism
since adjustment of the exponent requires a read-modify-wrtite at an
unfortunate time for pipelines.

In some cases involving exp2* on amd64 (A64), this change saves about
40 cycles or 30%. I think it is inherently only about 12 cycles faster
in these cases and the rest of the speedup is from partly-accidentally
avoiding compiler pessimizations (the construction of 2**k is now
manually scheduled for good results, and -O2 doesn't always mess this
up). In most cases on amd64 (A64) and i386 (A64) the speedup is about
20 cycles. The worst case that I found is expf on ia64 where this
change is a pessimization of about 10 cycles or 5%. The manual
scheduling for plain exp[f] is harder and not as tuned.

Details specific to expm1*:
- the saving is closer to 12 cycles than to 40 for expm1* on i386 (A64).
For some reason it is much larger for negative args.
- also convert to __FBSDID().

show more ...

Fix formatting, this is hard to explain, so I'll show one example.

- float ynf(int n, float x) /* wrapper ynf */
+float
+ynf(int n, float x) /* wrapper ynf */

This is because the __S

Fix formatting, this is hard to explain, so I'll show one example.

- float ynf(int n, float x) /* wrapper ynf */
+float
+ynf(int n, float x) /* wrapper ynf */

This is because the __STDC__ stuff was indented.

Reviewed by: md5

show more ...

Assume __STDC__, remove non-__STDC__ code.

Reviewed by: md5

$Id$ -> $FreeBSD$

Revert $FreeBSD$ to $Id$

Make the long-awaited change from $Id$ to $FreeBSD$

This will make a number of things easier in the future, as well as (finally!)
avoiding the Id-smashing problem which has plagued developers for so

Make the long-awaited change from $Id$ to $FreeBSD$

This will make a number of things easier in the future, as well as (finally!)
avoiding the Id-smashing problem which has plagued developers for so long.

Boy, I'm glad we're not using sup anymore. This update would have been
insane otherwise.

show more ...

Remove trailing whitespace.

J.T. Conklin's latest version of the Sun math library.

-- Begin comments from J.T. Conklin:
The most significant improvement is the addition of "float" versions
of the math functions that take float

J.T. Conklin's latest version of the Sun math library.

-- Begin comments from J.T. Conklin:
The most significant improvement is the addition of "float" versions
of the math functions that take float arguments, return floats, and do
all operations in floating point. This doesn't help (performance)
much on the i386, but they are still nice to have.

The float versions were orginally done by Cygnus' Ian Taylor when
fdlibm was integrated into the libm we support for embedded systems.
I gave Ian a copy of my libm as a starting point since I had already
fixed a lot of bugs & problems in Sun's original code. After he was
done, I cleaned it up a bit and integrated the changes back into my
libm.
-- End comments

Reviewed by: jkh
Submitted by: jtc

show more ...

Fix some rather obscene code that has ambiguous if...if...else...
constructs in it.

As usual, use a minimax polynomial that is specialized for float
precision. The new polynomial has degree 4 instead of 10, and a maximum
error of 2**-30.04 ulps instead of 2**-33.15. This doesn't a

As usual, use a minimax polynomial that is specialized for float
precision. The new polynomial has degree 4 instead of 10, and a maximum
error of 2**-30.04 ulps instead of 2**-33.15. This doesn't affect the
final error significantly; the maximum error was and is about 0.5015
ulps on i386 -O1, and the number of cases with an error of > 0.5 ulps
is increased from 13851 to 14407.

Note that the error is only this close to 0.5 ulps due to excessive
extra precision caused by compiler bugs on i386. The extra precision
could be obtained intentionally, and is useful for keeping the error
of the hyperbolic float functions below 1 ulp, since these functions
are implemented using expm1f. My recent change for scaling by 2**k
had the unintentional side effect of retaining extra precision for
longer, so callers of expm1f see errors of more like 0.0015 ulps than
0.5015 ulps, and for the hyperbolic functions this reduces the maximum
error from nearly about 2 ulps to about 0.75 ulps.

This is about 10% faster on i386 (A64). expm1* is still very slow,
but now the float version is actually significantly faster. The
algorithm is very sophisticated but not very good except on machines
with fast division.

show more ...