libm: fmod, remquo, remainder are cr
[www] / libm / index.html
index 5832554..dce0a91 100644 (file)
@@ -1,14 +1,12 @@
 <html><head><title>libm</title></head><body>
 <h2>libm</h2>
 
-<p>This page is about designing libm for the
-<a href="http://www.etalabs.net/musl/">musl</a> libc.
-<p>Writing the math code from scratch is a huge work
-so already existing code is used.
+<p>This page is about libm for the
+<a href="http://www.musl-libc.org/">musl</a> libc.
 
 <ul>
 <li><a href="#sources">Sources</a>
-<li><a href="#issues">Design issues</a>
+<li><a href="#rules">General rules</a>
 <li><a href="#representation">Representation</a>
 <li><a href="#ugly">Ugly</a>
 <li><a href="#implementations">libm implementations</a>
@@ -16,8 +14,13 @@ so already existing code is used.
 </ul>
 
 <h3><a name="sources" href="#sources">Sources</a></h3>
-<p>The math code is mostly from freebsd which in turn
-is based on <a href="http://www.netlib.org/fdlibm/">fdlibm</a>.
+<p>Writing math code from scratch is a huge work so already existing code is
+used. Several math functions are taken from the
+<a href="http://freebsd.org/">freebsd</a> libm and a few from the
+<a href="http://openbsd.org/">openbsd</a> libm implementations.
+Both of them are based on <a href="http://www.netlib.org/fdlibm/">fdlibm</a>.
+The freebsd libm seems to be the most well maintained and most correct version
+of fdlibm.
 <p>sources:
 <ul>
 <li>freebsd /lib/msun (<a href="http://svnweb.FreeBSD.org/base/head/lib/msun/">browse</a>)
@@ -30,91 +33,108 @@ cvs -d <a href="http://openbsd.org/anoncvs.html#CVSROOT">$CVSROOT</a> get src/li
 </pre>
 </ul>
 
-<h3><a name="issues" href="#issues">Design issues</a></h3>
-<p>The math code is available but there are still many design questions.
+<h3><a name="rules" href="#rules">General rules</a></h3>
 <ul>
-<li>Code organization, build system tweaks:
-<p>There are architecture specific code,
-long double representation specific code
-and complex code.
-<p>With respect to long double, bsd libm code can handle
-ld64, ld80 and ld128 (where ld64 means long double == double),
-for musl it's enough to support ld64 and ld80, but
-keeping ld128 shouldn't hurt much either.
-<p>So far musl had arch specific code and common code,
-it seems a new category is needed:
-long double specific code, which can be ld64, ld80 or ld128.
-These should be tied to the architecture in the build system
-somehow. In freebsd most of the long double code is
-common (using arch specific macros, ifdefs) and only a
-small part of the ld code is maintained separately for ld80 and ld128,
-but it's also possible to keep all ld code separately (with duplications)
-in ld80 and ld128.
-<p>Complex support is optional in c99 so maybe it should
-be optional in the musl build as well. It also can be
-kept separately in cmath/ or with the rest of the code in math/.
-<p>Pending questions:
+<li>Assumptions about floating-point representation and arithmetics
+(see <a href="http://port70.net/~nsz/c/c99/n1256.html#F.2">c99 annex F.2</a>):
 <ul>
-<li>Do we want ld128?
-<li>Should we try to use common code for ld80 and ld128?
-<li>How to do ld64: wrap double functions or alias them?
-<li>How to tie the ld* code to the arch in the build system?
-<li>Make complex optional?
-<li>Keep complex in math/ or cmath/?
-</ul>
-
-<li>Workarounds
-<p>The bsd libm code has many workarounds for various
-compiler issues. It's not clear what's the best way
-to handle the uglyness.
-<p>Pending questions:
-<ul>
-<li>Use STDC pragmas (eventhough gcc does not support them)?
-<li>Use volatile consistently to avoid evaluation precision and const folding issues?
-<li>Use special compiler flags against unwanted optimization (-frounding-math, -ffloat-store)?
-<li>Do inline/macro optimization for small functions? (isnan, isinf, signbit, creal, cimag,..)
-<li>In complex code prefer creal(), cimag() or a union to (un)pack re,im?
-</ul>
-
-<li>Code cleanups
-<p>The fdlibm code style is inconsistent with musl (and ugly)
-so it makes sense to clean it up, but keeping easy diffability
-can be useful. There are various inconsistencies and
-correctness issues in the code which might worth addressing
-(eg. reliance on signed overflow should be eliminated).
-<p>Pending questions:
-<ul>
-<li>Keep diffability with freebsd/openbsd code or reformat the code for clarity?
-<li>Keep e_, s_, k_  fdlibm source file name prefixes?
-<li>Should 0x1p0 float format be preferred over decimal format?
-<li>Should isnan, signbit,.. be preferred over inplace bithacks?
-<li>Is unpacking a double into 2 int32_t ok (signed int representation)?
-<li>Is unpacking the mantissa of ld80 into an int64_t ok?
+<li>float is ieee binary32
+<li>double is ieee binary64
+<li>long double is either ieee binary64 or little-endian 80bit extended precision (x87 fpu)
 </ul>
