Note
 Portions of stuff mentioned here are WIP. If this page doesn’t render properly, blame asciidoc.

Objectives

The goals of the project are to:

  • Write unit tests for the math.h, fenv.h, complex.h and float.h interfaces.

  • Implement support for floating-point exceptions for common architectures (x86/amd64) and for less common (sparc64, m68k).

  • Implement long double versions of math functions (software emulated for a start, i.e. not FPU optimized).

  • Write tools that will allow observability of the mathlib, such as measuring the error in floating-point computations or the speed of the algorithms.

  • Audit existing code, fix or document where due.

How to get you started

ATF tests

In order to checkout the tests you need Git:

git clone git://gitweb.dragonflybsd.org/~beket/mathlib.git

Subsequent updates are done with:

% cd mathlib
mathlib% git pull

There’s also a Gitweb interface for checking the code with your browser. From the same place you can download tar archives of the entire source tree, if for some reason you can’t/don’t want to afford git.

In order to run the tests you need Automated Testing Framework. NetBSD ships it by default. If you don’t run NetBSD and you have installed ATF by yourself with some non-default prefix, e.g. /usr/local, make sure that your PATH environmental variable is set properly.

You also need gmake. I’m not happy with it, but it was the only viable option to write portable Makefiles.

% cd mathlib
mathlib% gmake all
[...]
mathlib% gmake run
[...]
mathlib%

mathlib unified diffs

Here you can find (more or less) up to date diffs against NetBSD HEAD. The patches about fenv.h interface support have been committed to the official NetBSD CVS repository. We are now working on fenv.h for sparc64, long double versions for IEEE FP archs, namespace issues or comments to existing code.

You can apply the diffs with:

% patch -p0 < mathlib.diff

And then refer to NetBSD guide on how to build source and install new world. Or just follow these steps:

% su root -c 'mkdir /usr/obj; chmod g+w /usr/obj'
% cd /usr/src
% ./build.sh -O /usr/obj -T /usr/tools -U distribution sets
% su root -c 'cd /usr/src && ./build.sh -O /usr/obj -U install=/'

Error estimation of libm functions

Here you can find a utility to calculate the error in terms of ULPs for most of the libm functions. This includes double and long double precision, real and complex valued functions. The utility depends on the existence of gmp, mpfr and mpc libraries. gmp is an arbitrary precision arithmetic library, mpfr a multi-precision floating-point library and mpc same as mpfr but for complex numbers.

In order to build it, please do the following:

mathlib% cd ulps
mathlib/ulps% make

Only autoconf is needed (2.59 or later). Automake is not needed.

Example usage:

mathlib/ulps% ./ulps
ulps: [all | list | <function1> <function2> ...]
mathlib/ulps% ./ulps all
        FUNCTION     Max ULP    Min ULP    Avg ULP    skipped
[ 1/62] acos         0.5000     0.0000     0.0000           0
        acosl        0.5000     0.0000     0.0000           0
[ 2/62] acosh        0.5000     0.0000     0.0001           0
        acoshl       0.5000     0.0000     0.1638           0
[ 3/62] asin         0.0000     0.0000     0.0000           0
        asinl        0.5000     0.0000     0.0000           0
[ 4/62] asinh        0.5000     0.0000     0.0001           0
        asinhl       1.0000     0.0000     0.0824           0
[ 5/62] atan         0.0000     0.0000     0.0000           0
        atanl        0.5000     0.0000     0.0000           0
Percentage complete: 83.60
^C

Or individual function names may be provided, e.g.

mathlib/ulps% ./ulps cbrt fabs
        FUNCTION     Max ULP    Min ULP    Avg ULP    skipped
[ 1/ 2] cbrt         0.5000     0.0000     0.0909           0
[ 2/ 2] fabs         0.0000     0.0000     0.0000           0
mathlib/ulps%

The definition of ULP that we are using may be found in this document. Although reasonable voices have been raised regarding the fitness of this definition, its straight-forwardness from an implementation point of view and the fact that it allows us direct comparisons with good libm’s, like Solaris, made us stick to it for the moment.

If you have some exotic hardware, please do me a favor. Run the diagnostics and mail me back the results. I have uploaded mine here. More to come.

Note
 Incidentally, we discovered a bug in mpfr. If mpfr_tgamma() was called with an argument that would naturally cause an underflow, it would enter an infinite loop. Upstream developers were informed and problem was solved both in svn trunk and in 3.0 branch. See also this. Also, we discovered a few bugs in mpc library, affecting most notably the complex trigonometric functions. The symptoms are infinite loops on certain input ranges. If you experience such behaviour, skip that function until it is fixed upstream (we are wokring on it).

Performance measure

