<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
<pre>
|error| < 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,
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, frexp, ldexp,
-modf, nearbyint, nextafter, nexttoward, rint, round, scalbln, scalbn,
-sqrt, trunc.
+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)
+(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 sideeffect, other
+(FORCE_EVAL is used 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
(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
-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
<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:
</ul>
where m is the significand and se is the sign and exponent.
<p>
-ld128 is rare (eg. sparc64 with software emulation), it means
+hardware ld128 is rare (eg. sparc64 and aarch64 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;};
<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:
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
the x87 fpu precision setting.
(only happens in round to nearest rounding mode,
-but that's the most common one)
+but that's the default)
<p>
(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.
+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>Evaluation precision:
<p>
-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).
+C with FLT_EVAL_METHOD>=0 requires consistent evaluation
+precision, but gcc does not follow this on i386 and may
+store intermediate results and round them to double while keep
+other parts in extended precision.
+(And the precision gcc uses can be inconsistent:
+adding a printf to the code may spill 80bit float registers
+into memory changing the result of a nearby calculation).
So
<pre>
(a+b)==(a+b)
</pre>
-may be false when the two sides
+may be false with gcc and FLT_EVAL_METHOD!=0, 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.
-It goes away if the highest precision (long double) is used
-everywhere, but that can have a huge penalty).
+If the highest precision (long double) is used
+everywhere then there is no such issue, 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>,
both sides are guaranteed to be in double precision
when the comparision is done.
<p>
-(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 still be false)
-<p>
Unfortunately gcc does not respect the standard
and even if assingment or cast is used the result
may be kept in higher precision
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
+predictable)
<p>
The workaround for older gcc is to force the
compiler to store the intermediate results:
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
-<p>
-There is a way to check for it though using
-FLT_EVAL_METHOD, float_t and double_t.
-But it's probably easier to unconditionally
-use higher precision variables when that's
-what we want and don't depend on the implicit
-excess precision).
+then float_t and double_t may be used explicitly)
+
+<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 sets x87 precision to double precision by default
+where gcc is hosted and gcc cannot get the constants right for the target platform,
+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>
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)
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>
<pre>
if (x == 0.0 && 1/x < 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.
+is not ok: 1/x raises 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 (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>
-Another x87 gcc bug related to fp exceptions is that
-comparision operators (==, <, etc) don't raise invalid
-when an operand is nan
+Another x87 gcc bug related to fp exceptions is that in some cases
+relation operators (<, 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.
-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:
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.
+<p>
+For soft-float targets exceptions are problematic in c99.
+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.
+C11 fixed this by changing signal handler semantics: fenv is in unspecified state.
<li>Complex arithmetics
<p>
<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
-int32_t, uint32_t and u_int32_t
+int, int32_t, uint32_t and u_int32_t
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>
-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>
projects / www / blobdiff