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>
<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>.
</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
<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>
<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
(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),
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>
<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 (x87 issue):
+<li>Double rounding:
<p>
If a value rounded twice the result can be different
than rounding just once.
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 fpu setting.
+the x87 fpu precision setting.
(only happens in round to nearest rounding mode,
but that's the most common one)
<p>
-C99 annex F prohibits double rounding,
-but that's non-normative.
+(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.
<li>Wider exponent range (x87 issue):
<p>
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)
(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
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
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
(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
-<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>
-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>
-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 && 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.
+<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 (==, <, 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>
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>
<a href="http://gforge.inria.fr/projects/mpcheck/">mpcheck</a>
<li><a href="http://people.inf.ethz.ch/gonnet/FPAccuracy/Analysis.html">FPAccuracy</a>
<li><a href="http://www.math.utah.edu/~beebe/software/ieee/">beebe's ieee tests</a>
+<li><a href="http://www.math.utah.edu/pub/elefunt/">elefunt</a>
<li><a href="http://www.vinc17.org/research/testlibm/index.en.html">testlibm</a>
<li><a href="http://www.vinc17.org/research/fptest.en.html">fptest</a>
<li>tests in crlibm
projects / www / blobdiff