libm

This page is about libm for the musl libc.

Sources

Writing math code from scratch is a huge work so already existing code is used. Several math functions are taken from the freebsd libm and a few from the openbsd libm implementations. Both of them are based on fdlibm. The freebsd libm seems to be the most well maintained and most correct version of fdlibm.

sources:

General rules

Representation

Binary representation of floating point numbers matter because bit hacks are often needed in the math code. (in particular bit hacks are used instead of relational operations for nan and sign checks becuase relational operators raise invalid fp exception on nan and they treat -0.0 and +0.0 equally and more often than not these are not desired)

float and double bit manipulation can be handled in a portable way in c using union types:

(assuming the bits in the object representation of 32bit and 64bit unsigned ints map to the floating-point representation according to ieee-754, this is not always the case, eg. old arm floating-point accelerator (FPA) used mixed endian double representation, but musl does not support the old arm ABI)

long double bit manipulation is harder as there are various representations and some of them don't map to any unsigned integer type:

(and other non-standard formats are not supported)

In case of ld64 the bit manipulation is the same as with double and all long double math functions can be just wrappers around the corresponding double ones. (using symbol aliasing on the linker level is non-conformant since functions would not have unique address then)

ld80 is the most common long double on linux (i386 and x86_64 abi), it means 64bit significand with explicit msb (inconsistent with other ieee formats), 15bit exp, 1 sign bit. The m68k (and m88k) architecture uses the same format, but different endianness:

where m is the significand and se is the sign and exponent.

ld128 is rare (eg. sparc64 with software emulation), it means 113bit significand with implicit msb, 15bit exp, 1 sign bit:

There are other non-conformant long double types: eg. the old SVR4 abi for ppc uses 128 bit long doubles, but it's software emulated and traditionally implemented using two doubles (also called ibm long double as this is what ibm aix used on ppc). The ibm s390 supports the ieee 754-2008 compliant binary128 floating-point format, but previous ibm machines (S/370, S/360) used slightly different representation.

This variation shows the difficulty to consistently handle long double: the solution is to use ifdefs based on float.h and on the endianness and write different code for different architectures.

Ugly

The ugly parts of libm hacking.

Some notes are from: http://www.vinc17.org/research/extended.en.html

Useful info about floating-point in gcc: http://gcc.gnu.org/wiki/FloatingPointMath

libm implementations

libm tests

multiprecision libs (useful for tests)

other links

(page history, home)