From: nsz <nsz@port70.net>
Date: Sat, 3 Mar 2012 20:50:05 +0000 (+0100)
Subject: libm updates, link to prototype lib
X-Git-Url: http://nsz.repo.hu/git/?a=commitdiff_plain;h=d1ead004efc7440632d5c475f8736203ad37e286;p=www

libm updates, link to prototype lib
---

diff --git a/libm/index.html b/libm/index.html
index 296e769..611a08e 100644
--- a/libm/index.html
+++ b/libm/index.html
@@ -7,6 +7,7 @@
 so already existing code is used.
 
 <ul>
+<li><a href="#prototype">Prototype</a>
 <li><a href="#sources">Sources</a>
 <li><a href="#issues">Design issues</a>
 <li><a href="#representation">Representation</a>
@@ -15,6 +16,23 @@ so already existing code is used.
 <li><a href="#tests">libm tests</a>
 </ul>
 
+<h3><a name="prototype" href="#prototype">Prototype</a></h3>
+<ul>
+<li><a href="/git/?p=libm">browse</a> with gitweb
+<li>over git (recommended)
+<pre>
+git clone git://nsz.repo.hu:45100/repo/epoint
+</pre>
+<li>over http (slower)
+<pre>
+git clone http://nsz.repo.hu/repo/epoint
+</pre>
+<li>over ssh (passwd is anon)
+<pre>
+git clone ssh://anon@nsz.repo.hu:45022/repo/epoint
+</pre>
+</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>.
@@ -31,36 +49,17 @@ cvs -d <a href="http://openbsd.org/anoncvs.html#CVSROOT">$CVSROOT</a> get src/li
 </ul>
 
 <h3><a name="issues" href="#issues">Design issues</a></h3>
-<p>The math code is available but there are still many design questions.
+
 <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/.
+<li>Code organization
+<p>The long double code uses ifdefs and macro hacks to allow
+ld64, ld80 and ld128 implementation in the same file.
 <p>Pending questions:
 <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/?
+<li>What are the correct long double ifdefs?
+<li>Check fpu precision setting at runtime in long double code?
+<li>support c99 tgmath.h?
+<li>TODO: complex
 </ul>
 
 <li>Workarounds
@@ -70,26 +69,24 @@ 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?
+<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 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>The fdlibm code has code style and correctness issues
+which might worth addressing.
 <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>How much c99 is allowed? (inline, compound literals, 0x1p0,..)
+<li>Which one should be used 0f, 0l, 0.0, static const double zero = 0,.. ?
+<li>Can we use a polyeval(coeffs, x) function instead of writing out the poly?
+<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
 </ul>
 </ul>
 
@@ -104,7 +101,8 @@ in a portable way in c:
 <li>float is ieee binary32
 <li>double is ieee binary64
 </ul>
-(and old non-standard representations are not supported)
+(<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
@@ -123,8 +121,8 @@ various representations:
 (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)
+are just wrappers around the corresponding double ones
+(aliasing is non-conformant)
 <p>
 ld80 is the most common (i386, x86_64), it means
 64bit significand with explicit msb (inconsistent with other ieee formats),
@@ -136,7 +134,9 @@ ld128 is rare (sparc64 with software emulation), it means
 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)
+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.
 
 
 <h3><a name="ugly" href="#ugly">Ugly</a></h3>
@@ -144,7 +144,6 @@ m68k and m88k cpus use different ld80 than the intel ones)
 <p>Some notes are from:
 <a href="http://www.vinc17.org/research/extended.en.html">http://www.vinc17.org/research/extended.en.html</a>
 
-
 <ul>
 <li>Double rounding:
 <p>
@@ -159,7 +158,7 @@ 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)
+(freebsd and openbsd use double precision by default)
 <p>
 So x = a+b may give different results depending on
 the x87 fpu precision setting.
@@ -168,8 +167,8 @@ but that's the most common one)
 <p>
 (double rounding can happen with float vs double as well)
 <p>
-C99 annex F prohibits double rounding,
-but that's non-normative.
+<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.
 
 <li>Wider exponent range (x87 issue):
 <p>
@@ -190,6 +189,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)
@@ -218,8 +220,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
@@ -233,9 +234,9 @@ gcc 4.5 fixed it with '-fexcess-precision=standard'
 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
@@ -267,7 +268,8 @@ 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
@@ -279,33 +281,85 @@ and using the '-frounding-math' gcc flag.
 (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:
+but as a workaround 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
+<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>
-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.
+Arithmetics may set various floating point exception flags as a side effect.
+These can be queried and manipulated (fetestexcept, feraiseexcept,..).
 <p>
-int32_t is used most often which is wrong because of
-implementation defined signed int representation.
-<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 deliberatly
+but feraiseexcept is not available for some reason,
+then simple arithmetics can be be used just for its
+exception raising side effect
+(eg. 1/0.0 to raise divbyzero), however beaware
+of compiler optimizations (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 is that comparision operators (==, &lt;, etc)
+don't raise exceptions on 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
+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
+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.
 
 <li>Complex arithmetics
 <p>
@@ -339,6 +393,21 @@ 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 bitmanipulations:
+The bits of a double are unpacked into one of
+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.
+<p>
+In general signed int is not handled carefully
+in the libm code: scalbn even depends on signed int overflow.
 </ul>
 
 <h3><a name="implementations" href="#implementations">libm implementations</a></h3>