the freebsd fma code failed to raise underflow exception in some
cases in nearest rounding mode (affects fmal too) e.g.

  fma(-0x1p-1000, 0x1.000001p-74, 0x1p-1022)

and the inexact exception may be raised spuriously since the fenv
is not saved/restored around the exact multiplication algorithm
(affects x86 fma too).

another issue is that the underflow behaviour when the rounded result
is the minimal normal number is target dependent, ieee754 allows two
ways to raise underflow for inexact results: raise if the result before
rounding is in the subnormal range (e.g. aarch64, arm, powerpc) or if
the result after rounding with infinite exponent range is in the
subnormal range (e.g. x86, mips, sh).

to avoid all these issues the algorithm was rewritten with mostly int
arithmetics and float arithmetics is only used to get correct rounding
and raise exceptions according to the behaviour of the target without
any fenv.h dependency. it also unifies x86 and non-x86 fma.

fmaf is not affected, fmal need to be fixed too.

this algorithm depends on a_clz_64 and it required a few spurious
instructions to make sure underflow exception is raised in a particular
corner case. (normally FORCE_EVAL(tiny*tiny) would be used for this,
but on i386 gcc is broken if the expression is constant
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57245
and there is no easy portable fix for the macro.)

while the official elfv2 abi for "powerpc64le" sets power8 as the
baseline isa, we use it for both little and big endian powerpc64
targets and need to maintain compatibility with pre-power8 models. the
instructions for sqrt, fabs, and fma are in the baseline isa; support
for the rest is conditional via predefined isa-level macros.

patch by David Edelsohn.

these were introduced in z196 and not available in the baseline (z900)
ISA level. use __HTM__ as an alternate indicator for ISA level, since
gcc did not define __ARCH__ until 7.x.

patch by David Edelsohn.

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

this should increase performance and reduce code size on aarch64.

the compiled code was checked against using __builtin_* instead
of inline asm with gcc-6.2.0.

lrint is two instructions.

c with inline asm is used because it is safer than a pure asm
implementation, this prevents ll{rint,round} to be an alias
of l{rint,round} (because the types don't match) and depends
on gcc style inline asm support.

ceil, floor, round, trunc can either raise inexact on finite
non-integer inputs or not raise any exceptions. the new
implementation does not raise exceptions while the generic
c code does.

on aarch64, the underflow exception is signaled before rounding
(ieee 754 allows both before and after rounding, but it must be
consistent), the generic fma c code signals it after rounding
so using single instruction fixes a slight conformance issue too.

partly following freebsd rev 279491
https://svnweb.freebsd.org/base?view=revision&revision=279491
(musl had some of the fixes before freebsd).

the change should not matter much for j0f, y0f, but it improves
j1f and y1f in [2.5,~3.75] (that is [0x40200000,~0x40700000]).
near roots (e.g. around 3.8317 for j1f) there are still large
ulp errors.

dropped code that tried to raise inexact.

j is int32_t and thus j<<31 is undefined if j==1, so j is changed to
uint32_t locally as a quick fix, the generated code is not affected.

(this is a strict conformance fix, future c standard may allow 1<<31,
see DR 463.  the bug was inherited from freebsd fdlibm, the proper fix
is to use uint32_t for all bit hacks, but that requires more intrusive
changes.)

reported by Daniel Sabogal

there was a copy paste error that could cause large ulp errors
in atan2l, atanl, asinl and acosl on aarch64, mips64 and mipsn32.

(the implementation is from freebsd fdlibm, but the tail end
of the polynomial was wrong. 128 bit long double functions
are not yet tested so this went undetected.)

expf(-NAN) was treated as expf(-large) which unconditionally
returns +0, so special case +-NAN.
reported by Petr Hosek.

commit e4355bd6bec89688e8c739cd7b4c76e675643dca moved the math asm
from external source files to inline asm, but unfortunately, all
current releases of clang use the wrong inline asm constraint codes
for float and double ("w" and "P" instead of "t" and "w",
respectively). this patch adds detection for the bug in configure,
and, for now, just disables the affected asm on broken clang versions.