MenuBar/IconUsers.png "IconUsers.png")
MenuBar/gsoc_2010.png "gsoc_2010.png")
MenuBar/IconGlobe.png "IconGlobe.png")
Accurate and Portable Elementary Functions (A.P.E.F.)
Contents
Caution
This page is obsolete. This project is not part of the Scilab list of ideas anymore, since it has already been done in 2011.
Introduction
Scilab provides several elementary functions such as pow (i.e. x^y), sin, tan, exp, log, etc... These functions are provided by the mathematical library associated with the compiler. More specifically:
- on Windows, the elementary functions are provided by the Microsoft Visual C compiler,
- on Linux, the elementary functions are provided by GCC, the GNU Compiler Collection.
This implies that the results provided by the sine function, for example, are potentially not the same for Windows and for Linux. Still, for most common computations, the results are exactly the same. But, for some inputs x and some functions, the results are drastically different. This is especially the case when the function is ill-conditionned for some x.
GSOC 2011
Some work has been done in the context of the GSoC 2011. The work has been published in a dedicated git branch: http://cgit.scilab.org/scilab/log/?h=apef
See the list of reports:
To sum up, the main problem of the final result is the low performances of the library. In the context of the GSoC 2012, performance could be improved.
Some examples
In this section, we describe some examples of bugs which were detected with respect to accuracy of the elementary functions.
Bug #4048: The pow function
This bug is [http://bugzilla.scilab.org/show_bug.cgi?id=4048 #4048].
On Linux 32 bits, we get :
-->0.9999999999999999^-18014398509482000.0
ans =
54.598148 On Windows 32 bits, we get:
-->0.9999999999999999^-18014398509482000.0
ans =
7.3890561 The floating point representation is
fl(0.9999999999999999)= 9007199254740991 * 2^-53
And Wolfram Alpha produces :
7.389056098930664173180786636157853718020054020964384648...
The condition number for this case is 1.8D+16, which is a bad as can be.
It is due to a bug in the glibc on double precision: http://sources.redhat.com/bugzilla/show_bug.cgi?id=706
Bug #8617: The cos function is not accurate for some x on Linux 32 bits
This is the bug [http://bugzilla.scilab.org/show_bug.cgi?id=8617 #8617].
The test which fails in cos.tst is the comparison of the exact and computed value of cos(x), for x=1.570796326794896558, which is a 18 digits decimal representation of pi/2. The failure of this test indicates that the computed value is off by more than 1 ulp from the expected result. Actually, the error is much worse than a few ulps.
With Scilab 5.3.0 on Linux 32 bits Ubuntu 10.04:
-->format("e",25)
-->cos(1.570796326794896558D+00)
ans =
6.123031769111886291D-17 With Scilab 5.3.0 on Windows XP 32 bits:
-->format("e",25)
-->cos(1.570796326794896558D+00)
ans =
6.123233995736766036D-17 The exact value of cos(%pi/2) is :
cos(7074237752028440 * 2^(0-53+1))= 6.123233995736765886130329661375001464640377798836283... × 10^-17
as given by Wolfram Alpha. Hence, the relative error on each of these platforms is given by:
-->expected = 6.123233995736765886D-17; -->abs(cos(1.570796326794896558)-expected)/expected
which produces :
- Scilab 5.3.0 Linux 32 bits Ubuntu 10.04 : 3.3D-05
- Scilab 5.3.0 Windows XP 32 bits : 0
This shows that the cos function provided by Visual on Windows is exact on x=1.570796326794896558, while the cos function provided by gcc on this particular Linux 32 bits only has ~5 significant digits, which is not good at all. We may be satisfied by 15 digits, or even the exact result, but not less. The argument reduction algorithm in cos cannot justify such a gross innacuracy, since the value x=1.570796326794896558 is not large (the point x=1.e16 would be more questionable).
What do we need
The ultimate solution to fix these issues is to use a library which is correctly rounded for all inputs x. This is still a research topic for many functions, but is solved for some specific functions. On the other hand, if we were using the same library on all platforms, the computations should (at least partially) produce the same results on all platforms. The reason why this will solve only partially the problem is because compiling options or processor's specific features can lead to different floating point results, for example caused by the use of extended precision.
In the two next sections, we suggest libraries which could be connected to Scilab.
Fdlibm
Developed at Sun Microsystems, Inc.
FDLIBM (Freely Distributable LIBM) is a C math library for machines that support IEEE 754 floating-point arithmetic. In this release, only double precision is supported.
FDLIBM is intended to provide a reasonably portable (see assumptions below), reference quality (below one ulp for major functions like sin,cos,exp,log) math library (libm.a). For a copy of FDLIBM, please see
or
Fdlibm provides the following functions:
- double acos (double);
- double asin (double);
- double atan (double);
- double atan2 (double, double);
- double cos (double);
- double sin (double);
- double tan (double);
- double cosh (double);
- double sinh (double);
- double tanh (double);
- double exp (double);
- double frexp (double, int *);
- double ldexp (double, int);
- double log (double);
- double log10 (double);
- double modf (double, double *);
- double pow (double, double);
- double sqrt (double);
- double ceil (double);
- double fabs (double);
- double floor (double);
- double fmod (double, double);
- double erf (double);
- double erfc (double);
- double gamma (double);
- double hypot (double, double);
- int isnan (double);
- int finite (double);
- double j0 (double);
- double j1 (double);
- double jn (int, double);
- double lgamma (double);
- double y0 (double);
- double y1 (double);
- double yn (int, double);
- double acosh (double);
- double asinh (double);
- double atanh (double);
- double cbrt (double);
- double logb (double);
- double nextafter (double, double);
- double remainder (double, double);
- double scalb (double, int);
- int matherr (struct exception *);
- double significand (double);
- double copysign (double, double);
- int ilogb (double);
- double rint (double);
- double scalbn (double, int);
- double expm1 (double);
- double log1p (double);
Libmcr
http://www.sun.com/download/products.xml?id=41797765
A reference correctly-rounded library of basic double-precision transcendental elementary functions.
libmcr is offered as a reference library for selected double-precision transcendental elementary functions. The library is supplied as source code allowing the researcher to build it for a variety of machines.
The code is distributed under Sun's Public Source Code License. All library source code is distributed with the license detailing use and redistribution limitations.
The code was developed and is distributed and maintained by Sun's Floating Point Libraries group.
Included Functions:
- double exp(double);
- double log(double);
- double pow(double,double);
- double sin(double);
- double cos(double);
- double tan(double);
- double atan(double);
Crlibm
The purpose of this project is to offer a mathematical library (libm) with proven, IEEE-754 compliant, correct rounding in the four rounding modes, and performances comparable to standard libms.
CRlibm, an efficient and proven correctly-rounded mathematical library
CRlibm is a free mathematical library (libm) which provides:
- implementations of the double-precision C99 standard elementary functions,
- correctly rounded in the four IEEE-754 rounding modes,
- with a comprehensive proof of both the algorithms used and their implementation,
- sufficiently efficient in average time, worst-case time, and memory consumption to replace existing libms transparently.
CRlibm is distributed under the GNU Lesser General Public License (LGPL). Included in the distribution is an extensive documentation with the proof of each function (currently more than 100 pages), as well as all the Maple scripts used to develop the functions. This makes this library an excellent tutorial on software elementary function development.
The CRlibm library also includes a lightweight library for multiple precision, scslib (Software Carry Save Library). This library has been developed specifically to answer the needs of the CRlibm project: precision up to a few hundred bits, portability, compatibility with IEEE floating-point, performance comparable to or better than GMP, small footprint. It uses a data-structure which allows to avoid carry propagations during multiple-precision multiplications. Supported operations are essentially addition/subtraction and multiplication, and conversions. This library is independent from CRlibm.
The authors of the Crlibm project are
- Jean Michel Muller
- Florent de Dinechin
- Christoph Lauter
- Catherine Daramy Loirat
- Alex Bernier
- Matthieu Gallet
- David Defour
- Nicolas Gast
References
[1] MATLAB is a registered trademark of The MathWorks.
- [3] Visual C++ is a registered trademark of Microsoft Corporation.
[4] "Crlibm: l'évaluation des fonctions élémentaires avec arrondi correct", David Defour, 2005, http://lipforge.ens-lyon.fr/www/crlibm/documents/scilab.pdf
[5] "Computing with Floating Point. It's not Dark Magic, it's Science", de Dinechin, 2004.99999, http://lipforge.ens-lyon.fr/www/crlibm/documents/cern.pdf
Author
DIGITEO - Michaël Baudin