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

1. Objectives

The goals of the project are to:

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

  • Write tools for automatic code generation (i.e., automatic generation of known to be good {x, …, f(x, …)} pairs).

  • 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 now, i.e. not FPU optimized).

  • Write tools that will allow us 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.

2. Manually run the tests

In order to checkout the tests you need Git:

git clone https://ekamperi@bitbucket.org/ekamperi/mathlib.git

There’s also a web 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 or don’t want to afford git.

Subsequent updates are done with:

% cd mathlib
mathlib% git pull

In order for the tests to be ran, you need the Automated Testing Framework which NetBSD ships by default. If you don’t run NetBSD and you have installed ATF by yourself in some non-default prefix, e.g. /usr/local, make sure that your PATH environmental variable is properly set. The test suite has been rebased on top of ATF 0.12, but older versions may still work.

You also need gmake, autoconf and automake. In order to generate beautiful html reports you have to install libxslt. All of these packages, are provided by pkgsrc. We are not happy with these dependencies, but they allowed us to be portable and effective, in the limited timeframe of the Google Summer of Code event.

To build, run the tests and generate the html report:

% cd mathlib
mathlib% make all
[...]
mathlib% make run
[...]
mathilb% make run-html
[...]
Note The math/, complex/ and fenv/ subdirectories have an Atffile in them. There you can tweak the iterations parameter. This variable affects the stochastic tests inside the suite: the larger the value, the more exhaustive the tests and the more time they consume to complete.

3. Test run logs

We are running irregularly irregular automated tests against:

  • NetBSD

    • x86, amd64, vax, m68k and sparc64

  • FreeBSD

    • x86, amd64

  • Linux

    • x86, amd64

  • OpenSolaris

    • amd64

  • Darwin

    • 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.

4. Unified diffs

Here you can find (more or less) up to date diffs against NetBSD HEAD.

Note The patches regarding fenv.h support for x86, amd64 and for sparc64, have been committed to the official NetBSD source repository. Some long double versions and general cleanup here and here. The remquo{,f}() functions have been committed.

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=/'

5. 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 (the link has been broken). 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 ever-lasting loops on certain input ranges. If you experience such behaviour, skip that function until it is fixed upstream or worked around downstream.

6. Performance measurement

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.

7. 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.

8. Status | Todo

8.1. 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

8.2. fenv.h - data types

i387 amd64 sparc64 m68k

fexcept_t

YES

YES

YES

WIP

fenv_t

YES

YES

YES

WIP

8.3. 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. Currently we don’t.

Errno is already available. As for floating-point exceptions, currently only i386, amd64 and sparc64 support them. Question though is whether errno and/or floating-point evnironment are properly updated by libm during the various computations.

To shed some light on this, we forced the definition of MATH_ERRNO and MATH_ERREXCEPT, as if we did support them, and ran the test suite. Preliminary results for i386 errno and errexcept.

8.4. 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.

8.5. 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.

8.6. VAX

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

8.7. Build issues in tests (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/.

8.8. Notable references

9. 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.