bug list update (posix signal handler)

[www] / libm / index.html

diff --git a/libm/index.html b/libm/index.html
index e177d6c..6eb718c 100644 (file)
--- a/libm/index.html
+++ b/libm/index.html
@@ -35,7 +35,7 @@ cvs -d <a href="http://openbsd.org/anoncvs.html#CVSROOT">$CVSROOT</a> get src/li
 
 <h3><a name="rules" href="#rules">General rules</a></h3>
 <ul>
 
 <h3><a name="rules" href="#rules">General rules</a></h3>
 <ul>

-<li>Assumption about floating-point representation and arithmetics
+<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>float is ieee binary32
 (see <a href="http://port70.net/~nsz/c/c99/n1256.html#F.2">c99 annex F.2</a>):
 <ul>
 <li>float is ieee binary32


@@ -57,7 +57,7 @@ but the last bit may be wrong:
 <pre>
     |error| &lt; 1.5 ulp
 </pre>
 <pre>
     |error| &lt; 1.5 ulp
 </pre>

-should hold.
+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,
 (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,


@@ -69,16 +69,21 @@ 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.
 <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)
+(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
 <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 sideeffect, other
+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)
 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>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


@@ -88,13 +93,13 @@ macros should be protected by #ifdef
 (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.
 (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
+<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
 (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 ((float)0.1 is not ok, that style may lead
-to double rounding issues, but eg. 1.0 or 0.5 can be used instead of 1.0f or
-0.5f)
+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
 <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


@@ -111,6 +116,9 @@ functions they implement.
 <p>
 Binary representation of floating point numbers matter
 because bit hacks are often needed in the math code.
 <p>
 Binary representation of floating point numbers matter
 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 using
 union types:
 <p>
 float and double bit manipulation can be handled in a portable way in c using
 union types:


@@ -173,6 +181,8 @@ on the endianness and write different code for different architectures.
 <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>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:


@@ -187,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


@@ -199,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>


@@ -229,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>,


@@ -260,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:


@@ -278,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


@@ -292,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)


@@ -307,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>


@@ -358,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:


@@ -389,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>


@@ -429,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>