fnstsw does not wait for pending unmasked x87 floating-point exceptions
and it is the same as fstsw when all exceptions are masked which is the
only environment libc supports.

Some early x86_64 cpus (released before 2006) did not support sahf/lahf
instructions so they should be avoided (intel manual says they are only
supported if CPUID.80000001H:ECX.LAHF-SAHF[bit 0] = 1).

The workaround simplifies exp2l and expm1l because fucomip can be
used instead of the fucomp;fnstsw;sahf sequence copied from i386.

In fmodl and remainderl sahf is replaced by a simple bit test.

the idiomatic rounding of x is

  n = x + toint - toint;

where toint is either 1/EPSILON (x is non-negative) or 1.5/EPSILON
(x may be negative and nearest rounding mode is assumed) and EPSILON is
according to the evaluation precision (the type of toint is not very
important, because single precision float can represent the 1/EPSILON of
ieee binary128).

in case of FLT_EVAL_METHOD!=0 this avoids a useless store to double or
float precision, and the long double code became cleaner with
1/LDBL_EPSILON instead of ifdefs for toint.

__rem_pio2f and __rem_pio2 functions slightly changed semantics:
on i386 a double-rounding is avoided so close to half-way cases may
get evaluated differently eg. as sin(pi/4-eps) instead of cos(pi/4+eps)

The old code used the rounding idiom incorrectly:

  y = (double)(x + 0x1p52) - 0x1p52;

the cast is useless if FLT_EVAL_METHOD==0 and causes a second rounding
if FLT_EVAL_METHOD==2 which can give incorrect result in nearest rounding
mode, so the correct idiom is to add/sub a power-of-2 according to the
characteristics of double_t.

This did not cause actual bug because only i386 is affected where rint
is implemented in asm.

Other rounding functions use a similar idiom, but they give correct
results because they only rely on getting a neighboring integer result
and the rounding direction is fixed up separately independently of the
current rounding mode. However they should be fixed to use the idiom
correctly too.

previously the external definitions of these functions were omitted on
archs where long double is the same as double, since the code paths in
the math.h macros which would call them are unreachable. however, even
if they are unreachable, the definitions are still mandatory. omitting
them is invalid C, and in the case of a non-optimizing compiler, will
result in a link error.

This was not caught earlier because gcc incorrectly generates quiet
relational operators that never raise exceptions.

the previous commit was a no op in exp10l because LDBL_* macros
were implicitly 0 (the preprocessor does not warn about undefined
symbols).

__polevll, __p1evll and exp10l were provided on archs when long double
is the same as double. The first two were completely unused and exp10l
can be a wrapper around exp10.

modfl and sincosl were passing long double* instead of double*
to the wrapped double precision functions (on archs where long
double and double have the same size).
This is fixed now by using temporaries (this is not optimized
to a single branch so the generated code is a bit bigger).
Found by Morten Welinder.