nsz Git - www/blob - libm/index.html

   1 <html><head><title>libm</title></head><body>
   2 <h2>libm</h2>
   3
   4 <p>This page is about designing libm for the
   5 <a href="http://www.etalabs.net/musl/">musl</a> libc.
   6 <p>Writing the math code from scratch is a huge work
   7 so already existing code is used.
   8
   9 <ul>
  10 <li><a href="#prototype">Prototype</a>
  11 <li><a href="#sources">Sources</a>
  12 <li><a href="#issues">Design issues</a>
  13 <li><a href="#representation">Representation</a>
  14 <li><a href="#ugly">Ugly</a>
  15 <li><a href="#implementations">libm implementations</a>
  16 <li><a href="#tests">libm tests</a>
  17 </ul>
  18
  19 <h3><a name="prototype" href="#prototype">Prototype</a></h3>
  20 <ul>
  21 <li><a href="/git/?p=libm">browse</a> with gitweb
  22 <li>over git (recommended)
  23 <pre>
  24 git clone git://nsz.repo.hu:45100/repo/libm
  25 </pre>
  26 <li>over http (slower)
  27 <pre>
  28 git clone http://nsz.repo.hu/repo/libm
  29 </pre>
  30 <li>over ssh (passwd is anon)
  31 <pre>
  32 git clone ssh://anon@nsz.repo.hu:45022/repo/libm
  33 </pre>
  34 </ul>
  35
  36 <h3><a name="sources" href="#sources">Sources</a></h3>
  37 <p>The math code is mostly from freebsd which in turn
  38 is based on <a href="http://www.netlib.org/fdlibm/">fdlibm</a>.
  39 <p>sources:
  40 <ul>
  41 <li>freebsd /lib/msun (<a href="http://svnweb.FreeBSD.org/base/head/lib/msun/">browse</a>)
  42 <pre>
  43 svn checkout svn://svn.freebsd.org/base/head/lib/msun
  44 </pre>
  45 <li>openbsd /src/lib/libm (<a href="http://www.openbsd.org/cgi-bin/cvsweb/src/lib/libm/">browse</a>)
  46 <pre>
  47 cvs -d <a href="http://openbsd.org/anoncvs.html#CVSROOT">$CVSROOT</a> get src/lib/libm
  48 </pre>
  49 </ul>
  50
  51 <h3><a name="issues" href="#issues">Design issues</a></h3>
  52
  53 <ul>
  54 <li>Code organization
  55 <p>The long double code uses ifdefs and macro hacks to allow
  56 ld64, ld80 and ld128 implementation in the same file.
  57 <p>Pending questions:
  58 <ul>
  59 <li>What are the correct long double ifdefs?
  60 <li>Check fpu precision setting at runtime in long double code?
  61 <li>support c99 tgmath.h?
  62 <li>TODO: complex
  63 </ul>
  64
  65 <li>Workarounds
  66 <p>The bsd libm code has many workarounds for various
  67 compiler issues. It's not clear what's the best way
  68 to handle the uglyness.
  69 <p>Pending questions:
  70 <ul>
  71 <li>Use STDC pragmas (eventhough gcc does not support them)?
  72 <li>Where should we use volatile (to avoid evaluation precision and const folding issues)?
  73 <li>What tricks to use for proper exception raising? (huge*huge, etc)
  74 <li>Should signaling nans be considered?
  75 <li>Which compiler flags should be used (-frounding-math, -ffloat-store)?
  76 <li>TODO: inline/macro optimization for small functions (isnan, isinf, signbit, creal, cimag,..)
  77 </ul>
  78
  79 <li>Code cleanups
  80 <p>The fdlibm code has code style and correctness issues
  81 which might worth addressing.
  82 <p>Pending questions:
  83 <ul>
  84 <li>Should 0x1p0 float format be preferred over decimal format?
  85 <li>How much c99 is allowed? (inline, compound literals, 0x1p0,..)
  86 <li>Which one should be used 0f, 0l, 0.0, static const double zero = 0,.. ?
  87 <li>Can we use a polyeval(coeffs, x) function instead of writing out the poly?
  88 <li>TODO: use uint32_t instead of int32_t (to avoid signed int overflow and right shifts)
  89 <li>TODO: prefer isnan, signbit,.. over inplace bithacks
  90 </ul>
  91 </ul>
  92
  93 <h3><a name="representation" href="#representation">Representation</a></h3>
  94 <p>
  95 Binary representation of floating point numbers matter
  96 because lots of bithacks are needed in the math code.
  97 <p>
  98 float and double bit manipulation can be handled
  99 in a portable way in c:
 100 <ul>
 101 <li>float is ieee binary32
 102 <li>double is ieee binary64
 103 </ul>
 104 (<a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#F.2">C99 annex F</a>
 105 makes these mandatory)
 106 <p>
 107 The endianness may still vary, but that can be worked
 108 around by using a union with a single large enough
 109 unsigned int. (The only exception is probably arm/oabi fpa
 110 where the word order of a double is not consistent with the
 111 endianness [<a href="http://wiki.debian.org/ArmEabiPort#ARM_floating_points">debian wiki on arm</a>],
 112 but we probably won't support that)
 113 <p>
 114 long double bit manipulation is harder as there are
 115 various representations:
 116 <ul>
 117 <li>ld64: long double is the same as double (ieee binary64)
 118 <li>ld80: 80bit extended precision format [<a href="https://en.wikipedia.org/wiki/Extended_precision">wikipedia</a>]
 119 <li>ld128: quadruple precision, ieee binary128 [<a href="https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format">wikipedia</a>]
 120 </ul>
 121 (and other non-standard formats are not supported)
 122 <p>
 123 ld64 is easy to handle: all long double functions
 124 are just wrappers around the corresponding double ones
 125 (aliasing is non-conformant)
 126 <p>
 127 ld80 is the most common (i386, x86_64), it means
 128 64bit significand with explicit msb (inconsistent with other ieee formats),
 129 15bit exp, 1 sign bit.
 130 <p>
 131 ld128 is rare (sparc64 with software emulation), it means
 132 113bit significand with implicit msb, 15bit exp, 1 sign bit.
 133 <p>
 134 Endianness can vary (although on the supported i386 and x86_64 it is the same)
 135 and there is no large enough unsigned int to handle it.
 136 (In case of ld80 internal padding can vary as well, eg
 137 m68k and m88k cpus use different ld80 than the intel ones.)
 138 So there can be further variations in the binary representation
 139 than just ld80 and ld128.
 140
 141
 142 <h3><a name="ugly" href="#ugly">Ugly</a></h3>
 143 <p>The ugly parts of libm hacking.
 144 <p>Some notes are from:
 145 <a href="http://www.vinc17.org/research/extended.en.html">http://www.vinc17.org/research/extended.en.html</a>
 146
 147 <ul>
 148 <li>Double rounding:
 149 <p>
 150 If a value rounded twice the result can be different
 151 than rounding just once.
 152 <p>
 153 The C language allows arithmetic to be evaluated in
 154 higher precision than the operands have. If we
 155 use x87 fpu in extended precision mode it will round
 156 the results twice: round to 80bit when calculating
 157 and then round to 64bit when storing it, this can
 158 give different result than a single 64bit rounding.
 159 (on x86-linux the default fpu setting is to round the
 160 results in extended precision, this only affects x87 instructions, not see2 etc)
 161 (freebsd and openbsd use double precision by default)
 162 <p>
 163 So x = a+b may give different results depending on
 164 the x87 fpu precision setting.
 165 (only happens in round to nearest rounding mode,
 166 but that's the most common one)
 167 <p>
 168 (double rounding can happen with float vs double as well)
 169 <p>
 170 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#F.7.3">C99 annex F</a>
 171 prohibits double rounding, but that's non-normative.
 172
 173 <li>Wider exponent range (x87 issue):
 174 <p>
 175 Even if the fpu is set to double precision
 176 (which is not) the x87 registers use wider exponent
 177 range (mant:exp is 53:15 instead of 53:11 bits)
 178 so underflows (subnormals) may not be treated
 179 as expected. Rounding to double only occurs
 180 when a value is stored into memory.
 181 <p>
 182 Actually this beahviour is allowed by the ieee 754
 183 standard, but it can cause problems (see
 184 <a href="http://www.exploringbinary.com/php-hangs-on-numeric-value-2-2250738585072011e-308/">infamous php bug</a>)
 185
 186 <li>Evaluation precision:
 187 <p>
 188 C does not require consistent evaluation
 189 precision: the compiler may store intermediate
 190 results and round them to double while keep
 191 other parts in higher precision.
 192 (And the precision the compiler choses can
 193 be inconsistent: adding a printf to the code
 194 may change the result of a nearby calculation).
 195 So
 196 <pre>
 197 (a+b)==(a+b)
 198 </pre>
 199 may be false when the two sides
 200 are kept in different precision.
 201 (This is not an x87 specific problem, it matters whenever there
 202 is a higher precision fp type than the currently used one.
 203 It goes away if the highest precision (long double) is used
 204 everywhere, but that can have a huge penalty).
 205 <p>
 206 C99 has a way to control this (see
 207 <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>,
 208 <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
 209 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#6.3.1.8">6.3.1.8</a>):
 210 when a result is stored in a variable or a
 211 type cast is used, then it is guaranteed that
 212 the precision is appropriate to that type.
 213 So in
 214 <pre>
 215 (double)(a+b)==(double)(a+b)
 216 </pre>
 217 both sides are guaranteed to be in double precision
 218 when the comparision is done.
 219 <p>
 220 (This still does not solve the x87 double rounding
 221 issue though: eg if the left side is evaluated with
 222 sse2 and the right side with x87 extended precision setting
 223 and double rounding then the result may still be false)
 224 <p>
 225 Unfortunately gcc does not respect the standard
 226 and even if assingment or cast is used the result
 227 may be kept in higher precision
 228 (infamous <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323">gcc bug323</a>
 229 also see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36578">gcc bug36578</a>).
 230 gcc 4.5 fixed it with '-fexcess-precision=standard'
 231 (it is enabled by '-std=c99', but the default is
 232 '-fexcess-precision=fast')
 233 <p>
 234 The workaround for older gcc is to force the
 235 compiler to store the intermediate results:
 236 by using volatile double temporary variables
 237 or by '-ffloat-store' (which affects all
 238 intermediate results, but is not guaranteed
 239 by the gcc docs to always work).
 240 <p>
 241 (Sometimes the excess precision is good
 242 but it's hard to rely on it as it is optional
 243 <p>
 244 There is a way to check for it though using
 245 FLT_EVAL_METHOD, float_t and double_t.
 246 But it's probably easier to unconditionally
 247 use higher precision variables when that's
 248 what we want and don't depend on the implicit
 249 excess precision).
 250
 251 <li>Compiler optimizations:
 252 <p>
 253 Runtime and compile time semantics may be different
 254 gcc often does unwanted or even invalid compile
 255 time optimizations.
 256 (constant folding where floating point exception
 257 flags should be raised at runtime, or the result
 258 depends on runtime floating point environment,
 259 or constant expressions are just evaluated with
 260 different precision than at runtime).
 261 <p>
 262 C99 actually allows most of these optimizations
 263 but they can be turned off with STDC pragmas (see
 264 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#6.10.6">6.10.6</a>).
 265 Unfortunately <a href="http://gcc.gnu.org/c99status.html">gcc does not support these pragmas</a>.
 266 <p>
 267 FENV_ACCESS ON tells the compiler that the code wants
 268 to access the floating point environment (eg. set different rounding mode)
 269 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34678">gcc bug34678</a>).
 270 <p>
 271 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845">gcc bug37845</a> for FP_CONTRACT pragma
 272 which is relevant when the architecture has fma instruction, x86_64, i386, arm do not have it).
 273 <p>
 274 The workaround is again using named volatile
 275 variables for constants like
 276 <pre>
 277 static const volatile two52 = 0x1p52;
 278 </pre>
 279 and using the '-frounding-math' gcc flag.
 280 <p>
 281 (According the freebsd libm code gcc truncates
 282 long double const literals on i386.
 283 I haven't yet verified if this still the case,
 284 but as a workaround double-double arithmetics is used:
 285 initializing the long double constant from two doubles)
 286
 287 <li>ld80 vs ld128
 288 <p>
 289 The two representations are sufficiently different
 290 that treating them together is awkward.
 291 (Especially the explicit msb bit in ld80 can cause
 292 different behaviour).
 293 <p>
 294 In the freebsd libm code a few architecture specific
 295 macros and a union handle these issues, but the
 296 result is often less clear than treating ld80
 297 and ld128 separately.
 298
 299 <li>Signed zeros
 300 <p>Signed zeros can be tricky. They cannot be checked
 301 using the usual comparision operators (+0.0 == -0.0 is true),
 302 but they give different results in some cases
 303 (1/-0.0 == -inf) and can make similarly looking
 304 expressions different eg.
 305 (x + 0) is not the same as x, the former is +0.0 when x is -0.0
 306 <p>
 307 (To check for -0, the signbit macro can be used
 308 or the copysign function or bit manipulation.)
 309
 310 <li>Error handling
 311 <p>
 312 Arithmetics may set various floating point exception flags as a side effect.
 313 These can be queried and manipulated (fetestexcept, feraiseexcept,..).
 314 <p>
 315 So special care is needed
 316 when a library function wants to avoid changing
 317 the floating point status flags.
 318 eg. if one wants to check for -0 silently then
 319 <pre>
 320 if (x == 0.0 &amp;&amp; 1/x &lt; 0) { /* x is a -0 */ }
 321 </pre>
 322 is not ok: == raises invalid exception when x is nan,
 323 and even if we filter nans out 1/x will raise the
 324 divbyzero flag.
 325 <p>
 326 When a library wants to raise a flag deliberately
 327 but feraiseexcept is not available for some reason,
 328 then simple arithmetics can be be used just for their
 329 exception raising side effect
 330 (eg. 1/0.0 to raise divbyzero), however beaware
 331 of compiler optimizations (dead code elimination,..).
 332 <p>
 333 Unfortunately gcc does not always take fp exceptions into
 334 account: a simple x = 1e300*1e300; may not raise overflow
 335 exception at runtime, but get optimized into x = +inf.
 336 see compiler optimizations above.
 337 <p>
 338 Another x87 gcc bug related to fp exceptions is that
 339 comparision operators (==, &lt;, etc) don't raise invalid
 340 when an operand is nan
 341 (eventhough this is required by ieee + c99 annex F).
 342 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52451">gcc bug52451</a>).
 343 <p>
 344 The ieee standard defines signaling and quite nan
 345 floating-point numbers as well.
 346 The c99 standard only considers quiet nan, but it allows
 347 signaling nans to be supported as well.
 348 Without signaling nans x * 1 is equivalent to x,
 349 but if signaling nan is supported then the former
 350 raises an invalid exception.
 351 This may complicates things further if one wants to write
 352 portable fp math code.
 353 <p>
 354 A further libm design issue is the math_errhandling macro:
 355 it specifies the way math function errors can be checked
 356 which is either fp exceptions or the errno variable or both.
 357 The fp exception approach can be supported with enough care
 358 but errno is hard to support: certain library functions
 359 are implemented as a single asm instruction (eg sqrt),
 360 the only way to set errno is to query the fp exception flags
 361 and then set the errno variable based on that.
 362 So eventhough errno may be convenient in libm it is
 363 not the right thing to do.
 364
 365 <li>Complex arithmetics
 366 <p>
 367 gcc turns
 368 <pre>
 369 x*I
 370 </pre>
 371 into
 372 <pre>
 373 (x+0*I)*(0+I) = 0*x + x*I
 374 </pre>
 375 So if x=inf then the result is nan+inf*I instead of inf*I
 376 (an fmul instruction is generated for 0*x)
 377 <p>
 378 (There were various complex constant folding issues as well in gcc).
 379 <p>
 380 So a+b*I cannot be used to create complex numbers,
 381 instead some double[2] manipulation should be used to
 382 make sure there is no fmul in the generated code
 383 (freebsd uses a static inline cpack function)
 384 (It's not clear which is better: using union hacks
 385 everywhere in the complex code or creal, cimag, cpack)
 386 <p>
 387 Complex code is hard to make portable across compilers
 388 as there are no good compiler support, it is rarely used
 389 and there are various options available in the standard:
 390 complex support is optional, imaginary type
 391 is optional.
 392 <p>
 393 It's not clear how I (or _Complex_I) should be defined
 394 in complex.h
 395 (literal of the float imaginary unit is compiler specific,
 396 in gcc it can be 1.0fi).
 397
 398 <li>Signed int
 399 <p>
 400 The freebsd libm code has many inconsistencies
 401 (naming conventions, 0x1p0 notation vs decimal notation,..),
 402 one of them is the integer type used for bitmanipulations:
 403 The bits of a double are unpacked into one of
 404 int32_t, uint32_t and u_int32_t
 405 integer types.
 406 <p>
 407 int32_t is used most often which is wrong because of
 408 implementation defined signed int representation.
 409 <p>
 410 In general signed int is not handled carefully
 411 in the libm code: scalbn even depends on signed int overflow.
 412 </ul>
 413
 414 <h3><a name="implementations" href="#implementations">libm implementations</a></h3>
 415 <ul>
 416 <li>
 417 unix v7 libm by Robert Morris (rhm) around 1979<br>
 418 see: "Computer Approximations" by Hart & Cheney, 1968<br>
 419 used in: plan9, go
 420 <li>
 421 cephes by Stephen L. Moshier, 1984-1992<br>
 422 <a href="http://www.netlib.org/cephes/">http://www.netlib.org/cephes/</a><br>
 423 see: "Methods and Programs for Mathematical Functions" by Moshier 1989
 424 <li>
 425 fdlibm by K. C. Ng around 1993, updated in 1995?<br>
 426 (1993: copyright notice: "Developed at SunPro ..")<br>
 427 (1995: copyright notice: "Developed at SunSoft ..")<br>
 428 eg see: "Argument Reduction for Huge Arguments: Good to the Last Bit" by Ng 1992<br>
 429 <a href="http://www.netlib.org/fdlibm">http://www.netlib.org/fdlibm</a><br>
 430 used in: netbsd, openbsd, freebsd, bionic, musl
 431 <li>
 432 libultim by Abraham Ziv around 2001 (lgpl)<br>
 433 (also known as the ibm accurate portable math lib)<br>
 434 see: "Fast evaluation of elementary mathematical functions with correctly rounded last bit" by Ziv 1991<br>
 435 <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>
 436 used in: glibc
 437 <li>
 438 libmcr by K. C. Ng, 2004<br>
 439 (sun's correctly rounded libm)
 440 <li>
 441 mathcw by Nelson H. F. Beebe, 2010?<br>
 442 <a href="http://ftp.math.utah.edu/pub/mathcw/">http://ftp.math.utah.edu/pub/mathcw/</a><br>
 443 (no sources available?)
 444 <li>
 445 sleef by Naoki Shibata, 2010<br>
 446 (simd lib for evaluating elementary functions)<br>
 447 <a href="http://shibatch.sourceforge.net/">http://shibatch.sourceforge.net/</a>
 448 <li>
 449 crlibm by ens-lyon, 2004-2010<br>
 450 <a href="http://lipforge.ens-lyon.fr/www/crlibm/">http://lipforge.ens-lyon.fr/www/crlibm/</a><br>
 451 see: <a href="http://lipforge.ens-lyon.fr/frs/?group_id=8&amp;release_id=123">crlibm.pdf</a><br>
 452 see: "Elementary Functions" by Jean-Michel Muller 2005<br>
 453 see: "Handbook of Floating-Point Arithmetic" by (many authors from Lyon) 2009<br>
 454 <li>
 455 various other closed ones: intel libm, hp libm for itanium,..
 456 </ul>
 457
 458 <h3><a name="tests" href="#tests">libm tests</a></h3>
 459 <ul>
 460 <li><a href="http://www.netlib.org/fp/ucbtest.tgz">ucbtest.tgz</a>
 461 <li><a href="http://www.jhauser.us/arithmetic/TestFloat.html">TestFloat</a>
 462 <li><a href="http://cant.ua.ac.be/old/ieeecc754.html">ieeecc754</a>
 463 <li><a href="http://www.loria.fr/~zimmerma/mpcheck/">mpcheck</a>,
 464 <a href="http://gforge.inria.fr/projects/mpcheck/">mpcheck</a>
 465 <li><a href="http://people.inf.ethz.ch/gonnet/FPAccuracy/Analysis.html">FPAccuracy</a>
 466 <li><a href="http://www.math.utah.edu/~beebe/software/ieee/">beebe's ieee tests</a>
 467 <li><a href="http://www.math.utah.edu/pub/elefunt/">elefunt</a>
 468 <li><a href="http://www.vinc17.org/research/testlibm/index.en.html">testlibm</a>
 469 <li><a href="http://www.vinc17.org/research/fptest.en.html">fptest</a>
 470 <li>tests in crlibm
 471 </ul>
 472 <p>
 473 multiprecision libs (useful for tests)
 474 <ul>
 475 <li><a href="http://mpfr.org/">mpfr</a>
 476 <li><a href="http://www.multiprecision.org/index.php?prog=mpc">mpc</a>
 477 <li><a href="http://pari.math.u-bordeaux.fr/">pari</a>
 478 <li><a href="http://www.apfloat.org/apfloat/">apfloat</a>
 479 <li><a href="http://code.google.com/p/fastfunlib/">fastfunlib</a>
 480 <li>scs_lib in crlibm
 481 </ul>
 482
 483 <p>other links
 484 <ul>
 485 <li>ieee standard: http://754r.ucbtest.org/
 486 <li>extended precision issues: http://www.vinc17.org/research/extended.en.html
 487 <li>correctly rounded mult: http://perso.ens-lyon.fr/jean-michel.muller/MultConstant.html
 488 <li>table based sincos: http://perso.ens-lyon.fr/damien.stehle/IMPROVEDGAL.html
 489 <li>finding worst cases: http://perso.ens-lyon.fr/damien.stehle/WCLR.html
 490 <li>finding worst cases: http://perso.ens-lyon.fr/damien.stehle/DECIMALEXP.html
 491 <li>finding worst cases: http://www.loria.fr/equipes/spaces/slz.en.html
 492 <li>finding worst cases: http://perso.ens-lyon.fr/jean-michel.muller/Intro-to-TMD.htm
 493 <li>generating libm functions: http://lipforge.ens-lyon.fr/www/metalibm/
 494 <li>fast conversion to fixedpoint: http://stereopsis.com/sree/fpu2006.html
 495 <li>double-double, quad-double arithmetics: http://crd-legacy.lbl.gov/~dhbailey/mpdist/
 496 <li>papers by kahan: http://www.cs.berkeley.edu/~wkahan/
 497 <li>fp paper by goldberg: http://download.oracle.com/docs/cd/E19957-01/806-3568/ncg_goldberg.html
 498 <li>math functions in general: http://dlmf.nist.gov/
 499 </ul>
 500
 501 <p><small>(<a href="/git/?p=www">page history</a>, <a href="/">home</a>)</small>
 502 </body></html>