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="#sources">Sources</a>
  11 <li><a href="#issues">Design issues</a>
  12 <li><a href="#representation">Representation</a>
  13 <li><a href="#ugly">Ugly</a>
  14 <li><a href="#implementations">libm implementations</a>
  15 <li><a href="#tests">libm tests</a>
  16 </ul>
  17
  18 <h3><a name="sources" href="#sources">Sources</a></h3>
  19 <p>The math code is mostly from freebsd which in turn
  20 is based on <a href="http://www.netlib.org/fdlibm/">fdlibm</a>.
  21 <p>sources:
  22 <ul>
  23 <li>freebsd /lib/msun (<a href="http://svnweb.FreeBSD.org/base/head/lib/msun/">browse</a>)
  24 <pre>
  25 svn checkout svn://svn.freebsd.org/base/head/lib/msun
  26 </pre>
  27 <li>openbsd /src/lib/libm (<a href="http://www.openbsd.org/cgi-bin/cvsweb/src/lib/libm/">browse</a>)
  28 <pre>
  29 cvs -d <a href="http://openbsd.org/anoncvs.html#CVSROOT">$CVSROOT</a> get src/lib/libm
  30 </pre>
  31 </ul>
  32
  33 <h3><a name="issues" href="#issues">Design issues</a></h3>
  34 <p>The math code is available but there are still many design questions.
  35 <ul>
  36 <li>Code organization, build system tweaks:
  37 <p>There are architecture specific code,
  38 long double representation specific code
  39 and complex code.
  40 <p>With respect to long double, bsd libm code can handle
  41 ld64, ld80 and ld128 (where ld64 means long double == double),
  42 for musl it's enough to support ld64 and ld80, but
  43 keeping ld128 shouldn't hurt much either.
  44 <p>So far musl had arch specific code and common code,
  45 it seems a new category is needed:
  46 long double specific code, which can be ld64, ld80 or ld128.
  47 These should be tied to the architecture in the build system
  48 somehow. In freebsd most of the long double code is
  49 common (using arch specific macros, ifdefs) and only a
  50 small part of the ld code is maintained separately for ld80 and ld128,
  51 but it's also possible to keep all ld code separately (with duplications)
  52 in ld80 and ld128.
  53 <p>Complex support is optional in c99 so maybe it should
  54 be optional in the musl build as well. It also can be
  55 kept separately in cmath/ or with the rest of the code in math/.
  56 <p>Pending questions:
  57 <ul>
  58 <li>Do we want ld128?
  59 <li>Should we try to use common code for ld80 and ld128?
  60 <li>How to do ld64: wrap double functions or alias them?
  61 <li>How to tie the ld* code to the arch in the build system?
  62 <li>Make complex optional?
  63 <li>Keep complex in math/ or cmath/?
  64 </ul>
  65
  66 <li>Workarounds
  67 <p>The bsd libm code has many workarounds for various
  68 compiler issues. It's not clear what's the best way
  69 to handle the uglyness.
  70 <p>Pending questions:
  71 <ul>
  72 <li>Use STDC pragmas (eventhough gcc does not support them)?
  73 <li>Use volatile consistently to avoid evaluation precision and const folding issues?
  74 <li>Use special compiler flags against unwanted optimization (-frounding-math, -ffloat-store)?
  75 <li>Do inline/macro optimization for small functions? (isnan, isinf, signbit, creal, cimag,..)
  76 <li>In complex code prefer creal(), cimag() or a union to (un)pack re,im?
  77 </ul>
  78
  79 <li>Code cleanups
  80 <p>The fdlibm code style is inconsistent with musl (and ugly)
  81 so it makes sense to clean it up, but keeping easy diffability
  82 can be useful. There are various inconsistencies and
  83 correctness issues in the code which might worth addressing
  84 (eg. reliance on signed overflow should be eliminated).
  85 <p>Pending questions:
  86 <ul>
  87 <li>Keep diffability with freebsd/openbsd code or reformat the code for clarity?
  88 <li>Keep e_, s_, k_  fdlibm source file name prefixes?
  89 <li>Should 0x1p0 float format be preferred over decimal format?
  90 <li>Should isnan, signbit,.. be preferred over inplace bithacks?
  91 <li>Is unpacking a double into 2 int32_t ok (signed int representation)?
  92 <li>Is unpacking the mantissa of ld80 into an int64_t ok?
  93 </ul>
  94 </ul>
  95
  96 <h3><a name="representation" href="#representation">Representation</a></h3>
  97 <p>
  98 Binary representation of floating point numbers matter
  99 because lots of bithacks are needed in the math code.
 100 <p>
 101 float and double bit manipulation can be handled
 102 in a portable way in c:
 103 <ul>
 104 <li>float is ieee binary32
 105 <li>double is ieee binary64
 106 </ul>
 107 (and old non-standard representations are not supported)
 108 <p>
 109 The endianness may still vary, but that can be worked
 110 around by using a union with a single large enough
 111 unsigned int. (The only exception is probably arm/oabi fpa
 112 where the word order of a double is not consistent with the
 113 endianness [<a href="http://wiki.debian.org/ArmEabiPort#ARM_floating_points">debian wiki on arm</a>],
 114 but we probably won't support that)
 115 <p>
 116 long double bit manipulation is harder as there are
 117 various representations:
 118 <ul>
 119 <li>ld64: long double is the same as double (ieee binary64)
 120 <li>ld80: 80bit extended precision format [<a href="https://en.wikipedia.org/wiki/Extended_precision">wikipedia</a>]
 121 <li>ld128: quadruple precision, ieee binary128 [<a href="https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format">wikipedia</a>]
 122 </ul>
 123 (and other non-standard formats are not supported)
 124 <p>
 125 ld64 is easy to handle: all long double functions
 126 are aliased to the corresponding double one
 127 (long double is treated equivalently to double)
 128 <p>
 129 ld80 is the most common (i386, x86_64), it means
 130 64bit significand with explicit msb (inconsistent with other ieee formats),
 131 15bit exp, 1 sign bit.
 132 <p>
 133 ld128 is rare (sparc64 with software emulation), it means
 134 113bit significand with implicit msb, 15bit exp, 1 sign bit.
 135 <p>
 136 Endianness can vary (although on the supported i386 and x86_64 it is the same)
 137 and there is no large enough unsigned int to handle it.
 138 (In case of ld80 internal padding can vary as well, eg
 139 m68k and m88k cpus use different ld80 than the intel ones)
 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
 148 <ul>
 149 <li>Double rounding (x87 issue):
 150 <p>
 151 If a value rounded twice the result can be different
 152 than rounding just once.
 153 <p>
 154 The C language allows arithmetic to be evaluated in
 155 higher precision than the operands have. If we
 156 use x87 fpu in extended precision mode it will round
 157 the results twice: round to 80bit when calculating
 158 and then round to 64bit when storing it, this can
 159 give different result than a single 64bit rounding.
 160 (on x86-linux the default fpu setting is to round the
 161 results in extended precision, this only affects x87 instructions, not see2 etc)
 162 (afaik freebsd and openbsd use double precision by default)
 163 <p>
 164 So x = a+b may give different results depending on
 165 the fpu setting.
 166 (only happens in round to nearest rounding mode,
 167 but that's the most common one)
 168 <p>
 169 C99 annex F prohibits double rounding,
 170 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 So (a+b)==(a+b) may be false when the two sides
 192 are kept in different precision.
 193 (This is not an x87 specific problem, it matters whenever there
 194 is a higher precision fp type than the currently used one.
 195 It goes away if the highest precision (long double) is used
 196 everywhere, but that can have a huge penalty).
 197 <p>
 198 C99 has a way to control this (see
 199 <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>,
 200 <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
 201 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#6.3.1.8">6.3.1.8</a>):
 202 when a result is stored in a variable or a
 203 type cast is used, then it is guaranteed that
 204 the precision is appropriate to that type.
 205 So in (double)(a+b)==(double)(a+b) both
 206 sides are guaranteed to be in double precision
 207 when the comparision is done.
 208 <p>
 209 (This still does not solve the x87 double rounding
 210 issue though: eg if the left side is evaluated with
 211 sse2 and the right side with x87 extended precision setting
 212 and double rounding then the result may be false
 213 and still conformant to C99)
 214 <p>
 215 Unfortunately gcc does not respect the standard
 216 and even if assingment or cast is used the result
 217 may be kept in higher precision
 218 (infamous <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323">gcc bug323</a>
 219 also see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36578">gcc bug36578</a>).
 220 gcc 4.5 fixed it with '-fexcess-precision=standard'
 221 (it is enabled by '-std=c99', but the default is
 222 '-fexcess-precision=fast')
 223 <p>
 224 The workaround for older gcc is to force the
 225 compiler to store the intermediate results:
 226 by using volatile double temporary variables
 227 or by '-ffloat-store' (slow, and of course
 228 all intermediate results should be assigned
 229 to some variable, casting is not enough).
 230 <p>
 231 (Sometimes the excess precision is good
 232 but it's hard to rely on it as it is optional
 233 <p>
 234 There is a way to check for it though using
 235 FLT_EVAL_METHOD, float_t and double_t.
 236 But it's probably easier to unconditionally
 237 use higher precision variables when that's
 238 what we want and don't depend on the implicit
 239 excess precision).
 240
 241 <li>Compiler optimizations:
 242 <p>
 243 Runtime and compile time semantics may be different
 244 gcc often does unwanted or even invalid compile
 245 time optimizations.
 246 (constant folding where floating point exception
 247 flags should be raised at runtime, or the result
 248 depends on runtime floating point environment,
 249 or constant expressions are just evaluated with
 250 different precision than at runtime).
 251 <p>
 252 C99 actually allows most of these optimizations
 253 but they can be turned off with STDC pragmas (see
 254 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#6.10.6">6.10.6</a>).
 255 Unfortunately <a href="http://gcc.gnu.org/c99status.html">gcc does not support these pragmas</s>.
 256 <p>
 257 FENV_ACCESS ON tells the compiler that the code wants
 258 to access the floating point environment (eg. set different rounding mode)
 259 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34678">gcc bug34678</a>).
 260 <p>
 261 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845">gcc bug37845</a> for FP_CONTRACT pragma).
 262 <p>
 263 The workaround is again using named volatile
 264 variables for constants like
 265 <pre>
 266 static const volatile two52 = 0x1p52;
 267 </pre>
 268 and using the '-frounding-math' gcc flag.
 269 <p>
 270 (According the freebsd libm code gcc truncates
 271 long double const literals on i386.
 272 I haven't yet verified if this still the case,
 273 but as a workaround sometimes double-double arithmetics is used:
 274 initializing the long double constant from two doubles)
 275
 276 <li>ld80 vs ld128
 277 <p>
 278 The two representations are sufficiently different
 279 that treating them together is awkward.
 280 <p>
 281 In the freebsd libm code a few architecture specific
 282 macros and a union handle these issues, but the
 283 result is often less clear than treating ld80
 284 and ld128 separately.
 285
 286 <li>Signed int
 287 <p>
 288 The freebsd libm code has many inconsistencies
 289 (naming conventions, 0x1p0 notation vs decimal notation,..),
 290 one of them is the integer type used for bitmanipulations:
 291 The bits of a double are unpacked into one of
 292 int32_t, uint32_t and u_int32_t
 293 integer types.
 294 <p>
 295 int32_t is used most often which is wrong because of
 296 implementation defined signed int representation.
 297 <p>
 298 In general signed int is not handled carefully
 299 in the libm code: scalbn even depends on signed int overflow.
 300
 301 <li>Complex arithmetics
 302 <p>
 303 With gcc x*I turns into (x+0*I)*(0+I) = 0*x + x*I,
 304 so if x=inf then the result is nan+inf*I instead of inf*I
 305 (an fmul instruction is generated for 0*x)
 306 <p>
 307 (There were various complex constant folding issues as well).
 308 <p>
 309 So a+b*I cannot be used to create complex numbers,
 310 instead some double[2] manipulation should be used to
 311 make sure there is no fmul in the generated code
 312 (freebsd uses a static inline cpack function)
 313 (It's not clear which is better: using union hacks
 314 everywhere in the complex code or creal, cimag, cpack)
 315 <p>
 316 Complex code is hard to make portable across compilers
 317 as there are no good compiler support, it is rarely used
 318 and there are various options available in the standard:
 319 complex support is optional, imaginary type
 320 is optional.
 321 <p>
 322 It's not clear how I (or _Complex_I) should be defined
 323 in complex.h
 324 (literal of the float imaginary unit is compiler specific,
 325 in gcc it can be 1.0fi).
 326 </ul>
 327
 328 <h3><a name="implementations" href="#implementations">libm implementations</a></h3>
 329 <ul>
 330 <li>
 331 unix v7 libm by Robert Morris (rhm) around 1979<br>
 332 see: "Computer Approximations" by Hart & Cheney, 1968<br>
 333 used in: plan9, go
 334 <li>
 335 cephes by Stephen L. Moshier, 1984-1992<br>
 336 <a href="http://www.netlib.org/cephes/">http://www.netlib.org/cephes/</a><br>
 337 see: "Methods and Programs for Mathematical Functions" by Moshier 1989
 338 <li>
 339 fdlibm by K. C. Ng around 1993, updated in 1995?<br>
 340 (1993: copyright notice: "Developed at SunPro ..")<br>
 341 (1995: copyright notice: "Developed at SunSoft ..")<br>
 342 eg see: "Argument Reduction for Huge Arguments: Good to the Last Bit" by Ng 1992<br>
 343 <a href="http://www.netlib.org/fdlibm">http://www.netlib.org/fdlibm</a><br>
 344 used in: netbsd, openbsd, freebsd, bionic, musl
 345 <li>
 346 libultim by Abraham Ziv around 2001 (lgpl)<br>
 347 (also known as the ibm accurate portable math lib)<br>
 348 see: "Fast evaluation of elementary mathematical functions with correctly rounded last bit" by Ziv 1991<br>
 349 <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>
 350 used in: glibc
 351 <li>
 352 libmcr by K. C. Ng, 2004<br>
 353 (sun's correctly rounded libm)
 354 <li>
 355 mathcw by Nelson H. F. Beebe, 2010?<br>
 356 <a href="http://ftp.math.utah.edu/pub/mathcw/">http://ftp.math.utah.edu/pub/mathcw/</a><br>
 357 (no sources available?)
 358 <li>
 359 sleef by Naoki Shibata, 2010<br>
 360 (simd lib for evaluating elementary functions)<br>
 361 <a href="http://shibatch.sourceforge.net/">http://shibatch.sourceforge.net/</a>
 362 <li>
 363 crlibm by ens-lyon, 2004-2010<br>
 364 <a href="http://lipforge.ens-lyon.fr/www/crlibm/">http://lipforge.ens-lyon.fr/www/crlibm/</a><br>
 365 see: <a href="http://lipforge.ens-lyon.fr/frs/?group_id=8&amp;release_id=123">crlibm.pdf</a><br>
 366 see: "Elementary Functions" by Jean-Michel Muller 2005<br>
 367 see: "Handbook of Floating-Point Arithmetic" by (many authors from Lyon) 2009<br>
 368 <li>
 369 various other closed ones: intel libm, hp libm for itanium,..
 370 </ul>
 371
 372 <h3><a name="tests" href="#tests">libm tests</a></h3>
 373 <ul>
 374 <li><a href="http://www.netlib.org/fp/ucbtest.tgz">ucbtest.tgz</a>
 375 <li><a href="http://www.jhauser.us/arithmetic/TestFloat.html">TestFloat</a>
 376 <li><a href="http://cant.ua.ac.be/old/ieeecc754.html">ieeecc754</a>
 377 <li><a href="http://www.loria.fr/~zimmerma/mpcheck/">mpcheck</a>,
 378 <a href="http://gforge.inria.fr/projects/mpcheck/">mpcheck</a>
 379 <li><a href="http://people.inf.ethz.ch/gonnet/FPAccuracy/Analysis.html">FPAccuracy</a>
 380 <li><a href="http://www.math.utah.edu/~beebe/software/ieee/">beebe's ieee tests</a>
 381 <li><a href="http://www.vinc17.org/research/testlibm/index.en.html">testlibm</a>
 382 <li><a href="http://www.vinc17.org/research/fptest.en.html">fptest</a>
 383 <li>tests in crlibm
 384 </ul>
 385 <p>
 386 multiprecision libs (useful for tests)
 387 <ul>
 388 <li><a href="http://mpfr.org/">mpfr</a>
 389 <li><a href="http://www.multiprecision.org/index.php?prog=mpc">mpc</a>
 390 <li><a href="http://pari.math.u-bordeaux.fr/">pari</a>
 391 <li><a href="http://www.apfloat.org/apfloat/">apfloat</a>
 392 <li><a href="http://code.google.com/p/fastfunlib/">fastfunlib</a>
 393 <li>scs_lib in crlibm
 394 </ul>
 395
 396 <p>other links
 397 <ul>
 398 <li>ieee standard: http://754r.ucbtest.org/
 399 <li>extended precision issues: http://www.vinc17.org/research/extended.en.html
 400 <li>correctly rounded mult: http://perso.ens-lyon.fr/jean-michel.muller/MultConstant.html
 401 <li>table based sincos: http://perso.ens-lyon.fr/damien.stehle/IMPROVEDGAL.html
 402 <li>finding worst cases: http://perso.ens-lyon.fr/damien.stehle/WCLR.html
 403 <li>finding worst cases: http://perso.ens-lyon.fr/damien.stehle/DECIMALEXP.html
 404 <li>finding worst cases: http://www.loria.fr/equipes/spaces/slz.en.html
 405 <li>finding worst cases: http://perso.ens-lyon.fr/jean-michel.muller/Intro-to-TMD.htm
 406 <li>generating libm functions: http://lipforge.ens-lyon.fr/www/metalibm/
 407 <li>fast conversion to fixedpoint: http://stereopsis.com/sree/fpu2006.html
 408 <li>double-double, quad-double arithmetics: http://crd-legacy.lbl.gov/~dhbailey/mpdist/
 409 <li>papers by kahan: http://www.cs.berkeley.edu/~wkahan/
 410 <li>fp paper by goldberg: http://download.oracle.com/docs/cd/E19957-01/806-3568/ncg_goldberg.html
 411 <li>math functions in general: http://dlmf.nist.gov/
 412 </ul>
 413
 414 <p><small>(<a href="/git/?p=www">page history</a>, <a href="/">home</a>)</small>
 415 </body></html>