+(other long double representations may be supported in the future, until then
+long double math functions will be missing on non-supported platforms)
+<li>On soft-float architectures fenv should not be expected to work according to
+the c and ieee standards (ie. rounding modes and exceptions are not supported,
+the fenv functions are dummy ones)
+<li>Floating-point exception flags should be properly set in math functions
+according to c99 annex F, but without using fenv.h
+(eg. overflow flag can be raised by 0x1p900*0x1p900, because this works even
+without fenv support)
+<li>Most functions need to be precise only in nearest rounding mode.
+<li>Returned floating-point values should be correctly rounded in most cases,
+but the last bit may be wrong:
+<pre>
+    |error| &lt; 1.5 ulp
+</pre>
+should hold for most functions.
+(error is the difference between the exact result and the calculated
+floating-point value)
+(in theory correct rounding can be achieved but with big implementation cost,
+see <a href="http://lipforge.ens-lyon.fr/www/crlibm/">crlibm</a>)
+<li>At least the following functions must be correctly rounded:
+ceil, copysign, fabs, fdim, floor, fma, fmax, fmin, fmod, frexp, ldexp, logb,
+modf, nearbyint, nextafter, nexttoward, rint, remainder, remquo, round, scalbln,
+scalbn, sqrt, trunc.
+<li>Mathematical properties of functions should be as expected
+(monotonicity, range, symmetries).
+<li>If the FPU precision is altered then nothing is guaranteed to work.
+(ie. when long double does not have full 80bit precision on i386 then things
+may break, this also means that compiler must spill fpu registers at 80bits
+precision)
+<li>Signaling NaN is not supported
+<li>Quiet NaN is supported but all NaNs are treated equally without special
+attention to the internal representation of a NaN
+(eg. the sign of NaN may not be preserved).
+<li>Most gcc bug workarounds should be removed from the code
+(STRICT_ASSIGN macro is used when excessive precision is harmful and
+FORCE_EVAL when expressions must be evaluated for their side-effect, other
+usage of volatile is not justified, hacks around long double constants are
+not justified eventhough gcc can miscompile those with non-default FPU setting)
+<li>When excessive precision is not harmful, temporary variables
+should be float_t or double_t (so on i386 no superfluous store is
+generated)
+<li>Whenever fenv is accessed the FENV_ACCESS pragma of c99 should be used
+(eventhough gcc does not yet support it), and all usage of optional FE_
+macros should be protected by #ifdef
+<li>For bit manipulation of floating-point values an union should be used
+(eg. union {float f; uint32_t i;})
+<li>uint32_t and uint64_t should be used for bit manipulations.
+(eg signed int must not be used in bit shifts etc when it might invoke
+undefined or implementation defined behaviour).
+<li>POSIX namespace rules must be respected.
+<li>c99 hexfloat syntax (0x1.0p0) should be used when it makes the code
+clearer, but not in public header files
+(those should be c++ and ansi c compatible)
+<li>The 'f' suffix should be used for single precision values (0.1f) when the
+value cannot be exactly represented or the type of the arithmetics is important
+((float)0.1 is not ok, that style may lead to double rounding issues, but eg.
+1.0 or 0.5 may be used instead of 1.0f or 0.5f in some cases)
+<li>Prefer classification macros (eg. isnan) over inplace bit hacks.
+<li>For every math function there should be a c implementation.
+(a notable exception now is sqrtl, since most fpu has instruction for it
+and on soft-float architectures long double == double)
+<li>The c implementation of a long double function should use ifdefs with the
+LDBL_MANT_DIG etc constants from float.h for architecture specific
+implementations.
+<li>In the musl source tree math.h functions go to src/math, complex.h functions
+to src/complex and fenv.h functions to src/fenv. And files are named after the
+functions they implement.
 </ul>
 
 <h3><a name="representation" href="#representation">Representation</a></h3>
 <p>
 Binary representation of floating point numbers matter
