1 <html><head><title>libm</title></head><body>
4 <p>This page is about libm for the
5 <a href="http://www.musl-libc.org/">musl</a> libc.
8 <li><a href="#sources">Sources</a>
9 <li><a href="#rules">General rules</a>
10 <li><a href="#representation">Representation</a>
11 <li><a href="#ugly">Ugly</a>
12 <li><a href="#implementations">libm implementations</a>
13 <li><a href="#tests">libm tests</a>
16 <h3><a name="sources" href="#sources">Sources</a></h3>
17 <p>Writing math code from scratch is a huge work so already existing code is
18 used. Several math functions are taken from the
19 <a href="http://freebsd.org/">freebsd</a> libm and a few from the
20 <a href="http://openbsd.org/">openbsd</a> libm implementations.
21 Both of them are based on <a href="http://www.netlib.org/fdlibm/">fdlibm</a>.
22 The freebsd libm seems to be the most well maintained and most correct version
26 <li>freebsd /lib/msun (<a href="http://svnweb.FreeBSD.org/base/head/lib/msun/">browse</a>)
28 svn checkout svn://svn.freebsd.org/base/head/lib/msun
30 <li>openbsd /src/lib/libm (<a href="http://www.openbsd.org/cgi-bin/cvsweb/src/lib/libm/">browse</a>)
32 cvs -d <a href="http://openbsd.org/anoncvs.html#CVSROOT">$CVSROOT</a> get src/lib/libm
36 <h3><a name="rules" href="#rules">General rules</a></h3>
38 <li>Assumptions about floating-point representation and arithmetics
39 (see <a href="http://port70.net/~nsz/c/c99/n1256.html#F.2">c99 annex F.2</a>):
41 <li>float is ieee binary32
42 <li>double is ieee binary64
43 <li>long double is either ieee binary64 or little-endian 80bit extended precision (x87 fpu)
45 (other long double representations may be supported in the future, until then
46 long double math functions will be missing on non-supported platforms)
47 <li>On soft-float architectures fenv should not be expected to work according to
48 the c and ieee standards (ie. rounding modes and exceptions are not supported,
49 the fenv functions are dummy ones)
50 <li>Floating-point exception flags should be properly set in math functions
51 according to c99 annex F, but without using fenv.h
52 (eg. overflow flag can be raised by 0x1p900*0x1p900, because this works even
54 <li>Most functions need to be precise only in nearest rounding mode.
55 <li>Returned floating-point values should be correctly rounded in most cases,
56 but the last bit may be wrong:
60 should hold for most functions.
61 (error is the difference between the exact result and the calculated
63 (in theory correct rounding can be achieved but with big implementation cost,
64 see <a href="http://lipforge.ens-lyon.fr/www/crlibm/">crlibm</a>)
65 <li>At least the following functions must be correctly rounded:
66 ceil, copysign, fabs, fdim, floor, fma, fmax, fmin, frexp, ldexp, logb,
67 modf, nearbyint, nextafter, nexttoward, rint, round, scalbln, scalbn,
69 <li>Mathematical properties of functions should be as expected
70 (monotonicity, range, symmetries).
71 <li>If the FPU precision is altered then nothing is guaranteed to work.
72 (ie. when long double does not have full 80bit precision on i386 then things
73 may break, this also means that compiler must spill fpu registers at 80bits
75 <li>Signaling NaN is not supported
76 <li>Quiet NaN is supported but all NaNs are treated equally without special
77 attention to the internal representation of a NaN
78 (eg. the sign of NaN may not be preserved).
79 <li>Most gcc bug workarounds should be removed from the code
80 (STRICT_ASSIGN macro is used when excessive precision is harmful and
81 FORCE_EVAL when expressions must be evaluated for their side-effect, other
82 usage of volatile is not justified, hacks around long double constants are
83 not justified eventhough gcc can miscompile those with non-default FPU setting)
84 <li>Whenever fenv is accessed the FENV_ACCESS pragma of c99 should be used
85 (eventhough gcc does not yet support it), and all usage of optional FE_
86 macros should be protected by #ifdef
87 <li>For bit manipulation of floating-point values an union should be used
88 (eg. union {float f; uint32_t i;})
89 <li>uint32_t and uint64_t should be used for bit manipulations.
90 (eg signed int must not be used in bit shifts etc when it might invoke
91 undefined or implementation defined behaviour).
92 <li>POSIX namespace rules must be respected.
93 <li>c99 hexfloat syntax (0x1.0p0) should be used when it makes the code
94 clearer, but not in public header files
95 (those should be c++ and ansi c compatible)
96 <li>The 'f' suffix should be used for single precision values (0.1f) when the
97 value cannot be exactly represented or the type of the arithmetics is important
98 ((float)0.1 is not ok, that style may lead to double rounding issues, but eg.
99 1.0 or 0.5 may be used instead of 1.0f or 0.5f in some cases)
100 <li>Prefer classification macros (eg. isnan) over inplace bit hacks.
101 <li>For every math function there should be a c implementation.
102 (a notable exception now is sqrtl, since most fpu has instruction for it
103 and on soft-float architectures long double == double)
104 <li>The c implementation of a long double function should use ifdefs with the
105 LDBL_MANT_DIG etc constants from float.h for architecture specific
107 <li>In the musl source tree math.h functions go to src/math, complex.h functions
108 to src/complex and fenv.h functions to src/fenv. And files are named after the
109 functions they implement.
112 <h3><a name="representation" href="#representation">Representation</a></h3>
114 Binary representation of floating point numbers matter
115 because bit hacks are often needed in the math code.
116 (in particular bit hacks are used instead of relational operations for nan
117 and sign checks becuase relational operators raise invalid fp exception on nan
118 and they treat -0.0 and +0.0 equally and more often than not these are not desired)
120 float and double bit manipulation can be handled in a portable way in c using
123 <li>union {float f; uint32_t i;};
124 <li>union {double f; uint64_t i;};
126 (assuming the bits in the object representation of 32bit and 64bit unsigned ints
127 map to the floating-point representation according to ieee-754, this is not
128 always the case, eg. old
129 <a href="http://wiki.debian.org/ArmEabiPort#ARM_floating_points">arm floating-point accelerator</a>
130 (FPA) used mixed endian double representation, but musl does not support the old
133 long double bit manipulation is harder as there are various representations
134 and some of them don't map to any unsigned integer type:
136 <li>ld64: long double is the same as double (ieee binary64)
137 <li>ld80: 80bit extended precision format [<a href="https://en.wikipedia.org/wiki/Extended_precision">wikipedia</a>]
138 <li>ld128: quadruple precision, ieee binary128 [<a href="https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format">wikipedia</a>]
140 (and other non-standard formats are not supported)
142 In case of ld64 the bit manipulation is the same as with double
143 and all long double math functions can be just wrappers around the
144 corresponding double ones.
145 (using symbol aliasing on the linker level is non-conformant
146 since functions would not have unique address then)
148 ld80 is the most common long double on linux (i386 and x86_64 abi),
149 it means 64bit significand with explicit msb
150 (inconsistent with other ieee formats), 15bit exp, 1 sign bit.
151 The m68k (and m88k) architecture uses the same format, but different endianness:
153 <li>union {long double f; struct{uint64_t m; uint16_t se; uint16_t pad;} i;}; // x86
154 <li>union {long double f; struct{uint16_t se; uint16_t pad; uint64_t m;} i;}; // m68k
156 where m is the significand and se is the sign and exponent.
158 ld128 is rare (eg. sparc64 with software emulation), it means
159 113bit significand with implicit msb, 15bit exp, 1 sign bit:
161 <li>union {long double f; struct{uint16_t se; uint16_t hi; uint32_t mid; uint64_t lo;} i;};
164 There are other non-conformant long double types: eg. the old SVR4 abi for ppc
165 uses 128 bit long doubles, but it's software emulated and traditionally
167 <a href="https://en.wikipedia.org/wiki/Quadruple_precision#Double-double_arithmetic">two doubles</a>
168 (also called ibm long double as this is what ibm aix used on ppc).
169 The ibm s390 supports the ieee 754-2008 compliant binary128 floating-point
170 format, but previous ibm machines (S/370, S/360) used slightly different
173 This variation shows the difficulty to consistently handle
174 long double: the solution is to use ifdefs based on float.h and
175 on the endianness and write different code for different architectures.
177 <h3><a name="ugly" href="#ugly">Ugly</a></h3>
178 <p>The ugly parts of libm hacking.
179 <p>Some notes are from:
180 <a href="http://www.vinc17.org/research/extended.en.html">http://www.vinc17.org/research/extended.en.html</a>
185 If a value rounded twice the result can be different
186 than rounding just once.
188 The C language allows arithmetic to be evaluated in
189 higher precision than the operands have. If we
190 use x87 fpu in extended precision mode it will round
191 the results twice: round to 80bit when calculating
192 and then round to 64bit when storing it, this can
193 give different result than a single 64bit rounding.
194 (on x86-linux the default fpu setting is to round the
195 results in extended precision, this only affects x87 instructions, not see2 etc)
196 (freebsd and openbsd use double precision by default)
198 So x = a+b may give different results depending on
199 the x87 fpu precision setting.
200 (only happens in round to nearest rounding mode,
201 but that's the most common one)
203 (double rounding can happen with float vs double as well)
205 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#F.7.3">C99 annex F</a>
206 prohibits double rounding, but that's non-normative.
208 <li>Wider exponent range (x87 issue):
210 Even if the fpu is set to double precision
211 (which is not) the x87 registers use wider exponent
212 range (mant:exp is 53:15 instead of 53:11 bits)
213 so underflows (subnormals) may not be treated
214 as expected. Rounding to double only occurs
215 when a value is stored into memory.
217 Actually this beahviour is allowed by the ieee 754
218 standard, but it can cause problems (see
219 <a href="http://www.exploringbinary.com/php-hangs-on-numeric-value-2-2250738585072011e-308/">infamous php bug</a>)
221 <li>Evaluation precision:
223 C does not require consistent evaluation
224 precision: the compiler may store intermediate
225 results and round them to double while keep
226 other parts in higher precision.
227 (And the precision the compiler choses can
228 be inconsistent: adding a printf to the code
229 may change the result of a nearby calculation).
234 may be false when the two sides
235 are kept in different precision.
236 (This is not an x87 specific problem, it matters whenever there
237 is a higher precision fp type than the currently used one.
238 It goes away if the highest precision (long double) is used
239 everywhere, but that can have a huge penalty).
241 C99 has a way to control this (see
242 <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>,
243 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#6.3.1.5">6.3.1.5</a> and
244 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#6.3.1.8">6.3.1.8</a>):
245 when a result is stored in a variable or a
246 type cast is used, then it is guaranteed that
247 the precision is appropriate to that type.
250 (double)(a+b)==(double)(a+b)
252 both sides are guaranteed to be in double precision
253 when the comparision is done.
255 (This still does not solve the x87 double rounding
256 issue though: eg if the left side is evaluated with
257 sse2 and the right side with x87 extended precision setting
258 and double rounding then the result may still be false)
260 Unfortunately gcc does not respect the standard
261 and even if assingment or cast is used the result
262 may be kept in higher precision
263 (infamous <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323">gcc bug323</a>
264 also see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36578">gcc bug36578</a>).
265 gcc 4.5 fixed it with '-fexcess-precision=standard'
266 (it is enabled by '-std=c99', but the default is
267 '-fexcess-precision=fast')
269 The workaround for older gcc is to force the
270 compiler to store the intermediate results:
271 by using volatile double temporary variables
272 or by '-ffloat-store' (which affects all
273 intermediate results, but is not guaranteed
274 by the gcc docs to always work).
276 (Sometimes the excess precision is good
277 but it's hard to rely on it as it is optional
279 There is a way to check for it though using
280 FLT_EVAL_METHOD, float_t and double_t.
281 But it's probably easier to unconditionally
282 use higher precision variables when that's
283 what we want and don't depend on the implicit
286 <li>Compiler optimizations:
288 Runtime and compile time semantics may be different
289 gcc often does unwanted or even invalid compile
291 (constant folding where floating point exception
292 flags should be raised at runtime, or the result
293 depends on runtime floating point environment,
294 or constant expressions are just evaluated with
295 different precision than at runtime).
297 C99 actually allows most of these optimizations
298 but they can be turned off with STDC pragmas (see
299 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#6.10.6">6.10.6</a>).
300 Unfortunately <a href="http://gcc.gnu.org/c99status.html">gcc does not support these pragmas</a>.
302 FENV_ACCESS ON tells the compiler that the code wants
303 to access the floating point environment (eg. set different rounding mode)
304 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34678">gcc bug34678</a>).
306 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845">gcc bug37845</a> for FP_CONTRACT pragma
307 which is relevant when the architecture has fma instruction, x86_64, i386, arm do not have it).
309 The workaround is again using named volatile
310 variables for constants like
312 static const volatile two52 = 0x1p52;
314 and using the '-frounding-math' gcc flag.
316 According to the freebsd libm code gcc truncates long double
317 const literals on i386.
318 I assume this happens because freebsd uses 64bit long doubles by default
319 (double precision) and gcc incorrectly uses the precision setting of the
320 host platform instead of the target one, but i did not observe this on linux.
321 (as a workaround sometimes double-double arithmetics was used
322 to initialize long doubles on i386, but most of these should be
323 fixed in musl's math code now)
327 The two representations are sufficiently different
328 that treating them together is awkward.
329 (Especially the explicit msb bit in ld80 can cause
330 different behaviour).
332 In the freebsd libm code a few architecture specific
333 macros and a union handle these issues, but the
334 result is often less clear than treating ld80
335 and ld128 separately.
338 <p>Signed zeros can be tricky. They cannot be checked
339 using the usual comparision operators (+0.0 == -0.0 is true),
340 but they give different results in some cases
341 (1/-0.0 == -inf) and can make similarly looking
342 expressions different eg.
343 (x + 0) is not the same as x, the former is +0.0 when x is -0.0
345 (To check for -0, the signbit macro can be used
346 or the copysign function or bit manipulation.)
350 Arithmetics may set various floating point exception flags as a side effect.
351 These can be queried and manipulated (fetestexcept, feraiseexcept,..).
353 So special care is needed
354 when a library function wants to avoid changing
355 the floating point status flags.
356 eg. if one wants to check for -0 silently then
358 if (x == 0.0 && 1/x < 0) { /* x is a -0 */ }
360 is not ok: == raises invalid exception when x is nan,
361 and even if we filter nans out 1/x will raise the
364 When a library wants to raise a flag deliberately
365 but feraiseexcept is not available for some reason,
366 then simple arithmetics can be be used just for their
367 exception raising side effect
368 (eg. 1/0.0 to raise divbyzero), however beaware
369 of compiler optimizations (constant folding and dead code elimination,..).
371 Unfortunately gcc does not always take fp exceptions into
372 account: a simple x = 1e300*1e300; may not raise overflow
373 exception at runtime, but get optimized into x = +inf.
374 see compiler optimizations above.
376 Another x87 gcc bug related to fp exceptions is that in some cases
377 comparision operators (==, <, etc) don't raise invalid
378 when an operand is nan
379 (eventhough this is required by ieee + c99 annex F).
380 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52451">gcc bug52451</a>).
382 The ieee standard defines signaling and quiet nan
383 floating-point numbers as well.
384 The c99 standard only considers quiet nan, but it allows
385 signaling nans to be supported as well.
386 Without signaling nans x * 1 is equivalent to x,
387 but if signaling nan is supported then the former
388 raises an invalid exception.
389 This may complicate things further if one wants to write
390 portable fp math code.
392 A further libm design issue is the math_errhandling macro:
393 it specifies the way math function errors can be checked
394 which is either fp exceptions or the errno variable or both.
395 The fp exception approach can be supported with enough care
396 but errno is hard to support: certain library functions
397 are implemented as a single asm instruction (eg sqrt),
398 the only way to set errno is to query the fp exception flags
399 and then set the errno variable based on that.
400 So eventhough errno may be convenient, in libm it is
401 not the right thing to do.
403 For soft-float targets however errno seems to be the only option
404 (which means annex K cannot be fully supported, as it requires
405 the support of exception flags).
406 The problem is that at context switches the fpu status should
407 be saved and restored which is done by the kernel on hard-fp
408 architectures when the state is in an fpu status word.
409 In case of soft-fp emulation this must be done by the c runtime:
410 context switches between threads can be supported with thread local
411 storage of the exception state, but signal handlers may do floating-point
412 arithmetics which should not alter the fenv state.
413 Wrapping signal handlers is not possible/difficult for various
414 reasons and the compiler cannot know which functions will be used
415 as signal handlers, so the c runtime has no way to guarantee that
416 signal handlers do not alter the fenv.
418 <li>Complex arithmetics
426 (x+0*I)*(0+I) = 0*x + x*I
428 So if x=inf then the result is nan+inf*I instead of inf*I
429 (an fmul instruction is generated for 0*x)
431 (There were various complex constant folding issues as well in gcc).
433 So a+b*I cannot be used to create complex numbers,
434 instead some double[2] manipulation should be used to
435 make sure there is no fmul in the generated code
436 (freebsd uses a static inline cpack function)
437 (It's not clear which is better: using union hacks
438 everywhere in the complex code or creal, cimag, cpack)
440 Complex code is hard to make portable across compilers
441 as there are no good compiler support, it is rarely used
442 and there are various options available in the standard:
443 complex support is optional, imaginary type
446 It's not clear how I (or _Complex_I) should be defined
448 (literal of the float imaginary unit is compiler specific,
449 in gcc it can be 1.0fi).
453 The freebsd libm code has many inconsistencies
454 (naming conventions, 0x1p0 notation vs decimal notation,..),
455 one of them is the integer type used for bit manipulations:
456 The bits of a double are unpacked into one of
457 int, int32_t, uint32_t and u_int32_t
460 int32_t is used the most often which is not wrong in itself
461 but it is used incorrectly in many places.
463 int is a bit worse because unlike int32_t it is not guaranteed
464 to be 32bit two's complement representation. (but of course in
465 practice they are the same)
467 The issues found so far are left shift of negative integers
468 (undefined behaviour), right shift of negative integers
469 (implementation defined behaviour), signed overflow
470 (implementation defined behaviour), unsigned to signed conversion
471 (implementation defined behaviour).
473 It is easy to avoid these issues without performance impact,
474 but a bit of care should be taken around bit manipulations.
477 <h3><a name="implementations" href="#implementations">libm implementations</a></h3>
480 unix v7 libm by Robert Morris (rhm) around 1979<br>
481 see: "Computer Approximations" by Hart & Cheney, 1968<br>
484 cephes by Stephen L. Moshier, 1984-1992<br>
485 <a href="http://www.netlib.org/cephes/">http://www.netlib.org/cephes/</a><br>
486 see: "Methods and Programs for Mathematical Functions" by Moshier 1989
488 fdlibm by K. C. Ng around 1993, updated in 1995?<br>
489 (1993: copyright notice: "Developed at SunPro ..")<br>
490 (1995: copyright notice: "Developed at SunSoft ..")<br>
491 eg see: "Argument Reduction for Huge Arguments: Good to the Last Bit" by Ng 1992<br>
492 <a href="http://www.netlib.org/fdlibm">http://www.netlib.org/fdlibm</a><br>
493 used in: netbsd, openbsd, freebsd, bionic, musl
495 libultim by Abraham Ziv around 2001 (lgpl)<br>
496 (also known as the ibm accurate portable math lib)<br>
497 see: "Fast evaluation of elementary mathematical functions with correctly rounded last bit" by Ziv 1991<br>
498 <a href="http://web.archive.org/web/20050207110205/http://oss.software.ibm.com/mathlib/">http://web.archive.org/web/20050207110205/http://oss.software.ibm.com/mathlib/</a><br>
501 libmcr by K. C. Ng, 2004<br>
502 (sun's correctly rounded libm)
504 mathcw by Nelson H. F. Beebe, 2010?<br>
505 <a href="http://ftp.math.utah.edu/pub/mathcw/">http://ftp.math.utah.edu/pub/mathcw/</a><br>
506 (no sources available?)
508 sleef by Naoki Shibata, 2010<br>
509 (simd lib for evaluating elementary functions)<br>
510 <a href="http://shibatch.sourceforge.net/">http://shibatch.sourceforge.net/</a>
512 crlibm by ens-lyon, 2004-2010<br>
513 <a href="http://lipforge.ens-lyon.fr/www/crlibm/">http://lipforge.ens-lyon.fr/www/crlibm/</a><br>
514 see: <a href="http://lipforge.ens-lyon.fr/frs/?group_id=8&release_id=123">crlibm.pdf</a><br>
515 see: "Elementary Functions" by Jean-Michel Muller 2005<br>
516 see: "Handbook of Floating-Point Arithmetic" by (many authors from Lyon) 2009<br>
518 various other closed ones: intel libm, hp libm for itanium,..
521 <h3><a name="tests" href="#tests">libm tests</a></h3>
523 <li><a href="http://www.netlib.org/fp/ucbtest.tgz">ucbtest.tgz</a>
524 <li><a href="http://www.jhauser.us/arithmetic/TestFloat.html">TestFloat</a>
525 <li><a href="http://cant.ua.ac.be/old/ieeecc754.html">ieeecc754</a>
526 <li><a href="http://www.loria.fr/~zimmerma/mpcheck/">mpcheck</a>,
527 <a href="http://gforge.inria.fr/projects/mpcheck/">mpcheck</a>
528 <li><a href="http://people.inf.ethz.ch/gonnet/FPAccuracy/Analysis.html">FPAccuracy</a>
529 <li><a href="http://www.math.utah.edu/~beebe/software/ieee/">beebe's ieee tests</a>
530 <li><a href="http://www.math.utah.edu/pub/elefunt/">elefunt</a>
531 <li><a href="http://www.vinc17.org/research/testlibm/index.en.html">testlibm</a>
532 <li><a href="http://www.vinc17.org/research/fptest.en.html">fptest</a>
536 multiprecision libs (useful for tests)
538 <li><a href="http://mpfr.org/">mpfr</a>
539 <li><a href="http://www.multiprecision.org/index.php?prog=mpc">mpc</a>
540 <li><a href="http://pari.math.u-bordeaux.fr/">pari</a>
541 <li><a href="http://www.apfloat.org/apfloat/">apfloat</a>
542 <li><a href="http://code.google.com/p/fastfunlib/">fastfunlib</a>
543 <li>scs_lib in crlibm
548 <li>ieee standard: http://754r.ucbtest.org/
549 <li>extended precision issues: http://www.vinc17.org/research/extended.en.html
550 <li>correctly rounded mult: http://perso.ens-lyon.fr/jean-michel.muller/MultConstant.html
551 <li>table based sincos: http://perso.ens-lyon.fr/damien.stehle/IMPROVEDGAL.html
552 <li>finding worst cases: http://perso.ens-lyon.fr/damien.stehle/WCLR.html
553 <li>finding worst cases: http://perso.ens-lyon.fr/damien.stehle/DECIMALEXP.html
554 <li>finding worst cases: http://www.loria.fr/equipes/spaces/slz.en.html
555 <li>finding worst cases: http://perso.ens-lyon.fr/jean-michel.muller/Intro-to-TMD.htm
556 <li>generating libm functions: http://lipforge.ens-lyon.fr/www/metalibm/
557 <li>fast conversion to fixedpoint: http://stereopsis.com/sree/fpu2006.html
558 <li>double-double, quad-double arithmetics: http://crd-legacy.lbl.gov/~dhbailey/mpdist/
559 <li>papers by kahan: http://www.cs.berkeley.edu/~wkahan/
560 <li>fp paper by goldberg: http://download.oracle.com/docs/cd/E19957-01/806-3568/ncg_goldberg.html
561 <li>math functions in general: http://dlmf.nist.gov/
564 <p><small>(<a href="/git/?p=www">page history</a>, <a href="/">home</a>)</small>
projects / www / blob