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.
particularly good reason to do this, except that __strong_reference
does type checking, whereas __weak_reference does not.
On Alpha, the compiler won't accept a 'long double' parameter in
place of a 'double' parameter even thought the two types are
identical.
an invalid exception and return an NaN.
- If a long double has 113 bits of precision, implement fma in terms
of simple long double arithmetic instead of complicated double arithmetic.
- If a long double is the same as a double, alias fma as fmal.
identical to scalbnf, which is now aliased as ldexpf. Note that the
old implementations made the mistake of setting errno and were the
only libm routines to do so.
- Add nexttoward{,f,l} and nextafterl. On all platforms,
nexttowardl is an alias for nextafterl.
- Add fmal.
- Add man pages for new routines: fmal, nextafterl,
nexttoward{,f,l}, scalb{,l}nl.
Note that on platforms where long double is the same as double, we
generally just alias the double versions of the routines, since doing
so avoids extra work on the source code level and redundant code in
the binary. In particular:
ldbl53 ldbl64/113
fmal s_fma.c s_fmal.c
ldexpl s_scalbn.c s_scalbnl.c
nextafterl s_nextafter.c s_nextafterl.c
nexttoward s_nextafter.c s_nexttoward.c
nexttowardf s_nexttowardf.c s_nexttowardf.c
nexttowardl s_nextafter.c s_nextafterl.c
scalbnl s_scalbn.c s_scalbnl.c
sparc64's 128-bit long doubles.
- Define FP_FAST_FMAL for ia64.
- Prototypes for fmal, frexpl, ldexpl, nextafterl, nexttoward{,f,l},
scalblnl, and scalbnl.
- In scalbln and scalblnf, check the bounds of the second argument.
This is probably unnecessary, but strictly speaking, we should
report an error if someone tries to compute scalbln(x, INT_MAX + 1ll).
nexttowardl. These are not needed on machines where long doubles
look like IEEE-754 doubles, so the implementation only supports
the usual long double formats with 15-bit exponents.
Anything bizarre, such as machines where floating-point and integer
data have different endianness, will cause problems. This is the case
with big endian ia64 according to libc/ia64/_fpmath.h. Please contact
me if you managed to get a machine running this way.
that are intended to raise underflow and inexact exceptions.
- On systems where long double is the same as double, nextafter
should be aliased as nexttoward, nexttowardl, and nextafterl.
Symptoms of the problem included assembler warnings and
nondeterministic runtime behavior when a fe*() call that affects the
fpsr is closely followed by a float point op.
The bug (at least, I think it's a bug) is that gcc does not insert a
break between a volatile asm and a dependent instruction if the
volatile asm came from an inlined function. Volatile asms seem to be
fine in other circumstances, even without -mvolatile-asm-stop, so
perhaps the compiler adds the stop bits before inlining takes place.
The problem does not occur at -O0 because inlining is disabled, and it
doesn't happen at -O2 because -fschedule-insns2 knows better.
inputs. The trouble with replacing two floats with a double is that
the latter has 6 extra bits of precision, which actually hurts
accuracy in many cases. All of the constants are optimal when float
arithmetic is used, and would need to be recomputed to do this right.
Noticed by: bde (ucbtest)
results in a performance gain on the order of 10% for amd64 (sledge),
ia64 (pluto1), i386+SSE (Pentium 4), and sparc64 (panther), and a
negligible improvement for i386 without SSE. (The i386 port still
uses the hardware instruction, though.)
manpages. They are not very related, so separating them makes it
easier to add meaningful cross-references and extend some of the
descriptions.
- Move the part of math(3) that discusses IEEE 754 to the ieee(3)
manpage.
- Rearrange the list of functions into categories.
- Remove the ulps column. It was appropriate for only some
of the functions in the list, and correct for even fewer
of them.
- Add some new paragraphs, and remove some old ones about
NaNs that may do more harm than good.
- Document precisions other than double-precision.
flags, so they are not pure. Remove the __pure2 annotation from them.
I believe that the following routines and their float and long double
counterparts are the only ones here that can be __pure2:
copysign is* fabs finite fmax fmin fpclassify ilogb nan signbit
When gcc supports FENV_ACCESS, perhaps there will be a new annotation
that allows the other functions to be considered pure when FENV_ACCESS
is off.
Discussed with: bde
basically support this, subject to gcc's lack of FENV_ACCESS support.
In any case, the previous setting of math_errhandling to 0 is not
allowed by POSIX.
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.
- Make some minor rearrangements in the introduction.
- Mention the problem with argument reduction on i386.
- Add recently-implemented functions to the table.
- Un-document the error bounds that only apply to the old 4BSD math
library, and fill in the correct values where I know them. No
attempt has been made to document bounds lower than 1 ulp, although
smaller bounds are usually achievable in round-to-nearest mode.
/lib/{libm,libreadline}
/usr/lib/{libhistory,libopie,libpcap}
in preparation for doing the same thing to RELENG_5. HUGE amounts of
help for determining what to bump provided by kris.
Discussed on: freebsd-current
Approved by: re (not required for commit but something like this should be)
libc. The externally-visible effect of this is to add __isnanl() to
libm, which means that libm.so.2 can once again link against libc.so.4
when LD_BIND_NOW is set. This was broken by the addition of fdiml(),
which calls __isnanl().
- It was added to libc instead of libm. Hopefully no programs rely
on this mistake.
- It didn't work properly on large long doubles because its argument
was converted to type double, resulting in undefined behavior.
- Unlike the builtin relational operators, builtin floating-point
constants were not available until gcc 3.3, so account for this.[1]
- Apparently some versions of the Intel C Compiler fallaciously define
__GNUC__ without actually being compatible with the claimed gcc
version. Account for this, too.[2]
[1] Noticed by: Christian Hiris <4711@chello.at>
[2] Submitted by: Alexander Leidinger <Alexander@Leidinger.net>
isnormal() the hard way, rather than relying on fpclassify(). This is
a lose in the sense that we need a total of 12 functions, but it is
necessary for binary compatibility because we have never bumped libm's
major version number. In particular, isinf(), isnan(), and isnanf()
were BSD libc functions before they were C99 macros, so we can't
reimplement them in terms of fpclassify() without adding a dependency
on libc.so.5. I have tried to arrange things so that programs that
could be compiled in FreeBSD 4.X will generate the same external
references when compiled in 5.X. At the same time, the new macros
should remain C99-compliant.
The isinf() and isnan() functions remain in libc for historical
reasons; however, I have moved the functions that implement the macros
isfinite() and isnormal() to libm where they belong. Moreover,
half a dozen MD versions of isinf() and isnan() have been replaced
with MI versions that work equally well.
Prodded by: kris
builtins are available: HUGE_VAL, HUGE_VALF, HUGE_VALL, INFINITY,
and NAN. These macros now expand to floating-point constant
expressions rather than external references, as required by C99.
Other compilers will retain the historical behavior. Note that
it is not possible say, e.g.
#define HUGE_VAL 1.0e9999
because the above may result in diagnostics at translation time
and spurious exceptions at runtime. Hence the need for compiler
support for these features.
Also use builtins to implement the macros isgreater(),
isgreaterequal(), isless(), islessequal(), islessgreater(),
and isunordered() when such builtins are available.
Although the old macros are correct, the builtin versions
are much faster, and they avoid double-expansion problems.
These trivial implementations are about 25 times slower than
rint{,f}() on x86 due to the FP environment save/restore.
They should eventually be redone in terms of fegetround() and
bit fiddling.
These routines are specified in C99 for the sake of
architectures where an int isn't big enough to represent
the full range of floating-point exponents. However,
even the 128-bit long double format has an exponent smaller
than 15 bits, so for all practical purposes, scalbln() and
scalblnf() are aliases for scalbn() and scalbnf(), respectively.
kicking and screaming into the 1980's. This change converts most of
the markup from man(7) to mdoc(7) format, and I believe it removes or
updates everything that was flat out wrong. However, much work is
still needed to sanitize the markup, improve coverage, and reduce
overlap with other manpages. Some of the sections would better belong
in a philosophy_of_w_kahan.3 manpage, but they are informative and
remain at least as reminders of topics to cover.
Reviewed by: doc@, trhodes@
on all inputs of the form x.75, where x is an even integer and
log2(x) = 21. A similar problem occurred when rounding upward.
The bug involves the following snippet copied from rint():
i>>=1;
if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0);
The constant 0x100000 should be 0x200000. Apparently this case was
never tested.
It turns out that the bit manipulation is completely superfluous
anyway, so remove it. (It tries to simulate 90% of the rounding
process that the FPU does anyway.) Also, the special case of +-0 is
handled twice (in different ways), so remove the second instance.
Throw in some related simplifications from bde:
- Work around a bug where gcc fails to clip to float precision by
declaring two float variables as volatile. Previously, we
tricked gcc into generating correct code by declaring some
float constants as doubles.
- Remove additional superfluous bit manipulation.
- Minor reorganization.
- Include <sys/types.h> explicitly.
Note that some of the equivalent lines in rint() also appear to be
unnecessary, but I'll defer to the numerical analysts who wrote it,
since I can't test all 2^64 cases.
Discussed with: bde
It does not appear to be possible to cross-build arm from i386 at the
moment, and I have no ARM hardware anyway. Thus, I'm sure there are
bugs. I will gladly fix these when the arm port is more mature.
Reviewed by: standards@
features appear to work, subject to the caveat that you tell gcc you
want standard rather than recklessly fast behavior
(-mieee-with-inexact -mfp-rounding-mode=d).
The non-standard feature of delivering a SIGFPE when an application
raises an unmasked exception does not work, presumably due to a kernel
bug. This isn't so bad given that floating-point exceptions on the
Alpha architecture are not precise, so making them useful in userland
requires a significant amount of wizardry.
Reviewed by: standards@
