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>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>):
  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 for most functions.
  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, logb,
  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
  73 may break, this also means that compiler must spill fpu registers at 80bits
  74 precision)
  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
 106 implementations.
 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.
 110 </ul>
 111
 112 <h3><a name="representation" href="#representation">Representation</a></h3>
 113 <p>
 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)
 119 <p>
 120 float and double bit manipulation can be handled in a portable way in c using
 121 union types:
 122 <ul>
 123 <li>union {float f; uint32_t i;};
 124 <li>union {double f; uint64_t i;};
 125 </ul>
 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
 131 arm ABI)
 132 <p>
 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:
 135 <ul>
 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>]
 139 </ul>
 140 (and other non-standard formats are not supported)
 141 <p>
 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)
 147 <p>
 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:
 152 <ul>
 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
 155 </ul>
 156 where m is the significand and se is the sign and exponent.
 157 <p>
 158 ld128 is rare (eg. sparc64 with software emulation), it means
 159 113bit significand with implicit msb, 15bit exp, 1 sign bit:
 160 <ul>
 161 <li>union {long double f; struct{uint16_t se; uint16_t hi; uint32_t mid; uint64_t lo;} i;};
 162 </ul>
 163 <p>
 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
 166 implemented using
 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
 171 representation.
 172 <p>
 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.
 176
 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>
 181 <p>Useful info about floating-point in gcc:
 182 <a href="http://gcc.gnu.org/wiki/FloatingPointMath">http://gcc.gnu.org/wiki/FloatingPointMath</a>
 183
 184 <ul>
 185 <li>Double rounding:
 186 <p>
 187 If a value rounded twice the result can be different
 188 than rounding just once.
 189 <p>
 190 The C language allows arithmetic to be evaluated in
 191 higher precision than the operands have. If we
 192 use x87 fpu in extended precision mode it will round
 193 the results twice: round to 80bit when calculating
 194 and then round to 64bit when storing it, this can
 195 give different result than a single 64bit rounding.
 196 (on x86-linux the default fpu setting is to round the
 197 results in extended precision, this only affects x87 instructions, not see2 etc)
 198 (freebsd and openbsd use double precision by default)
 199 <p>
 200 So x = a+b may give different results depending on
 201 the x87 fpu precision setting.
 202 (only happens in round to nearest rounding mode,
 203 but that's the most common one)
 204 <p>
 205 (double rounding can happen with float vs double as well)
 206 <p>
 207 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#F.7.3">C99 annex F</a>
 208 prohibits double rounding, but that's non-normative.
 209
 210 <li>Wider exponent range (x87 issue):
 211 <p>
 212 Even if the fpu is set to double precision
 213 (which is not) the x87 registers use wider exponent
 214 range (mant:exp is 53:15 instead of 53:11 bits)
 215 so underflows (subnormals) may not be treated
 216 as expected. Rounding to double only occurs
 217 when a value is stored into memory.
 218 <p>
 219 Actually this beahviour is allowed by the ieee 754
 220 standard, but it can cause problems (see
 221 <a href="http://www.exploringbinary.com/php-hangs-on-numeric-value-2-2250738585072011e-308/">infamous php bug</a>)
 222
 223 <li>Evaluation precision:
 224 <p>
 225 C does not require consistent evaluation
 226 precision: the compiler may store intermediate
 227 results and round them to double while keep
 228 other parts in higher precision.
 229 (And the precision the compiler choses can
 230 be inconsistent: adding a printf to the code
 231 may change the result of a nearby calculation).
 232 So
 233 <pre>
 234 (a+b)==(a+b)
 235 </pre>
 236 may be false when the two sides
 237 are kept in different precision.
 238 (This is not an x87 specific problem, it matters whenever there
 239 is a higher precision fp type than the currently used one.
 240 It goes away if the highest precision (long double) is used
 241 everywhere, but that can have a huge penalty).
 242 <p>
 243 C99 has a way to control this (see
 244 <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>,
 245 <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
 246 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#6.3.1.8">6.3.1.8</a>):
 247 when a result is stored in a variable or a
 248 type cast is used, then it is guaranteed that
 249 the precision is appropriate to that type.
 250 So in
 251 <pre>
 252 (double)(a+b)==(double)(a+b)
 253 </pre>
 254 both sides are guaranteed to be in double precision
 255 when the comparision is done.
 256 <p>
 257 (This still does not solve the x87 double rounding
 258 issue though: eg if the left side is evaluated with
 259 sse2 and the right side with x87 extended precision setting
 260 and double rounding then the result may still be false)
 261 <p>
 262 Unfortunately gcc does not respect the standard
 263 and even if assingment or cast is used the result
 264 may be kept in higher precision
 265 (infamous <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323">gcc bug323</a>
 266 also see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=36578">gcc bug36578</a>).
 267 gcc 4.5 fixed it with '-fexcess-precision=standard'
 268 (it is enabled by '-std=c99', but the default is
 269 '-fexcess-precision=fast')
 270 <p>
 271 The workaround for older gcc is to force the
 272 compiler to store the intermediate results:
 273 by using volatile double temporary variables
 274 or by '-ffloat-store' (which affects all
 275 intermediate results, but is not guaranteed
 276 by the gcc docs to always work).
 277 <p>
 278 (Sometimes the excess precision is good
 279 but it's hard to rely on it as it is optional
 280 <p>
 281 There is a way to check for it though using
 282 FLT_EVAL_METHOD, float_t and double_t.
 283 But it's probably easier to unconditionally
 284 use higher precision variables when that's
 285 what we want and don't depend on the implicit
 286 excess precision).
 287
 288 <li>Compiler optimizations:
 289 <p>
 290 Runtime and compile time semantics may be different
 291 gcc often does unwanted or even invalid compile
 292 time optimizations.
 293 (constant folding where floating point exception
 294 flags should be raised at runtime, or the result
 295 depends on runtime floating point environment,
 296 or constant expressions are just evaluated with
 297 different precision than at runtime).
 298 <p>
 299 C99 actually allows most of these optimizations
 300 but they can be turned off with STDC pragmas (see
 301 <a href="http://repo.or.cz/w/c-standard.git/blob_plain/HEAD:/n1256.html#6.10.6">6.10.6</a>).
 302 Unfortunately <a href="http://gcc.gnu.org/c99status.html">gcc does not support these pragmas</a>.
 303 <p>
 304 FENV_ACCESS ON tells the compiler that the code wants
 305 to access the floating point environment (eg. set different rounding mode)
 306 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=34678">gcc bug34678</a>).
 307 <p>
 308 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845">gcc bug37845</a> for FP_CONTRACT pragma
 309 which is relevant when the architecture has fma instruction, x86_64, i386, arm do not have it).
 310 <p>
 311 The workaround is again using named volatile
 312 variables for constants like
 313 <pre>
 314 static const volatile two52 = 0x1p52;
 315 </pre>
 316 and using the '-frounding-math' gcc flag.
 317 <p>
 318 According to the freebsd libm code gcc truncates long double
 319 const literals on i386.
 320 I assume this happens because freebsd uses 64bit long doubles by default
 321 (double precision) and gcc incorrectly uses the precision setting of the
 322 host platform instead of the target one, but i did not observe this on linux.
 323 (as a workaround sometimes double-double arithmetics was used
 324 to initialize long doubles on i386, but most of these should be
 325 fixed in musl's math code now)
 326
 327 <li>ld80 vs ld128
 328 <p>
 329 The two representations are sufficiently different
 330 that treating them together is awkward.
 331 (Especially the explicit msb bit in ld80 can cause
 332 different behaviour).
 333 <p>
 334 In the freebsd libm code a few architecture specific
 335 macros and a union handle these issues, but the
 336 result is often less clear than treating ld80
 337 and ld128 separately.
 338
 339 <li>Signed zeros
 340 <p>Signed zeros can be tricky. They cannot be checked
 341 using the usual comparision operators (+0.0 == -0.0 is true),
 342 but they give different results in some cases
 343 (1/-0.0 == -inf) and can make similarly looking
 344 expressions different eg.
 345 (x + 0) is not the same as x, the former is +0.0 when x is -0.0
 346 <p>
 347 (To check for -0, the signbit macro can be used
 348 or the copysign function or bit manipulation.)
 349
 350 <li>Error handling
 351 <p>
 352 Arithmetics may set various floating point exception flags as a side effect.
 353 These can be queried and manipulated (fetestexcept, feraiseexcept,..).
 354 <p>
 355 So special care is needed
 356 when a library function wants to avoid changing
 357 the floating point status flags.
 358 eg. if one wants to check for -0 silently then
 359 <pre>
 360 if (x == 0.0 &amp;&amp; 1/x &lt; 0) { /* x is a -0 */ }
 361 </pre>
 362 is not ok: == raises invalid exception when x is nan,
 363 and even if we filter nans out 1/x will raise the
 364 divbyzero flag.
 365 <p>
 366 When a library wants to raise a flag deliberately
 367 but feraiseexcept is not available for some reason,
 368 then simple arithmetics can be be used just for their
 369 exception raising side effect
 370 (eg. 1/0.0 to raise divbyzero), however beaware
 371 of compiler optimizations (constant folding and dead code elimination,..).
 372 <p>
 373 Unfortunately gcc does not always take fp exceptions into
 374 account: a simple x = 1e300*1e300; may not raise overflow
 375 exception at runtime, but get optimized into x = +inf.
 376 see compiler optimizations above.
 377 <p>
 378 Another x87 gcc bug related to fp exceptions is that in some cases
 379 comparision operators (==, &lt;, etc) don't raise invalid
 380 when an operand is nan
 381 (eventhough this is required by ieee + c99 annex F).
 382 (see <a href="http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52451">gcc bug52451</a>).
 383 <p>
 384 The ieee standard defines signaling and quiet nan
 385 floating-point numbers as well.
 386 The c99 standard only considers quiet nan, but it allows
 387 signaling nans to be supported as well.
 388 Without signaling nans x * 1 is equivalent to x,
 389 but if signaling nan is supported then the former
 390 raises an invalid exception.
 391 This may complicate things further if one wants to write
 392 portable fp math code.
 393 <p>
 394 A further libm design issue is the math_errhandling macro:
 395 it specifies the way math function errors can be checked
 396 which is either fp exceptions or the errno variable or both.
 397 The fp exception approach can be supported with enough care
 398 but errno is hard to support: certain library functions
 399 are implemented as a single asm instruction (eg sqrt),
 400 the only way to set errno is to query the fp exception flags
 401 and then set the errno variable based on that.
 402 So eventhough errno may be convenient, in libm it is
 403 not the right thing to do.
 404 <p>
 405 For soft-float targets however errno seems to be the only option
 406 (which means annex K cannot be fully supported, as it requires
 407 the support of exception flags).
 408 The problem is that at context switches the fpu status should
 409 be saved and restored which is done by the kernel on hard-fp
 410 architectures when the state is in an fpu status word.
 411 In case of soft-fp emulation this must be done by the c runtime:
 412 context switches between threads can be supported with thread local
 413 storage of the exception state, but signal handlers may do floating-point
 414 arithmetics which should not alter the fenv state.
 415 Wrapping signal handlers is not possible/difficult for various
 416 reasons and the compiler cannot know which functions will be used
 417 as signal handlers, so the c runtime has no way to guarantee that
 418 signal handlers do not alter the fenv.
 419
 420 <li>Complex arithmetics
 421 <p>
 422 gcc turns
 423 <pre>
 424 x*I
 425 </pre>
 426 into
 427 <pre>
 428 (x+0*I)*(0+I) = 0*x + x*I
 429 </pre>
 430 So if x=inf then the result is nan+inf*I instead of inf*I
 431 (an fmul instruction is generated for 0*x)
 432 <p>
 433 (There were various complex constant folding issues as well in gcc).
 434 <p>
 435 So a+b*I cannot be used to create complex numbers,
 436 instead some double[2] manipulation should be used to
 437 make sure there is no fmul in the generated code
 438 (freebsd uses a static inline cpack function)
 439 (It's not clear which is better: using union hacks
 440 everywhere in the complex code or creal, cimag, cpack)
 441 <p>
 442 Complex code is hard to make portable across compilers
 443 as there are no good compiler support, it is rarely used
 444 and there are various options available in the standard:
 445 complex support is optional, imaginary type
 446 is optional.
 447 <p>
 448 It's not clear how I (or _Complex_I) should be defined
 449 in complex.h
 450 (literal of the float imaginary unit is compiler specific,
 451 in gcc it can be 1.0fi).
 452
 453 <li>Signed int
 454 <p>
 455 The freebsd libm code has many inconsistencies
 456 (naming conventions, 0x1p0 notation vs decimal notation,..),
 457 one of them is the integer type used for bit manipulations:
 458 The bits of a double are unpacked into one of
 459 int, int32_t, uint32_t and u_int32_t
 460 integer types.
 461 <p>
 462 int32_t is used the most often which is not wrong in itself
 463 but it is used incorrectly in many places.
 464 <p>
 465 int is a bit worse because unlike int32_t it is not guaranteed
 466 to be 32bit two's complement representation. (but of course in
 467 practice they are the same)
 468 <p>
 469 The issues found so far are left shift of negative integers
 470 (undefined behaviour), right shift of negative integers
 471 (implementation defined behaviour), signed overflow
 472 (implementation defined behaviour), unsigned to signed conversion
 473 (implementation defined behaviour).
 474 <p>
 475 It is easy to avoid these issues without performance impact,
 476 but a bit of care should be taken around bit manipulations.
 477 </ul>
 478
 479 <h3><a name="implementations" href="#implementations">libm implementations</a></h3>
 480 <ul>
 481 <li>
 482 unix v7 libm by Robert Morris (rhm) around 1979<br>
 483 see: "Computer Approximations" by Hart & Cheney, 1968<br>
 484 used in: plan9, go
 485 <li>
 486 cephes by Stephen L. Moshier, 1984-1992<br>
 487 <a href="http://www.netlib.org/cephes/">http://www.netlib.org/cephes/</a><br>
 488 see: "Methods and Programs for Mathematical Functions" by Moshier 1989
 489 <li>
 490 fdlibm by K. C. Ng around 1993, updated in 1995?<br>
 491 (1993: copyright notice: "Developed at SunPro ..")<br>
 492 (1995: copyright notice: "Developed at SunSoft ..")<br>
 493 eg see: "Argument Reduction for Huge Arguments: Good to the Last Bit" by Ng 1992<br>
 494 <a href="http://www.netlib.org/fdlibm">http://www.netlib.org/fdlibm</a><br>
 495 used in: netbsd, openbsd, freebsd, bionic, musl
 496 <li>
 497 libultim by Abraham Ziv around 2001 (lgpl)<br>
 498 (also known as the ibm accurate portable math lib)<br>
 499 see: "Fast evaluation of elementary mathematical functions with correctly rounded last bit" by Ziv 1991<br>
 500 <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 used in: glibc
 502 <li>
 503 libmcr by K. C. Ng, 2004<br>
 504 (sun's correctly rounded libm)
 505 <li>
 506 mathcw by Nelson H. F. Beebe, 2010?<br>
 507 <a href="http://ftp.math.utah.edu/pub/mathcw/">http://ftp.math.utah.edu/pub/mathcw/</a><br>
 508 (no sources available?)
 509 <li>
 510 sleef by Naoki Shibata, 2010<br>
 511 (simd lib for evaluating elementary functions)<br>
 512 <a href="http://shibatch.sourceforge.net/">http://shibatch.sourceforge.net/</a>
 513 <li>
 514 crlibm by ens-lyon, 2004-2010<br>
 515 <a href="http://lipforge.ens-lyon.fr/www/crlibm/">http://lipforge.ens-lyon.fr/www/crlibm/</a><br>
 516 see: <a href="http://lipforge.ens-lyon.fr/frs/?group_id=8&amp;release_id=123">crlibm.pdf</a><br>
 517 see: "Elementary Functions" by Jean-Michel Muller 2005<br>
 518 see: "Handbook of Floating-Point Arithmetic" by (many authors from Lyon) 2009<br>
 519 <li>
 520 various other closed ones: intel libm, hp libm for itanium,..
 521 </ul>
 522
 523 <h3><a name="tests" href="#tests">libm tests</a></h3>
 524 <ul>
 525 <li><a href="http://www.netlib.org/fp/ucbtest.tgz">ucbtest.tgz</a>
 526 <li><a href="http://www.jhauser.us/arithmetic/TestFloat.html">TestFloat</a>
 527 <li><a href="http://cant.ua.ac.be/old/ieeecc754.html">ieeecc754</a>
 528 <li><a href="http://www.loria.fr/~zimmerma/mpcheck/">mpcheck</a>,
 529 <a href="http://gforge.inria.fr/projects/mpcheck/">mpcheck</a>
 530 <li><a href="http://people.inf.ethz.ch/gonnet/FPAccuracy/Analysis.html">FPAccuracy</a>
 531 <li><a href="http://www.math.utah.edu/~beebe/software/ieee/">beebe's ieee tests</a>
 532 <li><a href="http://www.math.utah.edu/pub/elefunt/">elefunt</a>
 533 <li><a href="http://www.vinc17.org/research/testlibm/index.en.html">testlibm</a>
 534 <li><a href="http://www.vinc17.org/research/fptest.en.html">fptest</a>
 535 <li>tests in crlibm
 536 </ul>
 537 <p>
 538 multiprecision libs (useful for tests)
 539 <ul>
 540 <li><a href="http://mpfr.org/">mpfr</a>
 541 <li><a href="http://www.multiprecision.org/index.php?prog=mpc">mpc</a>
 542 <li><a href="http://pari.math.u-bordeaux.fr/">pari</a>
 543 <li><a href="http://www.apfloat.org/apfloat/">apfloat</a>
 544 <li><a href="http://code.google.com/p/fastfunlib/">fastfunlib</a>
 545 <li>scs_lib in crlibm
 546 </ul>
 547
 548 <p>other links
 549 <ul>
 550 <li>ieee standard: http://754r.ucbtest.org/
 551 <li>extended precision issues: http://www.vinc17.org/research/extended.en.html
 552 <li>correctly rounded mult: http://perso.ens-lyon.fr/jean-michel.muller/MultConstant.html
 553 <li>table based sincos: http://perso.ens-lyon.fr/damien.stehle/IMPROVEDGAL.html
 554 <li>finding worst cases: http://perso.ens-lyon.fr/damien.stehle/WCLR.html
 555 <li>finding worst cases: http://perso.ens-lyon.fr/damien.stehle/DECIMALEXP.html
 556 <li>finding worst cases: http://www.loria.fr/equipes/spaces/slz.en.html
 557 <li>finding worst cases: http://perso.ens-lyon.fr/jean-michel.muller/Intro-to-TMD.htm
 558 <li>generating libm functions: http://lipforge.ens-lyon.fr/www/metalibm/
 559 <li>fast conversion to fixedpoint: http://stereopsis.com/sree/fpu2006.html
 560 <li>double-double, quad-double arithmetics: http://crd-legacy.lbl.gov/~dhbailey/mpdist/
 561 <li>papers by kahan: http://www.cs.berkeley.edu/~wkahan/
 562 <li>fp paper by goldberg: http://download.oracle.com/docs/cd/E19957-01/806-3568/ncg_goldberg.html
 563 <li>math functions in general: http://dlmf.nist.gov/
 564 </ul>
 565
 566 <p><small>(<a href="/git/?p=www">page history</a>, <a href="/">home</a>)</small>
 567 </body></html>