------------------------------------------------------------------------
| QUAD-DOUBLE/DOUBLE-DOUBLE COMPUTATION PACKAGE                        |
|                                                                      |
| Yozo Hida        U.C. Berkeley               yozo@cs.berkeley.edu    |
| Xiaoye S. Li     Lawrence Berekely Natl Lab  xiaoye@nersc.gov        |
| David H. Bailey  Lawrence Berkeley Natl Lab  dhbailey@lbl.gov        |

Revised  2004-12-16  Copyright (c) 2004
Revised  2002-10-22  Copyright (c) 2002
Revised  2001-05-31  Copyright (c) 2001

I. Copyright and Legal Matters
II. Description
III. Directories and Files
IV. Instructions and System-Dependent Issues
V. C++ Usage
VI. Fortran Usage
VII. Note on systems with Intel X86-based processors.

I. Copyright and Legal Matters

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-AC03-76SF00098.

Copyright (c) 2003, The Regents of the University of California,
through Lawrence Berkeley National Laboratory (subject to receipt of
any required approvals from U.S. Dept. of Energy) 

All rights reserved. 

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are 
met: 

(1) Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer. 
(2) Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution. 
(3) Neither the name of Lawrence Berkeley National Laboratory,
U.S. Dept. of Energy nor the names of its contributors may be used to
endorse or promote products derived from this software without
specific prior written permission. 

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
  

II. Description

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 advanced features such as
operator-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).

QD-DD Library

The QD-DD library uses object oriented features and operator-function
overloading, which is available in C++ and Fortran-90, to simplify
end-user programming.  Here is a sample C++ calling program:

    #include <iostream>
    #include "qd.h"
    using namespace std;

    int main() {
      qd_real a, b, c;

      cout << "Enter a:";
      cin >> a;
      cout << "Enter b:";
      cin >> b;
      c = a + b;
      cout << "a + b = " << c << endl;
      cout << "sin(c) = " << sin(c) << endl;

      return 0;
    }

Note that one could replace `qd_real' with `double', and the same
program will run using IEEE double instead of quad-double numbers.
One could equally well specify `dd_real' instead of `qd_real' for
double-double precision computation.  Furthermore, many mixed mode
operations are also supported.  For example, following the
declarations

    double d;
    dd_real dd;
    qd_real qd;

all of the expressions d + dd, dd * qd, qd / d, etc., are legal, and
call routines optimized for speed in each case.  Mixed-mode operations
are not allowed with integer types -- these must first be converted to
double floats.

Here is a sample Fortran-90 program, equivalent to the above C++
program:

  program main
    use qdmodule
    implicit none
    type (qd_real) a, b, c

    write (6, *) 'Enter a:'
    call qdread (5, a)
    write (6, *) 'Enter b:'
    call qdread (5, b)
    c = a + b
    write (6, *) 'a + b ='
    call qdwrite (6, c)
    stop
  end program


III. 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.

         
IV. Instructions and System-Dependent Issues