-because lots of bithacks are needed in the math code.
+because bit hacks are often needed in the math code.
+(in particular bit hacks are used instead of relational operations for nan
+and sign checks becuase relational operators raise invalid fp exception on nan
+and they treat -0.0 and +0.0 equally and more often than not these are not desired)
 <p>
-float and double bit manipulation can be handled
-in a portable way in c:
+float and double bit manipulation can be handled in a portable way in c using
+union types:
 <ul>
-<li>float is ieee binary32
-<li>double is ieee binary64
+<li>union {float f; uint32_t i;};
+<li>union {double f; uint64_t i;};
 </ul>
-(and old non-standard representations are not supported)
-<p>
-The endianness may still vary, but that can be worked
-around by using a union with a single large enough
-unsigned int. (The only exception is probably arm/oabi fpa
-where the word order of a double is not consistent with the
-endianness [<a href="http://wiki.debian.org/ArmEabiPort#ARM_floating_points">debian wiki on arm</a>],
-but we probably won't support that)
-<p>
-long double bit manipulation is harder as there are
-various representations:
+(assuming the bits in the object representation of 32bit and 64bit unsigned ints
+map to the floating-point representation according to ieee-754, this is not
+always the case, eg. old
+<a href="http://wiki.debian.org/ArmEabiPort#ARM_floating_points">arm floating-point accelerator</a>
+(FPA) used mixed endian double representation, but musl does not support the old
+arm ABI)
+<p>
+long double bit manipulation is harder as there are various representations
+and some of them don't map to any unsigned integer type:
 <ul>
 <li>ld64: long double is the same as double (ieee binary64)
 <li>ld80: 80bit extended precision format [<a href="https://en.wikipedia.org/wiki/Extended_precision">wikipedia</a>]
@@ -122,31 +142,50 @@ various representations:
 </ul>
 (and other non-standard formats are not supported)
 <p>
-ld64 is easy to handle: all long double functions
-are aliased to the corresponding double one
-(long double is treated equivalently to double)
+In case of ld64 the bit manipulation is the same as with double
+and all long double math functions can be just wrappers around the
+corresponding double ones.
+(using symbol aliasing on the linker level is non-conformant
+since functions would not have unique address then)
 <p>
-ld80 is the most common (i386, x86_64), it means
-64bit significand with explicit msb (inconsistent with other ieee formats),
-15bit exp, 1 sign bit.
+ld80 is the most common long double on linux (i386 and x86_64 abi),
+it means 64bit significand with explicit msb
+(inconsistent with other ieee formats), 15bit exp, 1 sign bit.
+The m68k (and m88k) architecture uses the same format, but different endianness:
+<ul>
+<li>union {long double f; struct{uint64_t m; uint16_t se; uint16_t pad;} i;}; // x86
+<li>union {long double f; struct{uint16_t se; uint16_t pad; uint64_t m;} i;}; // m68k
+</ul>
+where m is the significand and se is the sign and exponent.
 <p>
-ld128 is rare (sparc64 with software emulation), it means
-113bit significand with implicit msb, 15bit exp, 1 sign bit.
+ld128 is rare (eg. sparc64 with software emulation), it means
+113bit significand with implicit msb, 15bit exp, 1 sign bit:
+<ul>
+<li>union {long double f; struct{uint16_t se; uint16_t hi; uint32_t mid; uint64_t lo;} i;};
+</ul>
 <p>
-Endianness can vary (although on the supported i386 and x86_64 it is the same)
-and there is no large enough unsigned int to handle it.
-(In case of ld80 internal padding can vary as well, eg
-m68k and m88k cpus use different ld80 than the intel ones)
-
+There are other non-conformant long double types: eg. the old SVR4 abi for ppc
+uses 128 bit long doubles, but it's software emulated and traditionally
+implemented using
+<a href="https://en.wikipedia.org/wiki/Quadruple_precision#Double-double_arithmetic">two doubles</a>
+(also called ibm long double as this is what ibm aix used on ppc).
+The ibm s390 supports the ieee 754-2008 compliant binary128 floating-point
+format, but previous ibm machines (S/370, S/360) used slightly different
+representation.
+<p>
+This variation shows the difficulty to consistently handle
+long double: the solution is to use ifdefs based on float.h and
+on the endianness and write different code for different architectures.
 
 <h3><a name="ugly" href="#ugly">Ugly</a></h3>
 <p>The ugly parts of libm hacking.
 <p>Some notes are from:
 <a href="http://www.vinc17.org/research/extended.en.html">http://www.vinc17.org/research/extended.en.html</a>
-
+<p>Useful info about floating-point in gcc:
+<a href="http://gcc.gnu.org/wiki/FloatingPointMath">http://gcc.gnu.org/wiki/FloatingPointMath</a>
 
 <ul>
-<li>Double rounding (x87 issue):
+<li>Double rounding:
 <p>
 If a value rounded twice the result can be different
 than rounding just once.
@@ -158,16 +197,30 @@ the results twice: round to 80bit when calculating
 and then round to 64bit when storing it, this can
 give different result than a single 64bit rounding.
 (on x86-linux the default fpu setting is to round the
-results in extended precision, this only affects x87 instructions, not see2 etc)
-(afaik freebsd and openbsd use double precision by default)
+results in extended precision, this only affects x87 instructions, not sse2 etc)
+(freebsd and openbsd use double precision by default)
 <p>
 So x = a+b may give different results depending on
