functions are only for compatibility with obsolete standards. They
shouldn't be used, so they shouldn't be optimized. Use the generic
versions instead.
This fixes scalbf() as a side effect. The optimized asm version left
garbage on the FP stack. I fixed the corresponding bug in the optimized
asm scalb() and scalbn() in 1996. NetBSD fixed it in scalb(), scalbn()
and scalbnf() in 1999 but missed fixing it in scalbf(). Then in 2005
the bug was reimplemented in FreeBSD by importing NetBSD's scalbf().
The generic versions have slightly different error handling:
- the asm versions blindly round the second parameter to a (floating
point) integer and proceed, while the generic versions return NaN
if this rounding changes the value. POSIX permits both behaviours
(these functions are XSI extensions and the behaviour for a bogus
non-integral second parameter is unspecified). Apart from this
and the bug in scalbf(), the behaviour of the generic versions seems
to be identical. (I only exhusatively tested
generic_scalbf(1.0F, anyfloat) == asm_scalb(1.0F, anyfloat). This
covers many representative corner cases involving NaNs and Infs but
doesn't test exception flags. The brokenness of scalbf() showed up
as weird behaviour after testing just 7 integer cases sequentially.)
The log message for 1.5 said that some small (one or two ulp) inaccuracies
were fixed, and a comment implied that the critical change is to switch
the rounding mode to to-nearest, with a switch of the precision to
extended at no extra cost. Actually, the errors are very large (ucbtest
finds ones of several hundred ulps), and it is the switch of the
precision that is critical.
Another comment was wrong about NaNs being handled sloppily.
floating-point arithmetic on i386. Now I'm going to make excuses
for why this code is kinda scary:
- To avoid breaking the ABI with 5.3-RELEASE, we can't change
sizeof(fenv_t). I stuck the saved mxcsr in some discontiguous
reserved bits in the existing structure.
- Attempting to access the mxcsr on older processors results
in an illegal instruction exception, so support for SSE must
be detected at runtime. (The extra baggage is optimized away
if either the application or libm is compiled with -msse{,2}.)
I didn't run tests to ensure that this doesn't SIGILL on older 486's
lacking the cpuid instruction or on other processors lacking SSE.
Results from running the fenv regression test on these processors
would be appreciated. (You'll need to compile the test with
-DNO_STRICT_DFL_ENV.) If you have an 80386, or if your processor
supports SSE but the kernel didn't enable it, then you're probably out
of luck.
Also, I un-inlined some of the functions that grew larger as a result
of this change, moving them from fenv.h to fenv.c.
fedisableexcept(), and fegetexcept(). These two sets of routines
provide the same functionality. I implemented the former as an
undocumented internal interface to make the regression test easier to
write. However, fe(enable|disable|get)except() is already part of
glibc, and I would like to avoid gratuitous differences. The only
major flaw in the glibc API is that there's no good way to report
errors on processors that don't support all the unmasked exceptions.
registers as volatile. Instructions that *wrote* to FP state were
already marked volatile, but apparently gcc has license to move
non-volatile asms past volatile asms. This broke amd64's feupdateenv
at -O2 due to a WAR conflict between fnstsw and fldenv there.
avoid easily avoidable loss of precision when |x| is nearly 1.
Extended (64-bit) precision only moves the meaning of "nearly" here.
This probably could be done better by splitting up the range into
|x| <= 0.5 and |x| > 0.5 like the C version. However, ucbtest
does't report any errors in this version. Perhaps the C version
should be used anyway. It's only 25% slower now on a P5, provided
the C version of sqrt() isn't used, and the C version could be
optimized better.
Errors checked by: ucbtest
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.
The fyl2xp1 instruction has such a limited range:
-(1 - (sqrt(2) / 2)) <= x <= sqrt(2) - 1
it's not worth trying to use it.
Also, I'm not sure fyl2xp1's extra precision will
matter once the result is converted from extended
real (80 bits) back to double real (64 bits).
Reviewed by: jkh
Submitted by: jtc
-- 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
