r355120 | dim | 2019-11-26 23:01:09 +0100 (Tue, 26 Nov 2019) | 32 lines

The fdlibm hypot() implementations shouldn't potentially left-shift
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
Differential Revision: https://reviews.freebsd.org/D22354

Copy head (r256279) to stable/10 as part of the 10.0-RELEASE cycle.

Approved by: re (implicit)
Sponsored by: The FreeBSD Foundation

Remove some unnecessary initializations.

Obtained from: DragonFlyBSD

Alias hypotl() and cabsl() for platforms where long double is the same
as double.

Use fabs[f]() instead of bit fiddling for setting absolute values.
This makes little difference in float precision, but in double
precision gives a speedup of about 30% on amd64 (A64 CPU) and i386
(A64). This depends on fabs[f]() being inline and efficient. The
bit fiddling (or any use of SET_HIGH_WORD(), which libm does too
much because it was best on old 32-bit machines) always causes
packing overheads and sometimes causes stalls in the packing, since
it operates on only part of a variable in the double precision case.
It apparently did cause stalls in a critical path here.

Use the expression fabs(x+0.0)-fabs(y+0.0) instead of
fabs(x+0.0)+fabs(y+0.0) when mixing NaNs. This improves
consistency of the result by making it harder for the compiler to reorder
the operands. (FP addition is not necessarily commutative because the
order of operands makes a difference on some machines iff the operands are
both NaNs.)

Use the expression fabs(x+0.0)+fabs(y+0.0) instad of a+b (where a is
|x| or |y| and b is |y| or |x|) when mixing NaN arg(s).

hypot*() had its own foot shooting for mixing NaNs -- it swaps the
args so that |x| in bits is largest, but does this before quieting
signaling NaNs, so on amd64 (where the result of adding NaNs depends
on the order) it gets inconsistent results if setting the quiet bit
makes a difference, just like a similar ia64 and i387 hardware comparison.
The usual fix (see e_powf.c 1.13 for more details) of mixing using
(a+0.0)+-(b+0.0) doesn't work on amd64 if the args are swapped (since
the rder makes a difference with SSE). Fortunately, the original args
are unchanged and don't need to be swapped when we let the hardware
decide the mixing after quieting them, but we need to take their
absolute value.

hypotf() doesn't seem to have any real bugs masked by this non-bug.
On amd64, its maximum error in 2^32 trials on amd64 is now 0.8422 ulps,
and on i386 the maximum error is unchanged and about the same, except
with certain CFLAGS it magically drops to 0.5 (perfect rounding).

Convert to __FBSDID().

Reduce diffs against vendor source (Sun fdlibm 5.3).

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

Assume __STDC__, remove non-__STDC__ code.

Reviewed by: md5

$Id$ -> $FreeBSD$

Use __ieee754_sqrt() instead of sqrt() internally. Similarly for the
float versions. Using sqrt() was inefficient.

Obtained from: NetBSD

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 long.

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

Remove trailing whitespace.

This commit was generated by cvs2svn to compensate for changes in r2116,
which included commits to RCS files with non-trunk default branches.

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