Quad Double computation package
Copyright (C) 2003-2012
================================================

Revision date:  13 Mar 2012

Authors:
Yozo Hida		U.C. Berkeley   			[email protected]
Xiaoye S. Li 		Lawrence Berkeley Natl Lab  	[email protected]
David H. Bailey  	Lawrence Berkeley Natl Lab  	[email protected]

C++ usage guide:
Alex Kaiser 		Lawrence Berkeley Natl Lab 	[email protected]

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