You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
17 KiB
17 KiB
Quad Double computation package
Copyright (C) 2003-2018
================================================
Revision date: 30 Oct 2018
Authors:
Yozo Hida U.C. Berkeley yozo@cs.berkeley.edu
Xiaoye S. Li Lawrence Berkeley Natl Lab xiaoye@nersc.gov
David H. Bailey Lawrence Berkeley Natl Lab dhbailey@lbl.gov
C++ usage guide:
Alex Kaiser Lawrence Berkeley Natl Lab adkaiser@lbl.gov
This work was supported by the Director, Office of Science, Division of Mathematical,
Information, and Computational Sciences of the U.S. Department of Energy under contract
number DE-AC02-05CH11231.
This work was supported by the Director, Office of Science, Division of Mathematical,
Information, and Computational Sciences of the U.S. Department of Energy under contract
numbers DE-AC03-76SF00098 and DE-AC02-05CH11231.
*** IMPORTANT NOTES:
See the file COPYING for modified BSD license information.
See the file INSTALL for installation instructions.
See the file NEWS for recent revisions.
See the file README.pdf for additional information. The file is mostly identical, but
but includes additional tables on constructors and constants.
Outline:
I. Introduction
II. Directories and Files
III. C++ Usage
IV. Fortran Usage
V. Note on x86-Based Processors (MOST systems in use today)
I. Introduction
This package provides numeric types of twice the precision of IEEE double (106 mantissa
bits, or approximately 32 decimal digits) and four times the precision of IEEE double (212
mantissa bits, or approximately 64 decimal digits). Due to features such as operator and
function overloading, these facilities can be utilized with only minor modifications to
conventional C++ and Fortran-90 programs.
In addition to the basic arithmetic operations (add, subtract, multiply, divide, square root),
common transcendental functions such as the exponential, logarithm, trigonometric and
hyperbolic functions are also included. A detailed description of the algorithms used is
available in the docs subdirectory (see docs/qd.ps). An abridged version of this paper,
which was presented at the ARITH-15 conference, is also available in this same directory
(see docs/arith15.ps).
II. Directories and Files
There are six directories and several files in the main directory of this distribution,
described below
src This contains the source code of the quad-double and double-double
library. This source code does not include inline functions,
which are found in the header files in the include directory.
include This directory contains the header files.
fortran This directory contains Fortran-90 files.
tests This directory contains some simple (not comprehensive) tests.
docs This directory contains two papers describing the algorithms.
config This directory contains various scripts used by the configure
script and the Makefile.
C++ Usage:
Please note that all commands refer to a Unix-type environment such as Mac OSX or Ubuntu
Linux using the bash shell.
A. Building
To build the library, first run the included configure script by typing
./configure
This script automatically generates makefiles for building the library and selects compilers
and necessary flags and libraries to include. If the user wishes to specify compilers or flags
they may use the following options.
CXX C++ compiler to use
CXXFLAGS C++ compiler flags to use
CC C compiler to use (for C demo program)
CFLAGS C compiler flags to use (for C demo program)
FC Fortran 90 compiler
FCFLAGS Fortran 90 compiler flags to use
FCLIBS Fortran 90 libraries needed to link with C++ code.
For example, if one is using GNU compilers, configure with:
./configure CXX=g++ FC=gfortran
The Fortran and C++ compilers must produce compatible binaries. On some systems
additional flags must be included to ensure that portions of the
library are not built with 32 and 64 bit object files. For example, on
64-Bit Mac OSX 10.6 (Snow Leopard) and 10.7 (Lion) the correct
configure line using GNU compilers is:
./configure CXX=g++ FC=gfortran FCFLAGS=-m64
To build the library, simply type
make
and the automatically generated makefiles will build the library including archive files.
To allow for easy linking to the library, the user may also wish to
install the archive files to a standard place. To do this type:
make install
This will also build the library if it has not already been built. Many systems, including Mac
and Ubuntu Linux systems, require administrator privileges to install the library at such
standard places. On such systems, one may type:
sudo make install
instead if one has sufficient access.
The directory "tests" contains programs for high precision quadrature and integer-relation
detection. To build such programs, type:
make demo
in the "tests" directory.
B. Linking
The simplest way to link to the library is to install it to a standard place as described above,
and use the -l option. For example
g++ compileExample.cpp -o compileExample -l qd
One can also use this method to build with make. A file called "compileExample.cpp" and the
associated makefile "makeCompileExample" illustrate the process.
A third alternative is to use a link script. If one types "make demo" in the test directory, the
output produced gives guidance as to how to build the files. By following the structure of
the compiling commands one may copy the appropriate portions, perhaps replacing the
filename with an argument that the user can include at link time. An example of such a
script is as follows:
g++ -DHAVE_CONFIG_H -I.. -I../include -I../include -O2 -MT $1.o -MD -MP -MF
.deps/qd_test.Tpo -c -o $1.o $1.cpp
mv -f .deps/$1.Tpo .deps/$1.Po
g++ -O2 -o $1 $1.o ../src/libqd.a Ðlm
To use it, make the link script executable and type:
./link.scr compileExample
Note that the file extension is not included because the script handles all extensions,
expecting the source file to have the extension ".cpp" .
Similarly, a script for compiling fortran programs may be constructed
as follows. In the fortran directory, type "make quadts". This
compiles the Fortran program tquadts.f, links with all necessary
library files, and produces the executable "quadts". As this is being
done, all flags and linked libraries are displayed. For instance, on
many Mac systems, presuming g++-4.0 was defined for C++ and gfortran
for F90, the following is output:
gfortran -O2 -ffree-form -c -o tquadts.o tquadts.f
g++-4.0 -O2 -Wall -o quadts tquadts.o second.o libarprec_f_main.a
libarprecmod.a ../src/libarprec.a
-L/usr/local/lib/gcc/i386-apple-darwin9.0.0/4.3.0
-L/usr/local/lib/gcc/i386-apple-darwin9.0.0/4.3.0/../../.. -lgfortranbegin
-lgfortran
Thus a general compile-link script (which could be saved in an
executable file named "complink.scr") is the following:
gfortran -O2 -ffree-form -c -o $1.o $1.f
g++-4.0 -O2 -Wall -o $1 $1.o second.o libarprec_f_main.a \
libarprecmod.a ../src/libarprec.a \
-L/usr/local/lib/gcc/i386-apple-darwin9.0.0/4.3.0 \
-L/usr/local/lib/gcc/i386-apple-darwin9.0.0/4.3.0/../../.. -lgfortranbegin \
-lgfortran
Note that if the .f90 suffix is used for Fortran-90 source files, the
-ffree-form flag can be omitted, and the first line above ends with
"$1.f90". Remember to type "chmod +x complink.scr". Then, for
instance, a program named "prog.f90" could be compiled and linked by
merely typing "./complink.scr prog".
C. Programming techniques
As much as possible, operator overloading is included to make basic programming as much
like using standard typed floating-point arithmetic. Changing many codes should be as
simple as changing type statements and a few other lines.
i. Constructors
To create dd_real and qd_real variables calculated to the proper precision, one must use
care to use the included constructors properly. Many computations in which variables are
not explicitly typed to multiple-precision may be evaluated with double-precision
arithmetic. The user must take care to ensure that this does not cause errors. In particular,
an expression such as 1.0/3.0 will be evaluated to double precision before assignment or
further arithmetic. Upon assignment to a multi-precision variable, the value will be zero
padded. This problem is serious and potentially difficult to debug. To avoid this, use the
included constructors to force arithmetic to be performed in the full precision requested.
For a table with descriptions, please see the included file README.pdf
ii. Included functions and Constants
Supported functions include assignment operators, comparisons, arithmetic and
assignment operators, and increments for integer types. Standard C math functions such as
exponentiation, trigonometric, logarithmic, hyperbolic, exponential and rounding functions
are included. As in assignment statements, one must be careful with implied typing of
constants when using these functions. Many codes need particular conversion for the power
function, which is frequently used with constants that must be explicitly typed for multi-
precision codes.
Many constants are included, which are global and calculated upon initialization. The
following list of constants is calculated for both the dd_real and qd_real classes separately.
Use care to select the correct value.
For a table with descriptions, please see the included file README.pdf
ii. Conversion of types
Static casts may be used to convert constants between types. One may also use constructors
to return temporary multi-precision types within expressions, but should be careful, as this
will waste memory if done repeatedly. For example:
qd_real y ;
y = sin( qd_real(4.0) / 3.0 ) ;
CÐstyle casts may be used, but are not recommended. Dynamic and reinterpret casts are
not supported and should be considered unreliable. Casting between multi-precision and
standard precision types can be dangerous, and care must be taken to ensure that programs
are working properly and accuracy has not degraded by use of a misplaced type-conversion.
D. Available precision, Control of Precision Levels,
The library provides greatly extended accuracy when compared to standard double
precision. The type dd_real provides for 106 mantissa bits, or about 32 decimal digits. The
type qd_real provides for 212 mantissa bits, or about 64 decimal digits.
Both the dd_real and qd_real values use the exponent from the highest double-precision
word for arithmetic, and as such do not extend the total range of values available. That
means that the maximum absolute value for either data type is the same as that of double-
precision, or approximately 10^308. The precision near this range, however, is greatly
increased.
To ensure that arithmetic is carried out with proper precision and accuracy, one must call
the function fpu_fix_start before performing any double-double or quad-double
arithmetic. This forces all arithmetic to be carried out in 64-bit double precision, not the 80-
bit precision that is found on certain compilers and interferes with the existing library.
unsigned int old_cw;
fpu_fix_start(&old_cw);
To return standard settings for arithmetic on oneÕs system, call the function Òfpu_fix_endÓ.
For example:
fpu_fix_end(&old_cw);
E. I/O
The standard I/O stream routines have been overloaded to be fully compatible with all
included data types. One may need to manually reset the precision of the stream to obtain
full output. For example, if 60 digits are desired, use:
cout.precision(60) ;
When reading values using cin, each input numerical value must start on a separate
line. Two formats are acceptable:
1. Write the full constant
3. Mantissa e exponent
Here are three valid examples:
1.1
3.14159 26535 89793
123.123123e50
When read using cin, these constants will be converted using full multi-precision accuracy.
IV. Fortran-90 Usage
NEW (2007-01-10): The Fortran translation modules now support the complex datatypes
"dd_complex" and "qd_complex".
Since the quad-double library is written in C++, it must be linked in with a C++ compiler (so
that C++ specific things such as static initializations are correctly handled). Thus the main
program must be written in C/C++ and call the Fortran 90 subroutine. The Fortran 90
subroutine should be called f_main.
Here is a sample Fortran-90 program, equivalent to the above C++ program:
subroutine f_main
use qdmodule
implicit none
type (qd_real) a, b
integer*4 old_cw
call f_fpu_fix_start(old_cw)
a = 1.d0
b = cos(a)**2 + sin(a)**2 - 1.d0
call qdwrite(6, b)
stop
end subroutine
This verifies that cos^2(1) + sin^2(1) = 1 to 64 digit accuracy.
Most operators and generic function references, including many mixed-mode type
combinations with double-precision (ie real*8), have been overloaded (extended) to work
with double-double and quad-double data. It is important, however, that users keep in
mind the fact that expressions are evaluated strictly according to conventional Fortran
operator precedence rules. Thus some subexpressions may be evaluated only to 15-digit
accuracy. For example, with the code
real*8 d1
type (dd_real) t1, t2
...
t1 = cos (t2) + d1/3.d0
the expression d1/3.d0 is computed to real*8 accuracy only (about 15 digits), since both d1
and 3.d0 have type real*8. This result is then converted to dd_real by zero extension before
being added to cos(t2). So, for example, if d1 held the value 1.d0, then the quotient d1/3.d0
would only be accurate to 15 digits. If a fully accurate double-double quotient is required,
this should be written:
real*8 d1
type (dd_real) t1, t2
...
t1 = cos (t2) + ddreal (d1) / 3.d0
which forces all operations to be performed with double-double arithmetic.
Along this line, a constant such as 1.1 appearing in an expression is evaluated only to real*4
accuracy, and a constant such as 1.1d0 is evaluated only to real*8 accuracy (this is
according to standard Fortran conventions). If full quad-double accuracy is required, for
instance, one should write
type (qd_real) t1
...
t1 = '1.1'
The quotes enclosing 1.1 specify to the compiler that the constant is to be converted to
binary using quad-double arithmetic, before assignment to t1. Quoted constants may only
appear in assignment statements such as this.
To link a Fortran-90 program with the C++ qd library, it is recommended to link with the
C++ compiler used to generate the library. The Fortran 90 interface (along with a C-style
main function calling f_main) is found in qdmod library. The qd-config script installed
during "make install" can be used to determine which flags to pass to compile and link your
programs:
"qd-config --fcflags" displays compiler flags needed to compile your Fortran files.
"qd-config --fclibs" displays linker flags needed by the C++ linker to link in all the
necessary libraries.
A sample Makefile that can be used as a template for compiling Fortran programs using
quad-double library is found in fortran/Makefile.sample.
F90 functions defined with dd_real arguments:
Arithmetic: + - * / **
Comparison tests: == < > <= >= /=
Others: abs, acos, aint, anint, asin, atan, atan2, cos, cosh, dble, erf,
erfc, exp, int, log, log10, max, min, mod, ddcsshf (cosh and sinh),
ddcssnf (cos and sin), ddranf (random number generator in (0,1)),
ddnrtf (n-th root), sign, sin, sinh, sqr, sqrt, tan, tanh
Similar functions are provided for qd_real arguments (with function
names qdcsshf, qdcssnf, qdranf and qdnrtf instead of the names in
the list above).
Input and output of double-double and quad-double data is done using the special
subroutines ddread, ddwrite, qdread and qdwrite. The first argument of these subroutines
is the Fortran I/O unit number, while additional arguments (as many as needed, up to 9
arguments) are scalar variables or array elements of the appropriate type. Example:
integer n
type (qd_real) qda, qdb, qdc(n)
...
call qdwrite (6, qda, qdb)
do j = 1, n
call qdwrite (6, qdc(j))
enddo
Each input values must be on a separate line, and may include D or E exponents. Double-
double and quad-double constants may also be specified in assignment statements by
enclosing them in quotes, as in
...
type (qd_real) pi
...
pi =
"3.14159265358979323846264338327950288419716939937510582097494459230"
...
Sample Fortran-90 programs illustrating some of these features are provided in the f90
subdirectory.
V. Note on x86-Based Processors (MOST systems in use today)
The algorithms in this library assume IEEE double precision floating point arithmetic. Since
Intel x86 processors have extended (80-bit) floating point registers, the round-to-double
flag must be enabled in the control word of the FPU for this library to function properly
under x86 processors. The following functions contains appropriate code to facilitate
manipulation of this flag. For non-x86 systems these functions do nothing (but still exist).
fpu_fix_start This turns on the round-to-double bit in the control word.
fpu_fix_end This restores the control flag.
These functions must be called by the main program, as follows:
int main() {
unsigned int old_cw;
fpu_fix_start(&old_cw);
... user code using quad-double library ...
fpu_fix_end(&old_cw);
}
A Fortran-90 example is the following:
subroutine f_main
use qdmodule
implicit none
integer*4 old_cw
call f_fpu_fix_start(old_cw)
... user code using quad-double library ...
call f_fpu_fix_end(old_cw)
end subroutine