High Precision Arithmetic Library

by Ivano Primi <ivprimi (a) libero.it>
Last Update: 2005-08-15

1. News

• 2005-08-12: HPAlib 1.6 is finally ready !

The High Precision Arithmetic (HPA) library implements a high precision floating point arithmetic together with a comprehensive set of support functions. The general areas covered by these functions include:

• Extended Precision Arithmetic
• Extended Precision Math Library
• Applications of High Precision Computation

The math library support includes evaluation of trigonometric, inverse trigonometric, hyperbolic, inverse hyperbolic, logarithm, and exponential functions at the same precision as the floating point math itself. The HPA library also supports high precision complex arithmetic and includes an Extended Precision Complex Math Library.

The last version (1.6) of the HPA library comes with a C++ wrapper, allowing to perform high precision computations with the same syntax of the normal code.

3. What does high precision mean ?

The C and C++ languages define three types of floating point numbers: a floating point number can be of type `float`, `double` or `long double`. When a programmer wants that a variable can store not only integer values but also values with a non-null fractional part, then he can define this variable as a `float`, `double` or `long double` variable. If you are still reading this page, probably you know about the usual binary representation of the numbers in the memory and in the CPU of a computer. You know that it is possible to represent only a limited range of numbers and that non integer numbers are often approximated. Actually, a non integer number could also have a non periodic and infinitely long fractional part, which has to be truncated to fit into the memory of a computer: even if today the computers have a big amount of memory, it is still finite ! Probably you know that the difference between the three kinds of floating point data types above introduced consists in the different `range` of representable numbers and in the accuracy level which these numbers are approximated with when they can not be exactly represented in binary notation. Likely you also know that each floating point data type has an associated data size, where size means the number of bytes in the memory of the computer which are needed to keep the current value of a variable of such data type. According to the IEEE754 standard, the size of a `float` variable is 4 bytes, the size of a `double` variable is 8 bytes and the size of a `long double` variable is at least 10 bytes, even if the ANSI rules for C and C++ only state that the size of a `long double` variable must be not less than that one of a `double` variable. Naturally, in C/C++ one can recover the exact size of each floating point data type through the `sizeof` operator. The general rule about the size of a floating point data type is that greater size implies a larger range of representable numbers and higher accuracy in the internal representation of the numbers.

Roughly speaking, the HPA library introduces a new floating point data type, which has a greater size (at least 16 bytes) and then allows to represent a larger range of numbers with an higher accuracy than the standard floating point data types.

The normal configuration of the library defines an extended floating point data type having a size of 16 bytes, which are so split:

• 1 bit for the sign (0 -> positive, 1 -> negative),
• 15 bits for the (biased) exponent,
• 112 bits (= 7 words of 16 bits length) of mantissa, corresponding about to 32 decimal digits of accuracy.

The range of representable numbers is then given by

```                    2^16384 > x > 2^[-16383]  or
```

or

```               1.19*10^4932 > x > 1.68*10^-[4932].
```

The HPA library also defines an extended precision complex type, used to manage computations involving complex numbers. From the point view of the HPA library, a complex number is simply a structure formed by two extended floating point numbers, representing respectively the real and the imaginary part of the complex number.

The accuracy provided by the extended floating point data type of the HPA library can be controlled by adjusting the size of its mantissa before compiling and installing the library. Normally, the size of the mantissa is equal to 7 words of 16 bits length, as told above. But it can be incremented by suitably setting the macro XDIM. XDIM defines the size of the mantissa of an extended floating point value in terms of words of 16 bits length. XDIM can take one of the following values: 7, 11, 15, 19, 23, 27 and 31. For instance, XDIM = 31 implies about 149 decimal digits of accuracy.

However, once the HPA library has been built and installed, the size of its extended floating point data type can not be changed, unless you rebuild and reinstall the library. So, the HPA library does not supply the dynamic resizing of the numbers according to the precision requested time by time by the current application and the actual number of binary digits needed to represent a given value (for instance, a number like 2.0 or 4.0 or 8.0 only needs one digit of mantissa to be represented !). If you are looking for similar features, then what you really want is an arbitrary precision arithmetic library and so you have to look elsewhere, for instance here.

4. Good qualities and weaknesses

