pow,powf(3),__ieee754_rem_pio2(f): Avoid negative integer left shift UB

A compiler clever enough to know that z is positive with a non-zero
biased exponent could, for example, optimize away the scal

pow,powf(3),__ieee754_rem_pio2(f): Avoid negative integer left shift UB

A compiler clever enough to know that z is positive with a non-zero
biased exponent could, for example, optimize away the scalbnf(z,n) in
pow() because behavior for left shift of negative values is undefined.
`n` is negative when y*log2(|x|) < -0.5. i.e. |x^y| < sqrt(0.5)

The intended behavior for operator<< in this code is to shift the two's
complement representation of the first operand.

In the pow() functions, the result is added to the IEEE 754 exponent of
z = 2^y'. n may be negative enough to underflow the biased IEEE 754
exponent below zero, which is manifested in the sign bit of j
(which would correspond to the IEEE 754 sign bit).

The conversion from uint32_t to int32_t for out-of-int32_t-range values
is implementation defined. The assumed behavior of interpreting the
uint32_t value as a two's complement representation of a signed value
is already assumed in many parts of the code, such as uses of
GET_FLOAT_WORD() with signed integers.

This code passes all the current tests, and makes some out of tree
fuzzing tests pass again rather than hit UB (detailed in the commentary
of the pull request).

Signed-off-by: Karl Tomlinson <karlt+@karlt.net>
Reviewed by: imp, steve kargl, dim
Pull Request: https://github.com/freebsd/freebsd-src/pull/1137

show more ...

lib/msun: Cleanup after $FreeBSD$ removal

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

PR: 276669
MFC after: 1 week

lib: Remove ancient SCCS tags.

Remove ancient SCCS tags from the tree, automated scripting, with two
minor fixup to keep things compiling. All the common forms in the tree
were removed with a perl s

lib: Remove ancient SCCS tags.

Remove ancient SCCS tags from the tree, automated scripting, with two
minor fixup to keep things compiling. All the common forms in the tree
were removed with a perl script.

Sponsored by: Netflix

show more ...

Remove $FreeBSD$: one-line .c pattern

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

Clean up libm use of the __ieee754_ prefix

This removes the __ieee754_ prefix from a number of the math functions.
msun/src/math_private.h contains the statement that

