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