The functions forming the HPA library are all implemented in a portable fashion in the C language. For values of XDIM not too large, they are enough quick. Even if there are no benchmarks, on an Intel Pentium(tm) I, 200 MHz, with 32 Mb of RAM (Operating system: NetBSD 2.0.2, compiler: gcc 3.3.3, optimization flags: -O3) the suite of test programs distributed together with the library runs in a few of seconds when XDIM = 7. For greater values of XDIM (>= 23 ) the overall performance of the routines of the library gets rapidly worse. For instance, when XDIM = 31 and on the same machine as before, the execution time for all the test programs is about 3 minutes !

Unfortunately, high precision arithmetic becomes too expensive when the mantissa of the extended floating point data type is very long. So, if you need an accuracy level much greater than the quadruple precision (XDIM = 7), then the arbitrary precision approach turns out more convenient. Again, if this is the case, take a look at http://swox.com/gmp/.

The HPA library (HPAlib) and its C++ wrapper have been extensively and deeply tested. At the moment there are only a few of known weaknesses. They are listed in the BUGS file distributed together with the source code of the library. I think that the routines of the library are quite reliable even if I can NOT give you any warranty about that ! Moreover, the Extended Precision Math Library is complete (or almost complete) both in its real and in its complex module: it includes trigonometric, inverse trigonometric, hyperbolic, inverse hyperbolic, logarithm, and exponential functions .

The documentation of the library is really detailed, even if there could be several grammatical or lexical errors: I am not a native English speaker !

The C++ wrapper makes very simple to pass from double precision computations to extended precision ones. In most cases, the unique changes you will have to make to your C++ source code will be:

1. Include the header file `<xreal.h>` (`<xcomplex.h>`) at the head of your source file(s),
2. Add the directive `using namespace HPA;`,
3. Replace the `double` (`complex`) data type with `xreal` (`xcomplex`).

For the future releases of the HPAlib I plan to speed up the execution times of the main routines of the library. Moreover, although the IEEE 754 standard for floating point hardware and software is assumed in the PC/Unix version of this library, the HPAlib is not fully compliant to the IEEE 754 standard in what concerns the implementation of the quadruple precision. The compliance to IEEE 754 could be one of the main purposes to get with the future releases of the library.

The HPA library comes from a branch of the source code of the CCMath library, which is a work by Daniel A. Atkinson <DanAtk (a) aol com>. Daniel A. Atkinson was very kind to release the code of the CCMath Library under GNU Lesser General Public License. This made possible that Ivano Primi could modify, complete and redistribute this source code under the same terms.

So, the HPA (abbreviation of High Precision Arithmetic) Library is copyrighted by Ivano Primi <ivprimi (a) libero it> and Daniel A. Atkinson (the author of the original code). As it is for the source code of the CCMath Library, the source code of the HPA library is released under the terms of the GNU Lesser General Public License, as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version.

The HPA library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.

You can contact me by paper mail writing to the next address

Ivano Primi via Colle Cannetacce 50/A C.A.P. 00038 - Valmontone (ROMA) Italy .

If you prefer the electronic mail you can write to the address

<ivprimi (a) libero it> .

The tar-gzipped archive with the source code of the HPA library can be downloaded from

The last stable release of the HPA library is given by the version 1.6 . Together with the source code, the archive contains the full documentation (in English) about the HPA library. The manual of the HPA library is available in the following formats:

• HTML,
• LaTeX,
• PDF,
• plain ASCII text file.

Permission is granted to copy, distribute and/or modify these documents under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of the license is always included in the section entitled "GNU Free Documentation License". You can also obtain a copy of the GNU Free Documentation License by writing to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.

The manual of the HPA library can also be browsed online here.

7. Acknowledgments

Firstly I feel like to thank Daniel A. Atkinson, since without its work the HPA library would not exist. Next I have to thank Aurelio Marinho Jargas (<verde (a) aurelio net>), author of txt2tags (http://txt2tags.sf.net), a free (GPL'ed) and wonderful text formatting and conversion tool, which I extensively used in writing the manual of the HPA library and this web page too !

Last but not least, I want to thank all the people till now involved in the Free Software community, starting from those ones directly involved in the GNU project (http://www.gnu.org). Without their great work, this little one would never be done.