-the fpu setting.
+the x87 fpu precision setting.
 (only happens in round to nearest rounding mode,
 but that's the most common one)
 <p>
-C99 annex F prohibits double rounding,
-but that's non-normative.
+(double rounding can happen with float vs double as well)
+<p>
+<a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#F.7.3">C99 annex F</a>
+prohibits double rounding, but that's non-normative.
+<p>
+Note that the value of the result can only be ruined by
+double rounding in nearest rounding mode, but the double
+rounding issue haunts directed rounding modes as well:
+raising the underflow flag might be omitted.
+On x86 with downward rounding
+<pre>
+(double)(0x1p-1070 + 0x1p-2000L)
+</pre>
+does not raise underflow (only inexact) eventhough the
+final result is an inexact subnormal.
+
 
 <li>Wider exponent range (x87 issue):
 <p>
@@ -188,6 +241,9 @@ C does not require consistent evaluation
 precision: the compiler may store intermediate
 results and round them to double while keep
 other parts in higher precision.
+(And the precision the compiler choses can
+be inconsistent: adding a printf to the code
+may change the result of a nearby calculation).
 So
 <pre>
 (a+b)==(a+b)
@@ -195,9 +251,12 @@ So
 may be false when the two sides
 are kept in different precision.
 (This is not an x87 specific problem, it matters whenever there
-is a higher precision fp type than the currently used one.
+is a higher precision fp type than the currently used one and
+FLT_EVAL_METHOD!=0.
 It goes away if the highest precision (long double) is used
 everywhere, but that can have a huge penalty).
+(clang uses sse by default on i386 with FLT_EVAL_METHOD==0,
+while gcc uses the 80bit x87 fp registers and FLT_EVAL_METHOD==2)
 <p>
 C99 has a way to control this (see
 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#5.1.2.3">5.1.2.3 example 4</a>,
@@ -216,8 +275,7 @@ when the comparision is done.
 (This still does not solve the x87 double rounding
 issue though: eg if the left side is evaluated with
 sse2 and the right side with x87 extended precision setting
-and double rounding then the result may be false
-and still conformant to C99)
+and double rounding then the result may still be false)
 <p>
 Unfortunately gcc does not respect the standard
 and even if assingment or cast is used the result
@@ -227,13 +285,18 @@ also see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36578">gcc bug3657
 gcc 4.5 fixed it with '-fexcess-precision=standard'
 (it is enabled by '-std=c99', but the default is
 '-fexcess-precision=fast')
+(An alternative solution would be if gcc spilled the
+registers with temporary results without rounding,
+storing the 80 bit registers entirely in memory
+which would make the behaviour under FLT_EVAL_METHOD==2
+mode more predictable)
 <p>
 The workaround for older gcc is to force the
 compiler to store the intermediate results:
 by using volatile double temporary variables
-or by '-ffloat-store' (slow, and of course
-all intermediate results should be assigned
-to some variable, casting is not enough).
+or by '-ffloat-store' (which affects all
+intermediate results, but is not guaranteed
+by the gcc docs to always work).
 <p>
 (Sometimes the excess precision is good
 but it's hard to rely on it as it is optional
@@ -245,6 +308,40 @@ use higher precision variables when that's
 what we want and don't depend on the implicit
 excess precision).
 
+<li>Float literals
+<p>
+The standard allows 1 ulp errors in the conversion
+of decimal floating-point literals into floating-point
+values (it only requires the same result for the same
+literal thus <tt>1.1 - 1.1</tt> is always 0,
+but <tt>1.1 - 11e-1</tt> maybe +-0x1p-52 or 0).
+<p>
+A reasonable compiler always use correctly rounded
+conversion according to the default (nearest) rounding
+mode, but there are exceptions:
+the x87 has builtin constants which are faster to load
+from hardware than from memory
+(and the hw has sticky bit set correctly for rounding).
+gcc can recognize these constants so an operation on
+<pre>
+3.141592653589793238462643383L
+</pre>
+can turn into code that uses the <tt>fldpi</tt>
+instruction instead of memory loads.
+The only issue is that fldpi depends on
+the current rounding mode at runtime
+so the result can indeed be 1 ulp off compared
+to the compile-time rounded value.
+<p>
+According to the freebsd libm code gcc truncates long double
+const literals on i386.
+I assume this happens because freebsd uses 64bit long doubles by default
+(double precision) and gcc incorrectly uses the precision setting of the
+host platform instead of the target one, but i did not observe this on linux.
+(as a workaround sometimes double-double arithmetics was used
+to initialize long doubles on i386, but most of these should be
+fixed in musl's math code now)
+
 <li>Compiler optimizations:
 <p>
 Runtime and compile time semantics may be different
@@ -259,13 +356,15 @@ different precision than at runtime).
 C99 actually allows most of these optimizations
 but they can be turned off with STDC pragmas (see
 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#6.10.6">6.10.6</a>).
-Unfortunately <a href="http://gcc.gnu.org/c99status.html">gcc does not support these pragmas</s>.
+Unfortunately <a href="http://gcc.gnu.org/c99status.html">gcc does not support these pragmas</a>
+nor clang (<a href="http://llvm.org/bugs/show_bug.cgi?id=8100">clang bug 8100</a>).
 <p>
 FENV_ACCESS ON tells the compiler that the code wants
 to access the floating point environment (eg. set different rounding mode)
 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34678">gcc bug34678</a>).
 <p>
-(see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845">gcc bug37845</a> for FP_CONTRACT pragma).
+(see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845">gcc bug37845</a> for FP_CONTRACT pragma
+which is relevant when the architecture has fma instruction, x86_64, i386, arm do not have it).
 <p>
 The workaround is again using named volatile
 variables for constants like
@@ -273,37 +372,99 @@ variables for constants like
 static const volatile two52 = 0x1p52;
 </pre>
 and using the '-frounding-math' gcc flag.
-<p>
-(According the freebsd libm code gcc truncates
-long double const literals on i386.
-I haven't yet verified if this still the case,
-but as a workaround sometimes double-double arithmetics is used:
-initializing the long double constant from two doubles)
 
 <li>ld80 vs ld128
 <p>
 The two representations are sufficiently different
 that treating them together is awkward.
+(Especially the explicit msb bit in ld80 can cause
+different behaviour).
 <p>
 In the freebsd libm code a few architecture specific
 macros and a union handle these issues, but the
 result is often less clear than treating ld80
 and ld128 separately.
 
-<li>Signed int
-<p>
-The freebsd libm code has many inconsistencies
-(naming conventions, 0x1p0 notation vs decimal notation,..),
-one of them is the integer type used for bitmanipulations:
-The bits of a double are unpacked into one of
-int32_t, uint32_t and u_int32_t
-integer types.
+<li>Signed zeros
+<p>Signed zeros can be tricky. They cannot be checked
+using the usual comparision operators (+0.0 == -0.0 is true),
+but they give different results in some cases
+(1/-0.0 == -inf) and can make similarly looking
+expressions different eg.
+(x + 0) is not the same as x, the former is +0.0 when x is -0.0
+<p>
+(To check for -0, the signbit macro can be used
+or the copysign function or bit manipulation.)
+
+<li>Error handling
 <p>
-int32_t is used most often which is wrong because of
-implementation defined signed int representation.
+Arithmetics may set various floating point exception flags as a side effect.
+These can be queried and manipulated (fetestexcept, feraiseexcept,..).
 <p>
-In general signed int is not handled carefully
-in the libm code: scalbn even depends on signed int overflow.
+So special care is needed
+when a library function wants to avoid changing
+the floating point status flags.
+eg. if one wants to check for -0 silently then
+<pre>
+if (x == 0.0 &amp;&amp; 1/x &lt; 0) { /* x is a -0 */ }
+</pre>
+is not ok: == raises invalid exception when x is nan,
+and even if we filter nans out 1/x will raise the
+divbyzero flag.
+<p>
+When a library wants to raise a flag deliberately
+but feraiseexcept is not available for some reason,
+then simple arithmetics can be be used just for their
+exception raising side effect
+(eg. 1/0.0 to raise divbyzero), however beaware
+of compiler optimizations (constant folding and dead code elimination,..).
+<p>
+Unfortunately gcc does not always take fp exceptions into
+account: a simple x = 1e300*1e300; may not raise overflow
+exception at runtime, but get optimized into x = +inf.
+see compiler optimizations above.
+<p>
+Another x87 gcc bug related to fp exceptions is that in some cases
+comparision operators (==, &lt;, etc) don't raise invalid
+when an operand is nan
+(eventhough this is required by ieee + c99 annex F).
+(see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52451">gcc bug52451</a>).
+<p>
+The ieee standard defines signaling and quiet nan
+floating-point numbers as well.
+The c99 standard only considers quiet nan, but it allows
+signaling nans to be supported as well.
+Without signaling nans x * 1 is equivalent to x,
+but if signaling nan is supported then the former
+raises an invalid exception.
+This may complicate things further if one wants to write
+portable fp math code.
+<p>
+A further libm design issue is the math_errhandling macro:
+it specifies the way math function errors can be checked
+which is either fp exceptions or the errno variable or both.
+The fp exception approach can be supported with enough care
+but errno is hard to support: certain library functions
+are implemented as a single asm instruction (eg sqrt),
+the only way to set errno is to query the fp exception flags
+and then set the errno variable based on that.
+So eventhough errno may be convenient, in libm it is
+not the right thing to do.
+<p>
+For soft-float targets however errno seems to be the only option
+(which means annex K cannot be fully supported, as it requires
+the support of exception flags).
+The problem is that at context switches the fpu status should
+be saved and restored which is done by the kernel on hard-fp
+architectures when the state is in an fpu status word.
+In case of soft-fp emulation this must be done by the c runtime:
+context switches between threads can be supported with thread local
+storage of the exception state, but signal handlers may do floating-point
+arithmetics which should not alter the fenv state.
+Wrapping signal handlers is not possible/difficult for various
+reasons and the compiler cannot know which functions will be used
+as signal handlers, so the c runtime has no way to guarantee that
+signal handlers do not alter the fenv.
 
 <li>Complex arithmetics
 <p>
@@ -337,6 +498,31 @@ It's not clear how I (or _Complex_I) should be defined
 in complex.h
 (literal of the float imaginary unit is compiler specific,
 in gcc it can be 1.0fi).
+
+<li>Signed int
+<p>
+The freebsd libm code has many inconsistencies
+(naming conventions, 0x1p0 notation vs decimal notation,..),
+one of them is the integer type used for bit manipulations:
+The bits of a double are unpacked into one of
+int, int32_t, uint32_t and u_int32_t
+integer types.
+<p>
+int32_t is used the most often which is not wrong in itself
+but it is used incorrectly in many places.
+<p>
+int is a bit worse because unlike int32_t it is not guaranteed
+to be 32bit two's complement representation. (but of course in
+practice they are the same)
+<p>
+The issues found so far are left shift of negative integers
+(undefined behaviour), right shift of negative integers
+(implementation defined behaviour), signed overflow
+(implementation defined behaviour), unsigned to signed conversion
+(implementation defined behaviour).
+<p>
+It is easy to avoid these issues without performance impact,
+but a bit of care should be taken around bit manipulations.
 </ul>
 
 <h3><a name="implementations" href="#implementations">libm implementations</a></h3>
@@ -392,6 +578,7 @@ various other closed ones: intel libm, hp libm for itanium,..
 <a href="http://gforge.inria.fr/projects/mpcheck/">mpcheck</a>
 <li><a href="http://people.inf.ethz.ch/gonnet/FPAccuracy/Analysis.html">FPAccuracy</a>
 <li><a href="http://www.math.utah.edu/~beebe/software/ieee/">beebe's ieee tests</a>
+<li><a href="http://www.math.utah.edu/pub/elefunt/">elefunt</a>
 <li><a href="http://www.vinc17.org/research/testlibm/index.en.html">testlibm</a>
 <li><a href="http://www.vinc17.org/research/fptest.en.html">fptest</a>
 <li>tests in crlibm