as a fairly faithful implementation of the algorithm found in
PTP Tang, "Table-driven implementation of the Expm1 function
in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 18,
211-222 (1992).
Over the last 18-24 months, the code has under gone significant
optimization and testing.
Reviewed by: bde
Obtained from: bde (most of the optimizations)
* Use integral numerical constants, and let the compiler do the
conversion to long double.
ld128/s_expl.c:
* Use integral numerical constants, and let the compiler do the
conversion to long double.
* Use the ENTERI/RETURNI macros, which are no-ops on ld128. This
however makes the ld80 and ld128 identical.
Reviewed by: bde (as part of larger diff)
* In the special case x = -Inf or -NaN, use a micro-optimization
to eliminate the need to access u.xbits.man.
* Fix an off-by-one for small arguments |x| < 0x1p-65.
ld128/s_expl.c:
* In the special case x = -Inf or -NaN, use a micro-optimization
to eliminate the need to access u.xbits.manh and u.xbits.manl.
* Fix an off-by-one for small arguments |x| < 0x1p-114.
Obtained from: bde
* Update the evaluation of the polynomial. This allows the removal
of the now unused variables t23 and t45.
ld128/s_expl.c:
* Update the evaluation of the polynomial and the intermediate
result t. This update allows several numerical constants to be
written as double rather than long double constants. Update
the constants as appropriate.
Obtained from: bde
and use macros to access the e component of the unions. This allows
the portions of the code in ld80 to be identical to the ld128 code.
Obtained from: bde
The names now coincide with the name used in PTP Tang's paper.
* Rename the variable from s to tbl to better reflect that
this is a table, and to be consistent with the naming scheme
in s_exp2l.c
Reviewed by: bde (as part of larger diff)
* Update Copyright years to include 2013.
ld128/s_expl.c:
* Correct and update Copyright years. This code originated from
the ld80 version, so it should reflect the same time period.
Reviewed by: bde (as part of larger diff)
* Use ENTERI/RETURNI to allow the use of FP_PE on i386 target.
Reviewed by: das (and bde a long time ago)
Approved by: das (mentor)
Obtained from: bde (polynomial coefficients)
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.
libc.a and libc_p.a. In addition, define isnan in libm.a and libm_p.a,
but not in libm.so.
This makes it possible to statically link executables using both isnan
and isnanf with libc and libm.
Tested by: kargl
MFC after: 1 week
table and the requirement on trailing zero bits.
* Remove the __aligned() compiler directives as these were found
to have a negative effect on the produced code.
Submitted by: bde
Approved by: das (mentor)
. Change the API for the LD80C by removing the explicit passing
of the sign bit. The sign can be determined from the last
parameter of the macro.
. On i386, load long double by bit manipulations to work around
at least a gcc compiler issue. On non-i386 ld80 architectures,
use a simple assignment.
* ld80/s_expl.c:
. Update the only consumer of LD80C.
Submitted by: bde
Approved by: das (mentor)
. Fix the threshold for expl(x) where |x| is small.
. Also update the previously incorrect comment to match the
new threshold.
* ld128/s_expl.c:
. Re-order logic in exceptional cases to match the logic used in
other long double functions.
. Fix the threshold for expl(x) where is |x| is small.
. Also update the previously incorrect comment to match the
new threshold.
Submitted by: bde
Approved by: das (mentor)
. Guard a comment from reformatting by indent(1).
. Re-order variables in declarations to alphabetical order.
. Remove a banal comment.
* ld128/s_expl.c:
. Add a comment to point to ld80/s_expl.c for implementation details.
. Move the #define of INTERVAL to reduce the diff with ld80/s_expl.c.
. twom10000 does not need to be volatile, so move its declaration.
. Re-order variables in declarations to alphabetical order.
. Add a comment that describes the argument reduction.
. Remove the same banal comment found in ld80/s_expl.c.
Reviewed by: bde
Approved by: das (mentor)
Also, update the comment to describe the choice of using
a high and low decomposition of 2^(i/INTERNVAL) for
0 <= i <= INTERVAL in preparation for an implementation of
expm1l.
* Move the #define of INTERVAL above the comment, because the
comment refers to INTERVAL.
Reviewed by: bde
Approved by: das (mentor)
only affects i386. The double case was intentionally left broken
as an optimization, but we are getting closer to supporting
applications and/or kernels that change the (FreeBSD i386) default
rounding precision from FP_PD to FP_PE and never change it back,
and this requires the STRICT_ALIGN()s that were added to support
FP_PE to actually work in all precisions.
* Remove an extraneous semicolon at the end of a macro that was
supposed to be function-like.
Submitted by: bde
Approved by: das (mentor)
they need to refer to static constants, which C99 does not allow for
extern inline functions.
While here, change a comment in e_rem_pio2f.c to mention the correct
number of bits.
Reviewed by: bde
MFC after: 1 week
compatibility with the INTERVALS macro used in the soon-to-be-commmitted
expm1l() and someday-to-be-committed log*l() functions.
Add a comment into ld128/s_expl.c noting at gcc issue that was
deleted when rewriting ld80/e_expl.c as ld128/s_expl.c.
Requested by: bde
Approved by: das (mentor)
. Remove a few #ifdefs that should have been removed in the initial
commit.
. Sort fpmath.h to its rightful place.
* ld128/s_expl.c:
. Replace EXPMASK with its actual value.
. Sort fpmath.h to its rightful place.
Requested by: bde
Approved by: das (mentor)
format. These implementations are based on
PTP Tang, "Table-driven implementation of the exponential function
in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15,
144-157 (1989).
PR: standards/152415
Submitted by: kargl
Reviewed by: bde, das
Approved by: das (mentor)
quotation. Also make sure we have the same amount of columns in each row as
the number of columns we specify in the head arguments.
Reviewed by: brueffer
correct sign when the remainder was 0.
Fix a separate bug in remquo alone, in which the remainder and
quotient were both off by a bit in certain cases involving subnormal
remainders.
The bugs affected all platforms except amd64 and i386, on which the
routines are implemented in assembly.
PR: 166463
Submitted by: Ilya Burylov
MFC after: 2 weeks
the function bodies require only 2 to 10 instructions. However, it
leads to application binaries that refer to a private ABI, namely, the
softfloat innards in libc. This could complicate future changes in
the implementation of the floating-point emulation layer, so it seems
best to have programs refer to the official fe* entry points in libm.
use softfloat.
Thanks to Ian Lepore for testing and debugging this patch. The fenv
regression tests pass (at least for Ian's arm chip) with this change.
- Handle cases where exp(x) would overflow, but ccosh(x) ~= exp(x) / 2
shouldn't.
- Use the ccosh(x) ~= exp(x) / 2 approximation to simplify the calculation
when x is large.
Similarly for csinh(). Also fixed the return value of csinh(-Inf +- 0i).
exp(x) scaled down by some factor, and the challenge is doing this
accurately when exp(x) would overflow. This change replaces all of
the tricks we've been using with common __ldexp_exp() and
__ldexp_cexp() routines that handle all the scaling.
bde plans to improve on this further by moving the guts of exp() into
k_exp.c and handling the scaling in a more direct manner. But the
current approach is simple and adequate for now.
library," since complex.h, tgmath.h, and fenv.h are also part of the
math library. Replace the outdated sentence with some references to
the other parts.
- Rename __kernel_log() to k_log1p().
- Move some of the work that was previously done in the kernel log into
the callers. This enables further refactoring to improve accuracy or
speed, although I don't recall the details.
- Use extra precision when adding the final scaling term, which improves
accuracy.
- Describe and work around compiler problems that break some of the
multiprecision calculations.
A fix for a small bug is also included:
- Add a special case for log*(1). This is needed to ensure that log*(1) == +0
instead of -0, even when the rounding mode is FE_DOWNWARD.
Submitted by: bde
no longer "fast" on sparc64. (It really wasn't to begin with, since
the old implementation was using long doubles, and long doubles are
emulated in software on sparc64.)
round-to-nearest mode when the result, rounded to twice machine
precision, was exactly halfway between two machine-precision
values. The essence of the fix is to simulate a "sticky bit" in
the pathological cases, which is how hardware implementations
break the ties.
MFC after: 1 month
fenv.h that are currently inlined.
The definitions are provided in fenv.c via 'extern inline'
declaractions. This assumes the compiler handles 'extern inline' as
specified in C99, which has been true under FreeBSD since 8.0.
The goal is to eventually remove the 'static' keyword from the inline
definitions in fenv.h, so that non-inlined references all wind up
pointing to the same external definition like they're supposed to.
I am deferring the second step to provide a window where
newly-compiled apps will still link against old math libraries.
(This isn't supported, but there's no need to cause undue breakage.)
Reviewed by: stefanf, bde
on i386-class hardware for sinl and cosl. The hand-rolled argument
reduction have been replaced by e_rem_pio2l() implementations. To
preserve history the following commands have been executed:
svn cp src/e_rem_pio2.c ld80/e_rem_pio2l.h
mv ${HOME}/bde/ld80/e_rem_pio2l.c ld80/e_rem_pio2l.h
svn cp src/e_rem_pio2.c ld128/e_rem_pio2l.h
mv ${HOME}/bde/ld128/e_rem_pio2l.c ld128/e_rem_pio2l.h
The ld80 version has been tested by bde, das, and kargl over the
last few years (bde, das) and few months (kargl). An older ld128
version was tested by das. The committed version has only been
compiled tested via 'make universe'.
Approved by: das (mentor)
Obtained from: bde
with r219571 and re-enable building of cbrtl.
Implement the long double version for the cube root function, cbrtl.
The algorithm uses Newton's iterations with a crude estimate of the
cube root to converge to a result.
Reviewed by: bde
Approved by: das
implementing accurate logarithms in different bases. This is based
on an approach bde coded up years ago.
This function should always be inlined; it will be used in only a few
places, and rudimentary tests show a 40% performance improvement in
implementations of log2() and log10() on amd64.
The kernel takes a reduced argument x and returns the same polynomial
approximation as e_log.c, but omitting the low-order term. The low-order
term is much larger than the rest of the approximation, so the caller of
the kernel function can scale it to the appropriate base in extra precision
and obtain a much more accurate answer than by using log(x)/log(b).
Explanation by Steve:
jn[f](n,x) for certain ranges of x uses downward recursion to compute
the value of the function. The recursion sequence that is generated is
proportional to the actual desired value, so a normalization step is
taken. This normalization is j0[f](x) divided by the zeroth sequence
member. As Bruce notes, near the zeros of j0[f](x) the computed value
can have giga-ULP inaccuracy. I found for the 1st zero of j0f(x) only
the leading decimal digit is correct. The solution to the issue is
fairly straight forward. The zeros of j0(x) and j1(x) never coincide,
so as j0(x) approaches a zero, the normalization constant switches to
j1[f](x) divided by the 2nd sequence member. The expectation is that
j1[f](x) is a more accurately computed value.
PR: bin/144306
Submitted by: Steven G. Kargl <kargl@troutmask.apl.washington.edu>
Reviewed by: bde
MFC after: 7 days
and one under lib/msun/amd64. This avoids adding the identifiers to the
.text section, and moves them to the .comment section instead.
Suggested by: bde
Approved by: rpaulo (mentor)
macro expand to __isnanf() instead of isnanf() for float arguments.
This change is needed because isnanf() isn't declared in strict POSIX
or C99 mode.
Compatibility note: Apps using isnan(float) that are compiled after
this change won't link against an older libm.
Reported by: Florian Forster <octo@verplant.org>
bottom of the manpages and order them consistently.
GNU groff doesn't care about the ordering, and doesn't even mention
CAVEATS and SECURITY CONSIDERATIONS as common sections and where to put
them.
Found by: mdocml lint run
Reviewed by: ru
argument for fnstsw. Explicitely specify sizes for the XMM control and
status word and X87 control and status words.
Reviewed by: das
Tested by: avg
MFC after: 2 weeks
The amd64-specific bits of msun use an undocumented constraint, which is
less likely to be supported by other compilers (such as Clang). Change
the code to use a more common machine constraint.
Obtained from: /projects/clangbsd/
FPA floating-point format is identical to the VFP format,
but is always stored in big-endian.
Introduce _IEEE_WORD_ORDER to describe the byte-order of
the FP representation.
Obtained from: Juniper Networks, Inc
#pragma STDC CX_LIMITED_RANGE ON
the "ON" needs to be in caps. gcc doesn't understand this pragma
anyway and assumes it is always on in any case, but icc supports
it and cares about the case.
conj() instead of using expressions like z * I. The latter is bad for
several reasons:
1. It is implemented using arithmetic, which is unnecessary, and can
generate floating point exceptions, contrary to the requirements on
these functions.
2. gcc implements complex multiplication using a formula that breaks
down for infinities, e.g., it gives INFINITY * I == nan + inf I.
- When y/x is huge, it's faster and more accurate to return pi/2
instead of pi - pi/2.
- There's no need for 3 lines of bit fiddling to compute -z.
- Fix a comment.
at compile time regardless of the dynamic precision, and there's
no way to disable this misfeature at compile time. Hence, it's
impossible to generate the appropriate tables of constants for the
long double inverse trig functions in a straightforward way on i386;
this change hacks around the problem by encoding the underlying bits
in the table.
Note that these functions won't pass the regression test on i386,
even with the FPU set to extended precision, because the regression
test is similarly damaged by gcc. However, the tests all pass when
compiled with a modified version of gcc.
Reported by: bde
- Adjust several constants for float precision. Some thresholds
that were appropriate for double precision were never changed
when these routines were converted to float precision. This
has an impact on performance but not accuracy. (Submitted by bde.)
- Reduce the degrees of the polynomials used. A smaller degree
suffices for float precision.
- In asinf(), use double arithmetic in part of the calculation to
avoid a corner case and some complicated arithmetic involving a
division and some buggy constants. This improves performance and
accuracy.
Max error (ulps):
asinf acosf atanf
before 0.925 0.782 0.852
after 0.743 0.804 0.852
As bde points out, it's cheaper for asin*() and acos*() to use
polynomials instead of rational functions, but that's a task for
another day.
spurious optimizations. gcc doesn't support FENV_ACCESS, so when it
folds constants, it assumes that the rounding mode is always the
default and floating point exceptions never matter.
1. architecture-specific files
2. long double format-specific files
3. bsdsrc
4. src
5. man
The original order was virtually the opposite of this.
This should not cause any functional changes at this time. The
difference is only significant when one wants to override, say, a
generic foo.c with a more specialized foo.c (as opposed to foo.S).
- fma(x, y, z) returns z, not NaN, if z is infinite, x and y are finite,
x*y overflows, and x*y and z have opposite signs.
- fma(x, y, z) doesn't generate an overflow, underflow, or inexact exception
if z is NaN or infinite, as per IEEE 754R.
- If the rounding mode is set to FE_DOWNWARD, fma(1.0, 0.0, -0.0) is -0.0,
not +0.0.
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.
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.)
e_rem_pio2.c:
This case goes up to about 2**20pi/2, but the comment about it said that
it goes up to about 2**19pi/2.
It went too far above 2**pi/2, giving a multiplier fn with 21 significant
bits in some cases. This would be harmful except for a numerical
accident. It happens that the terms of the approximation to pi/2,
when rounded to 33 bits so that multiplications by 20-bit fn's are
exact, happen to be rounded to 32 bits so multiplications by 21-bit
fn's are exact too, so the bug only complicates the error analysis (we
might lose a bit of accuracy but have bits to spare).
e_rem_pio2f.c:
The bogus comment in e_rem_pio2.c was copied and the code was changed
to be bug-for-bug compatible with it, except the limit was made 90
ulps smaller than necessary. The approximation to pi/2 was not
modified except for discarding some of it.
The same rough error analysis that justifies the limit of 2**20pi/2
for double precision only justifies a limit of 2**18pi/2 for float
precision. We depended on exhaustive testing to check the magic numbers
for float precision. More exaustive testing shows that we can go up
to 2**28pi/2 using a 53+25 bit approximation to pi/2 for float precision,
with a the maximum error for cosf() and sinf() unchanged at 0.5009
ulps despite the maximum error in rem_pio2f being ~0.25 ulps. Implement
this.
gives an average speedup of about 12 cycles or 17% for
9pi/4 < |x| <= 2**19pi/2 and a smaller speedup for larger x, and a
small speeddown for |x| <= 9pi/4 (only 1-2 cycles average, but that
is 4%).
Inlining this is less likely to bust caches than inlining the float
version since it is much smaller (about 220 bytes text and rodata) and
has many fewer branches. However, the float version was already large
due to its manual inlining of the branches and also the polynomial
evaluations.
__kernel_rem_pio2(). This simplifies analysis of aliasing and thus
results in better code for the usual case where __kernel_rem_pio2()
is not called. In particular, when __ieee854_rem_pio2[f]() is inlined,
it normally results in y[] being returned in registers. I couldn't
get this to work using the restrict qualifier.
In float precision, this saves 2-3% in most cases on amd64 and i386
(A64) despite it not being inlined in float precision yet. In double
precision, this has high variance, with an average gain of 2% for
amd64 and 0.7% for i386 (but a much larger gain for usual cases) and
some losses.
this function and its callers cosf(), sinf() and tanf() don't waste time
converting values from doubles to floats and back for |x| > 9pi/4.
All these functions were optimized a few years ago to mostly use doubles
internally and across the __kernel*() interfaces but not across the
__ieee754_rem_pio2f() interface.
This saves about 40 cycles in cosf(), sinf() and tanf() for |x| > 9pi/4
on amd64 (A64), and about 20 cycles on i386 (A64) (except for cosf()
and sinf() in the upper range). 40 cycles is about 35% for |x| < 9pi/4
<= 2**19pi/2 and about 5% for |x| > 2**19pi/2. The saving is much
larger on amd64 than on i386 since the conversions are not easy to
optimize except on i386 where some of them are automatic and others
are optimized invalidly. amd64 is still about 10% slower in cosf()
and tanf() in the lower range due to conversion overhead.
This also gives a tiny speedup for |x| <= 9pi/4 on amd64 (by simplifying
the code). It also avoids compiler bugs and/or additional slowness
in the conversions on (not yet supported) machines where double_t !=
double.
e_rem_pio2.c:
Float and double precision didn't work because init_jk[] was 1 too small.
It needs to be 2 larger than you might expect, and 1 larger than it was
for these precisions, since its test for recomputing needs a margin of
47 bits (almost 2 24-bit units).
init_jk[] seems to be barely enough for extended and quad precisions.
This hasn't been completely verified. Callers now get about 24 bits
of extra precision for float, and about 19 for double, but only about
8 for extended and quad. 8 is not enough for callers that want to
produce extra-precision results, but current callers have rounding
errors of at least 0.8 ulps, so another 1/2**8 ulps of error from the
reduction won't affect them much.
Add a comment about some of the magic for init_jk[].
e_rem_pio2.c:
Double precision worked in practice because of a compensating off-by-1
error here. Extended precision was asked for, and it executed exactly
the same code as the unbroken double precision.
e_rem_pio2f.c:
Float precision worked in practice because of a compensating off-by-1
error here. Double precision was asked for, and was almost needed,
since the cosf() and sinf() callers want to produce extra-precision
results, at least internally so that their error is only 0.5009 ulps.
However, the extra precision provided by unbroken float precision is
enough, and the double-precision code has extra overheads, so the
off-by-1 error cost about 5% in efficiency on amd64 and i386.
variations (e500 currently), this provides a gcc-level FPU emulation and is an
alternative approach to the recently introduced kernel-level emulation
(FPU_EMU).
Approved by: cognet (mentor)
MFp4: e500
fabs(), a conditional branch, and sign adjustments of 3 variables for
x < 0 when the branch is taken. In double precision, even when the
branch is perfectly predicted, this saves about 10 cycles or 10% on
amd64 (A64) and i386 (A64) for the negative half of the range, but
makes little difference for the positive half of the range. In float
precision, it also saves about 4 cycles for the positive half of the
range on i386, and many more cycles in both halves on amd64 (28 in the
negative half and 11 in the positive half for tanf), but the amd64
times for float precision are anomalously slow so the larger
improvement is only a side effect.
Previous commits arranged for the x < 0 case to be handled simply:
- one part of the rounding method uses the magic number 0x1.8p52
instead of the usual 0x1.0p52. The latter is required for large |x|,
but it doesn't work for negative x and we don't need it for large |x|.
- another part of the rounding method no longer needs to add `half'.
It would have needed to add -half for negative x.
- removing the "quick check no cancellation" in the double precision
case removed the need to take the absolute value of the quadrant
number.
Add my noncopyright in e_rem_pio2.c
FP-to-FP method to round to an integer on all arches, and convert this
to an int using FP-to-integer conversion iff irint() is not available.
This is cleaner and works well on at least ia64, where it saves 20-30
cycles or about 10% on average for 9Pi/4 < |x| <= 32pi/2 (should be
similar up to 2**19pi/2, but I only tested the smaller range).
After the previous commit to e_rem_pio2.c removed the "quick check no
cancellation" non-optimization, the result of the FP-to-integer
conversion is not needed so early, so using irint() became a much
smaller optimization than when it was committed.
An earlier commit message said that cos, cosf, sin and sinf were equally
fast on amd64 and i386 except for cos and sin on i386. Actually, cos
and sin on amd64 are equally fast to cosf and sinf on i386 (~88 cycles),
while cosf and sinf on amd64 are not quite equally slow to cos and sin
on i386 (average 115 cycles with more variance).
9pi/2 < |x| < 32pi/2 since it is only a small or negative optimation
and it gets in the way of further optimizations. It did one more
branch to avoid some integer operations and to use a different
dependency on previous results. The branches are fairly predictable
so they are usually not a problem, so whether this is a good
optimization depends mainly on the timing for the previous results,
which is very machine-dependent. On amd64 (A64), this "optimization"
is a pessimization of about 1 cycle or 1%; on ia64, it is an
optimization of about 2 cycles or 1%; on i386 (A64), it is an
optimization of about 5 cycles or 4%; on i386 (Celeron P2) it is an
optimization of about 4 cycles or 3% for cos but a pessimization of
about 5 cycles for sin and 1 cycle for tan. I think the new i386
(A64) slowness is due to an pipeline stall due to an avoidable
load-store mismatch (so the old timing was better), and the i386
(Celeron) variance is due to its branch predictor not being too good.
the the double to int conversion operation which is very slow on these
arches. Assume that the current rounding mode is the default of
round-to-nearest and use rounding operations in this mode instead of
faking this mode using the round-towards-zero mode for conversion to
int. Round the double to an integer as a double first and as an int
second since the double result is needed much earler.
Double rounding isn't a problem since we only need a rough approximation.
We didn't support other current rounding modes and produce much larger
errors than before if called in a non-default mode.
This saves an average about 10 cycles on amd64 (A64) and about 25 on
i386 (A64) for x in the above range. In some cases the saving is over
25%. Most cases with |x| < 1000pi now take about 88 cycles for cos
and sin (with certain CFLAGS, etc.), except on i386 where cos and sin
(but not cosf and sinf) are much slower at 111 and 121 cycles respectivly
due to the compiler only optimizing well for float precision. A64
hardware cos and sin are slower at 105 cycles on i386 and 110 cycles
on amd64.
the same as lrint() except it returns int instead of long. Though the
extern lrint() is fairly fast on these arches, it still takes about
12 cycles longer than the inline version, and 12 cycles is a lot in
applications where [li]rint() is used to avoid slow conversions that
are only a couple of times slower.
This is only for internal use. The libm versions of *rint*() should
also be inline, but that would take would take more header engineering.
Implementing irint() instead of lrint() also avoids a conflict with
the extern declaration of the latter.
on i386 (A64), 5 cycles on amd64 (A64), and 3 cycles on ia64). gcc
tends to generate very bad code for accessing floating point values
as bits except when the integer accesses have the same width as the
floating point values, and direct accesses to bit-fields (as is common
only for long double precision) always gives such accesses. Use the
expsign access method, which is good for 80-bit long doubles and
hopefully no worse for 128-bit long doubles. Now the generated code
is less bad. There is still unnecessary copying of the arg on amd64
and i386 and mysterious extra slowness on amd64.
