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