To build this library, follow the following steps.  In the future, the
autoconf utility will used to avoid some of these platform dependent
configuration issues.

  1. Edit the file "make.inc" to set compiler and platform specific options.
     For convenience, options for several common platforms (IBM SP R6000, 
     Sun USparc, Intel Linux, PPC Linux) are provided in the files
     "make.ibm", "make.sun", "make.x86", and "make.ppc".  Several options
     that affect the library build are specified in CC_OPTS and C_QD_OPTS 
     variables in the files.  The following CPP definitions are used in the
     library: 

         -DHAS_FMA
             If this option is specified in the CC_OPTS variable, then code
             suitable for systems with a fused multiply-accumulate instruction
             is generated.  This instruction exists in IBM Power3 and
             PowerPC processors.

         -Dx86  
             This option should be specified (in the CC_OPTS variable)
             for Intel x86 processors (and their clones).  This will disable 
             extended precision floating point registers, as they wreak havoc 
             with quad-double library routines.

         -DADD_UNDERSCORE
             If this is specified in the C_QD_OPTS variable, then a
             C wrapper function name will have an extra underscore (_)
             attached. This is required if the library is to be used by
             Fortran compilers that has this naming convention, such as
             Sun's f90 compiler.

         -DUPPER_CASE
             If this option is specified in the C_QD_OPTS variable, then a
             C wrapper function name will be capitalized.  This is required if
             the library is to be used by Fortran compilers with this naming
             convention, such as Cray's cf90 compiler.  Note that for the Cray
             systems, the flag -DCRAY_STRINGS is also required.

         -DCRAY_STRINGS
             This option (specified in C_QD_OPTS variable) is required when
             the library is to be used by cf90 Fortran compiler on the Cray 
             systems.  This allows for string passing, which has different 
             convention on the Cray systems.


  2. Edit the file make.config to set library options in the QD_OPTS variable:

     -DACCURATE
         If this flag is specified, then double-double and quad-double
         additions and subtractions will obey the IEEE floating-point model

               fl(a + b) = (1 + d) (a + b)

         where d is of the order eps.  Otherwise, addition and subtraction
         obeys a weaker, Cray-style error bound

               fl(a + b) = (1 + d1) a   +  (1 + d2) b

         where d1, d2 are of order eps.  The penalty is a substantial speed
         reduction.  Other operations are not affected.

     -DSLOPPY
         If this flag is specified, some routines become slightly less
         accurate but noticeably faster.  These mainly applies to
         multiplication and division.

     -DNO_INLINE
         This option disable any inline functions.  By default, many
         simple operations are inlined for benefit of speed; this option
         disables this.  Gives substantially faster compile time at the
         expense of execution speed.  Mostly used for debugging purposes.

  3. Type "make".  This will build the library "libqd.a".

  4. If Fortran 90 wrappers are needed, type "make wrapper", or else
     enter the fortran subdirectory and type "make".  This will
     compile the Fortran module files in the directory "fortran".  If
     problems are encountered, see the Fortran Usage section below.

  5. If you want the test programs to be compiled, type "make test".
     Note that depending on the compiler, pslq_test and or quadt_test
     may not compile. This is mainly due to lack of support for "Template"
     by the compilers.  However, the basic test qd_test should compile
     fine.  Run this test and make sure all test passes.


V. C++ Usage

In all C++ files that use double-double precision, include the header
file "dd.h".  In all C++ files that use quad-double precision, include
the header file "qd.h".  Note that some operators are defined as
inline functions in these header files, and the options -DHAS_FMA
-Dx86 -DSLOPPY -DNO_INLINE -DACCURATE affects these functions.  Thus
to compile any C++ file that uses double-double or quad-double
numbers, the appropriate flags must be passed to the compiler to
compile that file.

Writing C++ programs that utilize these routines is straightforward.
In many cases, it is only necessary to change the type statements for
variables that one wishes to define as quad-double or double-double.
Most operators, including the input-output operators << and >>, are
"overloaded", meaning that they are defined for the double-double and
quad-double datatypes, as well as combinations with type double.
Quad-double and double-double constants may be specified by means of
assignment statements, as in

  qd_real pi;
  ...
  pi = "3.14159265358979323846264338327950288419716939937510582097494459230";

Note that without the quotes, the constant on the right would be
truncated, and the assignment would not be fully precise.  A sample
test program that demonstrates many of these facilities is available
in the test suite.


VI. Fortran-90 Usage

A. Installation instructions

Unfortunately, there are differing conventions among computer vendors
as to how to handle a C function called from a Fortran program.  For
example, Sun's f90 compiler attaches an underscore to all function
names, while Cray's cf90 compiler makes the name all uppercase.  IBM's
compiler leaves these name alone.  To accomodate these different
conventions, two compiler flags are provided for compilation of the
library: -DADD_UNDERSCORE and -DUPPERCASE.  These will change the
names at compile time to satisfy these needs.

One needs to define the Fortran compiler in use by inserting a line
such as one of the following in the make.inc file (in the qd directory),
in the section "Compiler options":

F90 = xlf90 -O3 -qstrict [for IBM or Apple systems with the IBM xlf90 compiler]
F90 = pgf90 -Mfree  [for Intel systems with the Portland Group pgf90 compiler]

Fortran files needed for Fortran programs are included in the fortran
subdirectory, together with the application program fortran_test.f and
others.  Two module files, ddmod.f and qdmod.f, provide generic
functions and operator overloading from Fortran-90 programs.  These
can be compiled by typing "make" in the fortran subdirectory.

One can use either the F90 compiler or the C++ compiler to perform the
linking.  Either way, the two object files, ddmod.o and qdmod.o must
be linked with any Fortran program, together with the libqd.a file in
the qd directory, as well as F90 and C++ libraries.  If one uses the
F90 compiler to link, then the C++ libraries must be linked; if one
uses the C++ compiler to link, then the F90 libraries must be linked.

If one uses the F90 compiler to link, do the following to determine
which C++ libraries to use: first type the short C++ program

 #include "stdio.h"
 main() { printf("Hello world\n"); }