We have performance probes for fenv.h interface, e.g.

mathlib/etc% ./proffenv
feclearexcept() FE_ALL_EXCEPT:   256 msecs for 1000000 iterations
feclearexcept()        random:   289 msecs for 1000000 iterations
     fegetenv()              :   233 msecs for 1000000 iterations
feraiseexcept()              :   457 msecs for 1000000 iterations
     fesetenv()    FE_DFL_ENV:   241 msecs for 1000000 iterations
     fesetenv()        random:   483 msecs for 1000000 iterations
  feupdateenv()    FE_DFL_ENV:   291 msecs for 1000000 iterations
  feupdateenv()        random:   533 msecs for 1000000 iterations
mathlib/etc%

Since recently, we are able to measure the performance of math.h and complex.h function, with respect to their input. Code is here. Graphs are here. You can generate your own, if you have gnuplot installed, via:

% cd mathlib/etc
mathlib/etc% make gen-graphs
[...]

Most algorithms exhibit constant-time behaviour, but there are some interesting cases, e.g., compare NetBSD cos() against, OpenSolaris cos(). (Our graph maybe interpreted on the basis of iterative argument reduction for very large inputs. It seems that Sun’s developers managed to short circuit it.)

This tool will allow us to identify pathological cases in existing implementations and protect us from introducing such when we will implement the long double versions.

C code generator

For every math function, we test it against a small set of good {x, …, f(x, …)} pairs. These data are generated automatically by a utility, that resides here. It uses the gmp/mpfr libraries. The output is valid C code that can be used without any modifications. This is an example.

Automated tests

We are running irregularly irregular automated tests against:

  • NetBSD

    • x86, amd64, vax, m68k and sparc64

  • FreeBSD

    • x86, amd64

  • Linux

    • x86, amd64

  • OpenSolaris

    • amd64

If you have some {OS, platform} combination that is not listed in here, I would be grateful if you did:

~% cd mathlib
mathlib% make run-html

And mailed me back the results.html file.

Status | Todo

fenv.h - functions

i387 amd64 sparc64 m68k

feupdateenv

YES

YES

YES

WIP

feclearexcept

YES

YES

YES

WIP

fegetexceptflag

YES

YES

YES

WIP

feraiseexcept

YES

YES

YES

WIP

fesetexceptflag

YES

YES

YES

WIP

fetestexcept

YES

YES

YES

WIP

fegetround

YES

YES

YES

WIP

fesetround

YES

YES

YES

WIP

fegetenv

YES

YES

YES

WIP

feholdexcept

YES

YES

YES

WIP

fesetenv

YES

YES

YES

WIP

fenv.h - data types

i387 amd64 sparc64 m68k

fexcept_t

YES

YES

YES

WIP

fenv_t

YES

YES

YES

WIP
Note
 sparc64 fenv.h interface has been written on a x86 host (cross-compiled). Therefore, the only guarantee at the moment is that it builds fine. If you have a sparc64 machine, please consider applying the patch and running the fenv/ test suite.

Error signaling

There are two ways for libm to signal errors during computations. Either via errno or by raising floating-point exceptions. POSIX expects that at least one is supported. We should expose our capabilities by defining math_errhandling and MATH_ERRNO|MATH_ERREXCEPT.

Errno is readily available. As for floating-point exceptions, currently only i386 and amd64 support them. Question is whether errno and/or floatint-point evnironment is properly manipulated by libm.

Preliminary results for i386 errno and errexcept.

Namespace issues

  • Defining _POSIX_C_SOURCE or _XOPEN_SOURCE prior to the inclusion of math.h hides copysign() function. Possibly others too.

  • signbit() uses a gcc builtin function that only works with floats and doubles. Roll out our own just like FreeBSD does.

Precision loss of periodic functions in x86

If you take a look at the ulp reports you will notice that, except for OpenSolaris, all libm’s perform very poorly when it comes to sin or cos functions (and possibly others). The cause appears to be the limited precision of Pi (66 bits), that is used during the argument reduction process.

This topic is further discussed in this thread of Intel’s forum.

VAX

We don’t have nextafter() or nextbefore() implementations for VAX, which blocks us from calculating the ULPs of the various math functions.

Build issues in tests (almost done)

For every function we write test cases for, we check its float, double and long double version, in the same test program. On the other hand, NetBSD misses a lot of *l() functions, which blocks test programs from building. Therefore, we need to use autoconf and surround the hot spots with #ifdef blocks. Just like we do in ulps/.

Notable references

Contact

For any questions and/or suggestions, please feel free to contact me via e-mail. My address is ekamperi at gmail dot com. You can also find me via IRC in Freenode network at #netbsd-code. My nick is Beket.