/*
* ieee style elementar

Clean up libm use of the __ieee754_ prefix

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 3f70824172feb82ea3dcdb3866b54fe0eb7cd890
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

show more ...

pow(3): Workaround possible signed shift Undefined Behavior.

j is int32_t and thus j<<31 is undefined if j==1.

Hinted by: muusl-lib (git 688d3da0f1730daddbc954bbc2d27cc96ceee04c)
Discussed with: fr

pow(3): Workaround possible signed shift Undefined Behavior.

j is int32_t and thus j<<31 is undefined if j==1.

Hinted by: muusl-lib (git 688d3da0f1730daddbc954bbc2d27cc96ceee04c)
Discussed with: freebsd-numerics (kargl)

show more ...

Add a macro nan_mix() and use it to get NaN results that are (bitwise)
independent of the precision in most cases. This is mainly to simplify
checking for errors. r176266 did this for e_pow[f].c us

Add a macro nan_mix() and use it to get NaN results that are (bitwise)
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.

show more ...

msun: add ld80/ld128 powl, cpow, cpowf, cpowl from openbsd

This corresponds to the latest status (hasn't changed in 9+
years) from openbsd of ld80/ld128 powl, and source cpowf, cpow,
cpowl (the comp

msun: add ld80/ld128 powl, cpow, cpowf, cpowl from openbsd

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

show more ...

Reduce diff between msun/src/e_pow.c and msun/src/e_powf.c.

Remove unnecessary casts, use integer literal constants instead of
floating point constants where possible, and introduce three const
stat

Reduce diff between msun/src/e_pow.c and msun/src/e_powf.c.

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

show more ...

lib/msun: remove trailing whitespace from e_pow.c

Submitted by: Steve Kargl
MFC after: 1 week

Revert r286421, now that the fix for LLVM PR24343 is in.

Apply an atrocious hack to libm's pow() implementation, to work around
LLVM PR 24343. This fixes the graphics/png build on i386, which
otherwise failed during its tests because of floating point exc

Apply an atrocious hack to libm's pow() implementation, to work around
LLVM PR 24343. This fixes the graphics/png build on i386, which
otherwise failed during its tests because of floating point exceptions,
leading to >5000 skipped other ports.

show more ...

IFC @r269962

Submitted by: Anish Gupta (akgupt3@gmail.com)

Merge head from 7/28

Sync to HEAD@r269943.

* Update the comments to agree with commit r226595.

* While here, fix a nearby typo in a comment.

IFC @230489 (pending review).

IFC @ r226824

Per IEEE754r, pow(1, y) is 1 even if y is NaN, and pow(-1, +-Inf) is 1.

MFC after: 2 weeks

Use the expression (x+0.0)-(y+0.0) instead of x+y when mixing NaN arg(s).
This uses 2 tricks to improve consistency so that more serious problems
aren't hidden in simple regression tests by noise for

Use the expression (x+0.0)-(y+0.0) instead of x+y when mixing NaN arg(s).
This uses 2 tricks to improve consistency so that more serious problems
aren't hidden in simple regression tests by noise for the NaNs:

- for a signaling NaN, adding 0.0 generates the invalid exception and
converts to a quiet NaN, and doesn't have too many effects for other
types of args (it converts -0 to +0 in some rounding modes, but that
hopefully doesn't change the result after adding the NaN arg). This
avoids some inconsistencies on i386 and ia64. On these arches, the
result of an operation on 2 NaNs is apparently the largest or the
smallest of the NaNs as bits (consistently largest or smallest for
each arch, but the opposite). I forget which way the comparison
goes and if the sign bit affects it. The quiet bit is is handled
poorly by not always setting it before the comparision or ignoring
it. Thus if one of the args was originally a signaling NaN and the
other was originally a quiet NaN, then the result depends too much
on whether the signaling NaN has been quieted at this point, which
in turn depends on optimizations and promotions. E.g., passing float
signaling NaNs to double functions must quiet them on conversion;
on i387, loading a signaling NaN of type float or double (but not
long double) into a register involves a conversion, so it quiets
signaling NaNs, so if the addition has 2 register operands than it
only sees quiet NaNs, but if the addition has a memory operand then
it sees a signaling NaN iff it is in the memory operand.

- subtraction instead of addition is used to avoid a dubious optimization
in old versions of gcc. For SSE operations, mixing of NaNs apparently
always gives the target operand. This is not as good as the i387
and ia64 behaviour. It doesn't mix NaNs at all, and makes addition
not quite commutative. Old versions of gcc sometimes rewrite x+y
to y+x and thus give different results (in bits) for NaNs. gcc-3.3.3
rewrites x+y to y+x for one of pow() and powf() but not the other,
so starting from float NaN args x and y, powf(x, y) was almost always
different from pow(x, y).

These tricks won't give consistency of 2-arg float and double functions
with long double ones on amd64, since long double ones use the i387
which has different semantics from SSE.

Convert to __FBSDID().

show more ...

Reduce diffs against vendor source (Sun fdlibm 5.3).

Fixed the sign of the result in some overflow and underflow cases (ones
where the exponent is an odd integer and the base is negative).

Obtained from: fdlibm-5.3

Sun finally released a new version

Fixed the sign of the result in some overflow and underflow cases (ones
where the exponent is an odd integer and the base is negative).

Obtained from: fdlibm-5.3

Sun finally released a new version of fdlibm just a coupe of weeks
ago. It only fixes 3 bugs (this one, another one in pow() that we
already have (rev.1.9), and one in tan(). I've learned too much about
powf() lately, so this fix was easy to merge. The patch is not verbatim,
because our base version has many differences for portability and I
didn't like global renaming of an unrelated variable to keep it separate
from the sign variable. This patch uses a new variable named sn for
the sign.

show more ...

e_pow.c:
Fixed pow(x, y) when x is very close to -1.0 and y is a very large odd
integer. E.g., pow(-1.0 - pow(2.0, -52.0), 1.0 + pow(2.0, 52.0)) was
0.0 instead of being very close to -exp(1.0).

PR

e_pow.c:
Fixed pow(x, y) when x is very close to -1.0 and y is a very large odd
integer. E.g., pow(-1.0 - pow(2.0, -52.0), 1.0 + pow(2.0, 52.0)) was
0.0 instead of being very close to -exp(1.0).

PR: 39236
Submitted by: Stephen L Moshier <steve@moshier.net>

e_powf.c:
Apply the same patch although it is just cosmetic because odd integers
large enough to cause the problem are too large to be precisely represented
as floats.

MFC after: 1 week

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