libm: fmod, remquo, remainder are cr
[www] / libm / index.html
index e65d7f1..dce0a91 100644 (file)
@@ -1,42 +1,26 @@
 <html><head><title>libm</title></head><body>
 <h2>libm</h2>
 
 <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>
 
 <ul>
-<li><a href="#branch">Branch</a>
 <li><a href="#sources">Sources</a>
 <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>
 <li><a href="#tests">libm tests</a>
 </ul>
 
 <li><a href="#representation">Representation</a>
 <li><a href="#ugly">Ugly</a>
 <li><a href="#implementations">libm implementations</a>
 <li><a href="#tests">libm tests</a>
 </ul>
 
-<h3><a name="branch" href="#branch">Branch</a></h3>
-<p>musl branch for math hacks:
-<ul>
-<li><a href="/git/?p=musl">browse</a> with gitweb
-<li>over git (recommended)
-<pre>
-git clone git://nsz.repo.hu:45100/repo/musl
-</pre>
-<li>over http (slower)
-<pre>
-git clone http://nsz.repo.hu/repo/musl
-</pre>
-<li>over ssh (passwd is anon)
-<pre>
-git clone ssh://anon@nsz.repo.hu:45022/repo/musl
-</pre>
-</ul>
-
 <h3><a name="sources" href="#sources">Sources</a></h3>
 <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>)
 <p>sources:
 <ul>
 <li>freebsd /lib/msun (<a href="http://svnweb.FreeBSD.org/base/head/lib/msun/">browse</a>)
@@ -49,59 +33,108 @@ cvs -d <a href="http://openbsd.org/anoncvs.html#CVSROOT">$CVSROOT</a> get src/li
 </pre>
 </ul>
 
 </pre>
 </ul>
 
-<h3><a name="issues" href="#issues">Design issues</a></h3>
-
-<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>In long double code fpu is assumed to be set to extended precision
-<p>Pending questions:
+<h3><a name="rules" href="#rules">General rules</a></h3>
 <ul>
 <ul>
-<li>Use STDC pragmas (eventhough gcc does not support them)?
-<li>Where should we use volatile (to avoid evaluation precision and const folding issues)?
-<li>What tricks to use for proper exception raising? (huge*huge, etc)
-<li>Should signaling nans be considered?
-<li>Which compiler flags should be used (-frounding-math, -ffloat-store)?
-<li>TODO: inline/macro optimization for small functions (isnan, isinf, signbit, creal, cimag,..)
-</ul>
-
-<li>Code cleanups
-<p>The fdlibm code has code style and correctness issues
-which might worth addressing.
-<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>
 <ul>
-<li>Can we use a polyeval(coeffs, x) function instead of writing out the poly?
-<li>TODO: prefer 0x1p0 format over decimal one
-<li>TODO: use 1.0f instead of (float)1 and const float one = 1;
-<li>TODO: use uint32_t instead of int32_t (to avoid signed int overflow and right shifts)
-<li>TODO: prefer isnan, signbit,.. over inplace bithacks
+<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>
 </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
 </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>
 <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>
 <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>
 </ul>
-(<a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#F.2">C99 annex F</a>
-makes these mandatory)
-<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>]
 <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>]
@@ -109,29 +142,47 @@ various representations:
 </ul>
 (and other non-standard formats are not supported)
 <p>
 </ul>
 (and other non-standard formats are not supported)
 <p>
-ld64 is easy to handle: all long double functions
-are just wrappers around the corresponding double ones
-(aliasing is non-conformant)
+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>
 <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>
 <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>
 <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.)
-So there can be further variations in the binary representation
-than just ld80 and ld128.
-
+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>
 
 <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:
 
 <ul>
 <li>Double rounding:
@@ -146,7 +197,7 @@ 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
 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)
+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
 (freebsd and openbsd use double precision by default)
 <p>
 So x = a+b may give different results depending on
@@ -158,6 +209,18 @@ but that's the most common one)
 <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>
 <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>
 
 <li>Wider exponent range (x87 issue):
 <p>
@@ -188,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
 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).
 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>,
 <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>,
@@ -219,6 +285,11 @@ 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')
 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:
 <p>
 The workaround for older gcc is to force the
 compiler to store the intermediate results:
@@ -237,6 +308,40 @@ use higher precision variables when that's
 what we want and don't depend on the implicit
 excess precision).
 
 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
 <li>Compiler optimizations:
 <p>
 Runtime and compile time semantics may be different
@@ -251,7 +356,8 @@ 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>).
 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</a>.
+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)
 <p>
 FENV_ACCESS ON tells the compiler that the code wants
 to access the floating point environment (eg. set different rounding mode)
@@ -266,12 +372,6 @@ variables for constants like
 static const volatile two52 = 0x1p52;
 </pre>
 and using the '-frounding-math' gcc flag.
 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 double-double arithmetics is used:
-initializing the long double constant from two doubles)
 
 <li>ld80 vs ld128
 <p>
 
 <li>ld80 vs ld128
 <p>
@@ -317,27 +417,27 @@ 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
 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 (dead code elimination,..).
+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>
 <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
+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>
 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 quite nan
+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.
 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 complicates things further if one wants to write
+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:
 portable fp math code.
 <p>
 A further libm design issue is the math_errhandling macro:
@@ -348,8 +448,23 @@ 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.
 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
+So eventhough errno may be convenient, in libm it is
 not the right thing to do.
 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>
 
 <li>Complex arithmetics
 <p>
@@ -388,16 +503,26 @@ in gcc it can be 1.0fi).
 <p>
 The freebsd libm code has many inconsistencies
 (naming conventions, 0x1p0 notation vs decimal notation,..),
 <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:
+one of them is the integer type used for bit manipulations:
 The bits of a double are unpacked into one of
 The bits of a double are unpacked into one of
-int32_t, uint32_t and u_int32_t
+int, int32_t, uint32_t and u_int32_t
 integer types.
 <p>
 integer types.
 <p>
-int32_t is used most often which is wrong because of
-implementation defined signed int representation.
+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>
 <p>
-In general signed int is not handled carefully
-in the libm code: scalbn even depends on signed int overflow.
+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>
 </ul>
 
 <h3><a name="implementations" href="#implementations">libm implementations</a></h3>