GNU/Linux Decimal Floating Point with GCC, libdfp, and Glibc
Digital computers were first implemented to solve problems in science and engineering. In this scenario, binary floating point (FP), with its conventional "double" data type containing 53-bits of binary precision or 16-bits of decimal precision, is sufficient for high computational accuracy. Most physical constants are not known to more than 4-5 decimal places and thus 16-decimal-place precision allows the majority of algorithms to experience an acceptable level of round-off or other errors.
It may be of note to indicate that the 53-bits of binary floating point precision was chosen by William Cahan, the father of the IEEE754 standard, because it was a value that would permit programmers to create algorithms without having to hire expensive error analysts. The IEEE754 standard indeed does allow most of us to write math routines without concern for the complex details of error prediction.
However, financial mathematics is a different ballgame. Significant figures can become very large and even exceed the 16-decimal-digit limit of conventional floating point. Moreover, most decimal numbers are not exactly representable in binary FP. The common example is $0.33. In fact, any fractional number that is not expressible in powers of 2, the binary base, cannot be exactly represented in binary FP. These issues can lead to unacceptable round off errors in financial calculations.
To circumvent these problems, standards for decimal floating point numbers were developed. Unfortunately, whereas binary FP is handled in hardware for virtually all processors, only high-end mainframes, such as IBM z-series, contain hardware for decimal FP. As a consequence, on desktop workstations or even server farms using common processors, decimal FP calculations must be done with software. This means much slower execution but due to the nature of financial versus scientific computation a loss of speed is usually not a hindrance.
GNU/Linux has the ability to implement decimal floating point operations and programmers should be encouraged to write any financial software by using decimal FP. Indeed, many institutions and even governments mandate the use of decimal FP for financial calculations.
Basic decimal FP was introduced into GNU/Linux with GCC-4.3 via the libdecnumber package. To obtain complete decimal FP mathematics an external library, libdfp, is required, along with a version of glibc above 2.10. libdfp uses GCC libdecnumber as a backend but also brings in the standard C mathematics routines and other functions (see below) in decimal floating point format.
I will now outline how to install libdfp, which was generously donated to the Free Software Foundation by IBM, and how to create C programs using both gcc and libdfp.
Decimal floating point is implemented with three new C data types:
_Decimal32 _Decimal64 _Decimal128
The details of the decimal FP format is beyond the scope of this article but can be found in many places, such as these links:
However, I will say that there are two major storage types for decimal FP: the BID (binary integer decimal) and the DPD (densely packed decimal). Only the BID type is relevant on Intel-based machines and is the default type for x86_64 GCC. Hence, only the BID type will be discussed here. The DPD type is exclusive to hardware-based decimal FP on IBM mainframes.
GCC, via the built-in accessory libgcc library, provides basic decimal FP arithmetic, comparison, and conversion operations. These are outlined here:
Although such basic operations may be sufficient to create useful financial software, there are some major limitations.
First, there are no direct conversion routines from numeric data in string format to the decimal FP types. Only indirect conversions, from the standard C binary types, such as int, long int, and double, to decimal FP, are available. Considering the error problems mentioned in the beginning of this article, these conversions may be unacceptable.
Second, there are no operations for direct output of decimal FP values, such as with fprintf(). Again, only indirect output by first converting decimal FP to standard C binary types is available.
Third, only basic arithmetic operations are available. The more complex math functions which are needed for many financial or scientific routines are unavailable.
To circumvent these limitations (if it is deemed necessary) an external library is needed. GNU/Linux has, courtesy of IBM, the libdfp package:
libdfp provides direct conversion functions from string to decimal FP, and, by using glibc "hooks," also provides direct output of decimal FP values using fprintf().
To view the many decimal FP functions that are available in libdfp just examine the ieee754r and the stdlib directories within the libdfp source code. These functions mirror those of glibc and include primarily advanced math and string/wide-char conversions. Unfortunately, at this time, no other listing or manuals of these functions are given. However, all of these functions are formatted identically, are named identically, and behave identically to their standard C function counterparts, with the exception, of course, that they operate on decimal FP data types.
Together with the libgcc arithmetic, conversion, and comparison functions, these additional libdfp functions allow the decimal FP format to be fully incorporated into C programs for any purpose.
What follows will describe how to install libdfp from source code and then how to use libdfp routines in conjunction with the built-in GCC types and operations.
Building and Installing Libdfp
Of course, many GNU/Linux distros may supply a pre-compiled version of libdfp and any user is free to use it. But some users may desire to build it from source for themselves.
Building libdfp is straightforward and follows the standard GNU/Linux "configure, make, and make install" sequence.
The first step is to download the latest source from the GitHub repository and unzip to a directory:
cd /some/directory; unzip /path/to/libdfp-master.zip
Next, the source code is configured with the desired options. A complete configure command would be something like the following. See the output of "./configure --help" for all available options.
$CFLAGS $LDFLAGS ./configure --prefix=/usr/local --with-cpu=haswell
An important configure option is "--with-cpu" as this allows the software to be optimized by using the instruction sets specific to a particular processor, i.e. SSE2, SSE4.2, AVX, etc. Not all processors may be supported but the builder should check the output of configure to see what has been selected.
Note that for configure to succeed, gcc must have been compiled with the "--enable-decimal-float" option. I would assume that most recent GNU/Linux distros have done this but I cannot be sure. The DIY builder should ensure that this is the case.
Note also that the CFLAGS and LDFLAGS variables should be set with the desired compile and linking flags.
After the configuration, just build and install:
make make check make install
Using Libdfp and GCC
Once libdfp is installed a user can create programs using decimal FP just as easily as any other C program.
We need first to examine the headers required for all software that intends to use libdfp. It is important to understand that libdfp headers are designed to extend the use of the existing system headers. libdfp headers will be picked up first during a program build and these will in turn include the standard system headers. In order for this extension scheme to work properly it is necessary to specify the following compile flags during any software build:
-std=gnu99 -I/path/to/dfp/include/ -D__STDC_WANT_DEC_FP__=1The variable can also be included within the C source file as follows:
#define __STDC_WANT_DEC_FP__
As an example, if libdfp was installed with "--prefix=/usr/local" then the include path will be "/usr/local/include/dfp."
Also note that "-std=gnu99" is an essential compiler flag for a correct build of any software that uses libdfp.
By using these compiler flags, when we include a file in our program source such as:
#include <math.h>
the file dfp/math.h will be picked up first which will then include the system "math.h" header. Ditto for other header files that are relevant to libdfp.
A simple program is now shown to calculate the amortisation of a loan using decimal FP values. This simple program will illustrate all the concepts needed to create more complex decimal FP programs with libdfp.
First the program code is given and then is followed with a line-by-line description. The C source file is also available here.
1 /*********************************************
2 gcc -std=gnu99 $CFLAGS -Wall -Wextra -I/usr/local/include/dfp -o dfp-amort dfp-amort.c -ldfp
3
4 ./dfp-amort
5
6 **********************************************/
7
8 #define _GNU_SOURCE
9 #define __STDC_WANT_DEC_FP__
10 #include <fenv.h>
11 #include <math.h>
12 #include <stdlib.h>
13 #include <stdio.h>
14 extern int __bid_fixddsi(_Decimal64 a);
15
16 int main() {
17
18 _Decimal64 P, pmt, interest, ann_rate, num_periods, per_rate, Ppmt;
19 int i;
20
21 // set rounding mode to round half up
22 fe_dec_setround(FE_DEC_TONEARESTFROMZERO);
23
24 ann_rate = 0.0329DD;
25 per_rate = ann_rate / 12.0DD;
26
27 P = 60000.0DD; // P = strtod64("60000.0", NULL);
28 num_periods = 60.0DD;
29
30 // The formula used to determine payment of an annuity (loan):
31 // pmt = P * per_rate/(1-(1+per_rate)^-n)
32 //
33 pmt = P * per_rate/(1.0DD - powd64((1.0DD + per_rate), -num_periods));
34 pmt = quantized64(pmt, 1.00DD);
35
36 fprintf(stdout, "Payment = %10.2DF\n\n", pmt);
37
38 fprintf(stdout, "%10s%12s%12s%16s%12s\n", "Period" , "Principal", "Interest", "Pmt on Princ", "Check");
39
40 i = 1;
41 interest = per_rate * P;
42 interest = quantized64(interest, 1.00DD);
43 Ppmt = pmt - interest;
44 fprintf(stdout, "%10d%12.2Df%12.2Df%16.2Df%12.2Df\n", i, P, interest, Ppmt, (interest + Ppmt));
45
46 for(i=2; i<=__bid_fixddsi(num_periods); i++) {
47 P = P - Ppmt;
48 interest = per_rate * P;
49 interest = quantized64(interest, 1.00DD);
50 Ppmt = pmt - interest;
51 fprintf(stdout, "%10d%12.2Df%12.2Df%16.2Df%12.2Df\n", i, P, interest, Ppmt, (interest + Ppmt));
52 }
53
54 }
This program calculates a loan amortisation schedule for a given principle value and a desired number of payment periods. The annual interest rate is defined as a constant within the program.
Programming in decimal floating point is essentially the same as with binary floating point except for one major difference. Decimal floating point (DECFP) numbers are not normalized. Hence, DECFP numbers can be equal but differ in their "quantum exponent," which is probably best understood as being the significant figures. For example, we can specify two constants that are equal in magnitude but contain different quantum exponents:
_Decimal64 X = 1000.00DD; // quantum exponent = -2 _Decimal64 Y = 1000.000DD; // quantum exponent = -3
Although X = Y they will behave differently in printing and also will produce different quantum exponents in the results of math operations.
Conceptually, each DECFP number is stored as a triplet: (sign, coefficient, quantum exponent). In the above example, X = (1, 100000, -2) and Y = (1, 1000000, -3).
The idea is simple although somewhat involved and cannot be covered thoroughly in this short article. But because of its importance to DECFP programming, the reader is referred to the standard references given in the libfp document:
What follows is a line-by-line description of the loan amortization code which illustrates some of the principles of DECFP programming.
Lines 1 - 6 show the GCC command and options that are used to compile this code. Here, I have named the program "dfp-amort" and line 4 shows how to invoke it. I have chosen not to have any input options to keep things as simple as possible, but command-line options would be identical to ordinary C programming.
Lines 8 - 13 state what header files are to be included. Also, the necessary variable "__STDC_WANT_DEC_FP__" is defined. As mentioned earlier, these header files are actually libdfp headers that will in turn include the corresponding system headers. This is accomplished by specifying "-D__STDC_WANT_DEC_FP__=1" in the GCC command or, as shown here, by defining it in the source file.
Line 14 is necessary because we are using one of the GCC built-in arithmetic functions. GCC functions are not defined in any header file but are linked in the accessory library libgcc. Hence we must declare these as external functions to avoid any "implicit function definition" warnings or errors.
As already indicated above, a listing of these libgcc functions is found in the GCC manual:
Lines 18 - 19 declare the variables used. Most are of type _Decimal64 except for the standard C int type.
Lines 21-22 set the DECFP rounding mode. The default mode is round-to-nearest with ties even, i.e. the same as binary FP. To make these monetary calculation behave as is commonly expected, the rounding is changed to round-to-nearest with ties upward. (Consult the "fenv.h" in the libdfp source for the defined constants.)
Lines 24 - 28 initialize the variables. Note that _Decimal64 constants are specified with the "DD" suffix. Also note I have indicated in the comment to line 24 the use of the strtod64() function from libdfp to accomplish the same initialization. This function format follows the standard C strtod() function.
As indicated above, all libdfp functions follow the naming and formatting of their standard C counterparts. This differs from other decimal FP libraries, such as the Intel library, and makes them very easy to apply.
Lines 30 - 32 shows the math formula used to compute the period payments based on the principal, period interest, and the number of periods desired. Line 33 encodes this formula using the powd64() function from libdfp.
Line 34 invokes the quantized64() function from libdfp to reduce the quantum exponent of the result to match that of US currency. It also, if necessary, takes care of rounding. This quantized64() function is one of the new math functions introduced by libdfp.
Line 36 outputs the period payment amount, a _Decimal64 value, using the standard fprintf() function from glibc. This is possible because libdfp employs "glib-hooks" to extend the fprintf() capabilities. To output a _Decimal64 value the "D" length modifier is used within the format string. A basic example of this is shown next:
fprintf(stdout, "A _Decimal64 value: %Df\n", some_dec64_variable);
Line 38 outputs a table header string and lines 40 - 52 computes the amortisation values, i.e. the period interest payment, the payment on the principle, and the remaining principle. A "check" column is included to verify that the sum of the interest and the payment on principle is the same for each period.
Note that in line 46 the upper limit of the for loop variable "i" is given by converting a _Decimal64 value to a signed integer with the __bid_fixddsi() function from libgcc.
Although only the _Decimal64 type is used here, the other _Decimal32 and _Decimal128 types are used in the same manner with the appropriate changes to the constant suffix, function name, and fprintf() length modifier.
Now follows the output of this C program using a principal value of $60,000.00 and with the number of periods being 60:
./dfp-amort
Payment = 1085.87
Period Principal Interest Pmt on Princ Check
1 60000.00 164.50 921.37 1085.87
2 59078.63 161.97 923.90 1085.87
3 58154.73 159.44 926.43 1085.87
4 57228.30 156.90 928.97 1085.87
5 56299.33 154.35 931.52 1085.87
6 55367.81 151.80 934.07 1085.87
7 54433.74 149.24 936.63 1085.87
8 53497.11 146.67 939.20 1085.87
9 52557.91 144.10 941.77 1085.87
10 51616.14 141.51 944.36 1085.87
11 50671.78 138.93 946.94 1085.87
12 49724.84 136.33 949.54 1085.87
13 48775.30 133.73 952.14 1085.87
14 47823.16 131.12 954.75 1085.87
15 46868.41 128.50 957.37 1085.87
16 45911.04 125.87 960.00 1085.87
17 44951.04 123.24 962.63 1085.87
18 43988.41 120.60 965.27 1085.87
19 43023.14 117.96 967.91 1085.87
20 42055.23 115.30 970.57 1085.87
21 41084.66 112.64 973.23 1085.87
22 40111.43 109.97 975.90 1085.87
23 39135.53 107.30 978.57 1085.87
24 38156.96 104.61 981.26 1085.87
25 37175.70 101.92 983.95 1085.87
26 36191.75 99.23 986.64 1085.87
27 35205.11 96.52 989.35 1085.87
28 34215.76 93.81 992.06 1085.87
29 33223.70 91.09 994.78 1085.87
30 32228.92 88.36 997.51 1085.87
31 31231.41 85.63 1000.24 1085.87
32 30231.17 82.88 1002.99 1085.87
33 29228.18 80.13 1005.74 1085.87
34 28222.44 77.38 1008.49 1085.87
35 27213.95 74.61 1011.26 1085.87
36 26202.69 71.84 1014.03 1085.87
37 25188.66 69.06 1016.81 1085.87
38 24171.85 66.27 1019.60 1085.87
39 23152.25 63.48 1022.39 1085.87
40 22129.86 60.67 1025.20 1085.87
41 21104.66 57.86 1028.01 1085.87
42 20076.65 55.04 1030.83 1085.87
43 19045.82 52.22 1033.65 1085.87
44 18012.17 49.38 1036.49 1085.87
45 16975.68 46.54 1039.33 1085.87
46 15936.35 43.69 1042.18 1085.87
47 14894.17 40.83 1045.04 1085.87
48 13849.13 37.97 1047.90 1085.87
49 12801.23 35.10 1050.77 1085.87
50 11750.46 32.22 1053.65 1085.87
51 10696.81 29.33 1056.54 1085.87
52 9640.27 26.43 1059.44 1085.87
53 8580.83 23.53 1062.34 1085.87
54 7518.49 20.61 1065.26 1085.87
55 6453.23 17.69 1068.18 1085.87
56 5385.05 14.76 1071.11 1085.87
57 4313.94 11.83 1074.04 1085.87
58 3239.90 8.88 1076.99 1085.87
59 2162.91 5.93 1079.94 1085.87
60 1082.97 2.97 1082.90 1085.87
Libdfp, in conjunction with the built-in GCC routines, provides complete decimal floating point capabilities to the programmer.
Hopefully, this discussion has adequately illustrated how to employ decimal floating point with C programming on GNU/Linux systems.