libm: code instead of branch

[www] / libm / index.html

diff --git a/libm/index.html b/libm/index.html
index 296e769..313773e 100644 (file)
--- a/libm/index.html
+++ b/libm/index.html
@@ -7,6 +7,7 @@
 so already existing code is used.
 
 <ul>
 so already existing code is used.
 
 <ul>

+<li><a href="#code">Code</a>

 <li><a href="#sources">Sources</a>
 <li><a href="#issues">Design issues</a>
 <li><a href="#representation">Representation</a>
 <li><a href="#sources">Sources</a>
 <li><a href="#issues">Design issues</a>
 <li><a href="#representation">Representation</a>


@@ -15,6 +16,24 @@ so already existing code is used.
 <li><a href="#tests">libm tests</a>
 </ul>
 
 <li><a href="#tests">libm tests</a>
 </ul>
 

+<h3><a name="code" href="#code">Code</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>
 <p>The math code is mostly from freebsd which in turn
 is based on <a href="http://www.netlib.org/fdlibm/">fdlibm</a>.
 <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,65 +50,32 @@ 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>
 </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/.
-<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/?
-</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.
 
 <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:
 <ul>
 <li>Use STDC pragmas (eventhough gcc does not support them)?
 <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
 </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>
 <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>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

 </ul>
 </ul>
 
 </ul>
 </ul>
 


@@ -104,7 +90,8 @@ in a portable way in c:
 <li>float is ieee binary32
 <li>double is ieee binary64
 </ul>
 <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
 <p>
 The endianness may still vary, but that can be worked
 around by using a union with a single large enough


@@ -123,8 +110,8 @@ various representations:
 (and other non-standard formats are not supported)
 <p>
 ld64 is easy to handle: all long double functions
 (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),
 <p>
 ld80 is the most common (i386, x86_64), it means
 64bit significand with explicit msb (inconsistent with other ieee formats),


@@ -136,7 +123,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
 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>
 
 
 <h3><a name="ugly" href="#ugly">Ugly</a></h3>


@@ -144,7 +133,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>
 
 <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>
 <ul>
 <li>Double rounding:
 <p>


@@ -159,7 +147,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)
 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.
 <p>
 So x = a+b may give different results depending on
 the x87 fpu precision setting.


@@ -168,8 +156,8 @@ but that's the most common one)
 <p>
 (double rounding can happen with float vs double as well)
 <p>
 <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>
 
 <li>Wider exponent range (x87 issue):
 <p>


@@ -190,6 +178,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.
 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)
 So
 <pre>
 (a+b)==(a+b)


@@ -218,8 +209,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
 (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
 <p>
 Unfortunately gcc does not respect the standard
 and even if assingment or cast is used the result


@@ -233,9 +223,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
 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
 <p>
 (Sometimes the excess precision is good
 but it's hard to rely on it as it is optional


@@ -267,7 +257,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>
 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
 <p>
 The workaround is again using named volatile
 variables for constants like


@@ -279,33 +270,86 @@ 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,
 (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.
 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.
 
 <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>
 <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>
 <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 (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 (==, &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
+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>
 
 <li>Complex arithmetics
 <p>


@@ -339,6 +383,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).
 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>
 </ul>
 
 <h3><a name="implementations" href="#implementations">libm implementations</a></h3>