then compile it with the command "gcc -v test.cc" (or replace "gcc"
with the actual name of the C++ compiler being used).  By examining
the output, one can see what flags to use.  For instance, on an Apple
system with the IBM C++ compiler xlC and the IBM F90 compiler xlf90,
one obtains the list

 -L/opt/ibmcmp/xlsmp/1.4/lib  -L/opt/ibmcmp/vacpp/6.0/lib \
 -L/usr/lib/gcc/darwin/3.3  -L/usr/lib/gcc/darwin \
 -L/usr/libexec/gcc/darwin/ppc/3.3/../../.. \
 -lxlopt -lxl -libmc++ -lstdc++ -lgcc -lSystem -lstdc++

Thus the lines

LD = xlf90 -L/opt/ibmcmp/xlsmp/1.4/lib  -L/opt/ibmcmp/vacpp/6.0/lib \
 -L/usr/lib/gcc/darwin/3.3  -L/usr/lib/gcc/darwin \
 -L/usr/libexec/gcc/darwin/ppc/3.3/../../.. \
 -lxlopt -lxl -libmc++ -lstdc++ -lgcc -lSystem -lstdc++

must be included in the Makefile file (in the fortran subdirectory),
immediately after the line "QD_LIB = ../lib$(QD).a".

If one uses the C++ compiler to link, do the following to determine
which flags to use: first type the short F90 program

 program test
 write (6, *) 'hello world'
 stop
 end

the compile it with the command "f90 -v test.f" (or replace "f90" by
the actual name of the F90 compiler being used).  By examining the
output, one can see what flags to use.  For instance, on an Apple
system with xlC and xlf90 as before, one obtains the list

 -L/opt/ibmcmp/xlsmp/1.4/lib -L/opt/ibmcmp/xlf/8.1/lib \
 -L/usr/lib/gcc/darwin/3.3 -L/usr/lib/gcc/darwin \
 -L/usr/libexec/gcc/darwin/ppc/3.3/../../.. \
 -lxlf90 -lxlopt -lxlomp_ser -lxl -lxlfmath -lm -lc -lgcc -lSystem

Thus the lines

CPLUSPLUS_LIB = -L/opt/ibmcmp/xlsmp/1.4/lib -L/opt/ibmcmp/xlf/8.1/lib \
  -L/usr/lib/gcc/darwin/3.3 -L/usr/lib/gcc/darwin \
  -L/usr/libexec/gcc/darwin/ppc/3.3/../../.. \
  -lxlf90 -lxlopt -lxlomp_ser -lxl -lxlfmath -lm -lc -lgcc -lSystem

must be included in the Makefile file in the fortran subdirectory,
immediately after the line "QD_LIB = ../lib$(QD).a".

The sample Fortran programs, including fortran_test.f, may be compiled
by typing "make" while in the fortran subdirectory.  These sample
application files include tquaderq.f, tquaderq2d.f, tquadgsq.f
tquadtsq.f and tquadtsq2d.f, which perform 1-D or 2-D numerical
integration, using either the "error function" algorithm, the
"tanh-sinh" algorithm or Gaussian quadrature.

B. Fortran programming

The lines "use qdmodule" must be included at the beginning of each
routine that uses the quad-double type, and "use ddmodule" must be
included at the beginning of each routine that uses the double-double
type.  Double-double and quad-double variables are then specified by
means of Fortran-90 type statements, as in

   type (qd_real) a, b, c

A sample program is:

  program main
  use qdmodule 
  implicit none
  type (qd_real) a, b
  a = 1.d0
  b = cos(a)**2 + sin(a)**2 - 1.d0
  call qdwrite(6, b)
  stop
  end 

This verifies that cos^2(1) + sin^2(1) = 1 to 64 digit accuracy.  If
you are using a system with an Intel X86-based processor, three
additional lines are required -- see Section VII.

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 mp_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.

C. Functions defined with double-double and quad-double arguments

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).

D. Input and output of double-double and quad-double data

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.


VII. Note on Intel x86 Processors

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 flag -Dx86 will fill in the following
functions with appropriate code to facilitate manipulation of this
flag.  If the flag is not set, the functions do nothing (but still
exist).

x86_fix_start    This turns on the round-to-double bit in the control word.
x86_fix_end      This restores the control flag.

These functions must be called by the main program, as follows:

    int main() {
      unsigned short old_cw;
      x86_fix_start(&old_cw);

      ... user code using quad-double library ...

      x86_fix_end(&old_cw);
    }

A Fortran-90 example is the following:

    program foo
    use qdmodule
    implicit none
    integer*2 old_cw

    call x86_fix_start(old_cw)

      ... user code using quad-double library

    call x86_fix_end(old_cw)
    end program
