A Fortran-90 Based Multiprecision System

yejet

README File for the Fortran-90 Multiprecision System

Contact:

David H. Bailey
NERSC, Lawrence Berkeley Lab
Mail Stop 50B-2239
Berkeley, CA 94720
Email: dhbailey@lbl.gov

Update as of 2005-03-11

--------------------------------------------

I. Copyright and Legal Matters
II. Description
III. Directories and Files
IV. Makefiles and System-Dependent Issues
V. Fortran Usage

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.

The software included in this package is a revision of software
earlier written while the author was an employee of NASA, in the
Numerical Aerospace Systems Division at the NASA Ames Research Center.

Copyright (c) 2004 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 contains a library of Fortran-90 multiprecision
computation routines.  It includes Fortran-90 translation modules,
which permit one to use the multiprecision library by making only
minor changes to existing Fortran-90 programs.  In most cases it is
only necessary to change the type declarations of certain variables,
plus two or three other very minor changes.  All usual arithmetic
operations are naturally extended, as are many transcendental
functions.  See "Fortran Usage", section V below for details.

A recent addition (Sept 2004) are module with array syntax
definitions, thus permitting one to use Fortran-90 array syntax
involving multiprecision variables.  These modules are contained in
file mpmod90v.f.  By default, they are NOT compiled or linked by the
Makefile.

These programs have recently been revamped to be in conformance with
the ARPREC package, which is available in the directory
http://crd.lbl.gov/~dhbailey/mpdist.  The ARPREC package is written in
C++ and provides the same high-level programming feature for C++
users, but also includes Fortran-90 translation modules similar to
those available in the MPFUN90 package.  Users are free to use
whichever package is best suited for their system -- a Fortran-90
multiprecision program that works with MPFUN90 should also work with
ARPREC, and vice versa.  Note however that, in the wake of the recent
updating of the MPFUN90 library, user programs with low-level calls,
which is not a recommended practice, may no longer work.

The "f90" directory in the MPFUN90 distribution package includes several
sample application programs, which may be run as is, or may be used as
templates for other programs.  These include programs for: (1) PSLQ
(integer relation detection), (2) quadrature (numerical integration),
(3) polynomial root-finding and (4) test programs.  The "toolkit"
directory in the MPFUN90 package contains code for the Experimental
Mathematician's Toolkit.  This program performs many multiprecision
computations via a simple interactive interface -- see the example
commands in the initial prompt.  The PSLQ, quadrature, root-finding,
and toolkit programs are exactly the same as the corresponding files
in the ARPREC package -- these programs are completely compatible, at
the Fortran-90 user level, with both packages.

III. Directories and Files

The following files are included in the mpfun90.tar.gz file:

Directory 'f90':

Makefile      Makefile for the f90 directory programs.
mpfun90.f     Main MPFUN90 library file.
testmp90.f    Test program for MPFUN-90 library file.
mpmod90.f     MPFUN90 translator modules (to enable high-level coding).
tmpmod90.f    Test program for translator modules.
tmpmod.out    Reference output of tmpmod90.f.
mpmodx90.f    Translator modules for using the advanced routines (>1000 digits).
tpslq1.f      One-level PSLQ (integer relation) program.
tpslq2.f      Two-level PSLQ (integer relation) program.
tpslq3.f      Three-level PSLQ (integer relation) program.
tpslqm1.f     One-level multi-pair PSLQ (integer relation) program.
tpslqm2.f     Two-level multi-pair PSLQ (integer relation) program.
tpslqm3.f     Three-level multi-pair pslq (integer relation) program.
tquadgs.f     Gaussian quadrature (numerical integration) program.
tquaderf.f    Error function quadrature (numerical integration) program.
tpslqts.f     Tanh-Sinh quadrature (numerical integration) program.
mpmodm90.f    Extra library file needed for tpslq3.f and tpslqm3.f.
mpmodx90.f    Extra library file needed for tquaderf.f and quadts.f.
mpmod90v.f    MPFUN90 array modules (to enable Fortran-90 array syntax).
second.f      Timing routine.
mpfun90.tex   Technical paper describing the mpfun90 package.
mpfun90.ps    PostScript of paper.

Directory 'toolkit':

Makefile      Makefile for the toolkit directory programs.
mpfun90.f     Main MPFUN90 library file.
mpmod90.f     MPFUN90 translator modules.
mathinit.f    Toolkit initialization program.  Must be compiled and
              executed (to produce some .dat files) before running mathtool.
mathtool.f    Main Toolkit program.
globdata.f    Defines some global data.
pslqsub.f     PSLQ subroutines.
quadsub.f     Quadrature subroutines.
rootsub.f     Root-finding subroutines.
zetapz.f      Multi-zeta function subroutines.
second.f      Timing routine.

IV. Makefiles and System-Dependent Issues

The two Makefiles are set up for an Apple or IBM system, using the IBM
xlf90 Fortran-90 comipler.  These Makefile will need to be changed for
other systems.  In many cases, it is simply necessary to change the
name of the compiler, and make a couple of simple changes to the
compiler flags.  Typical flags that may need to be set include a flag
to specify free-form Fortran source code format (instead of fixed
column format), and a flag to specify that a backslash (\) in a
literal string does not mean an escape character (this latter feature
is needed for the mathtool.f program).

On some systems, the 'second.f' file may need to be changed, typically
by deleting the underscore (_) after ETIME.  Other definitions of a
"second" (CPU time) function may be used.

These codes have been tested quite thoroughly, but a few bugs may
remain.  If you encounter any, please let me know and I will fix them
as soon as possible.

V. Fortran Usage

A. Introduction

The basic concept of the MPFUN90 package is to extend the Fortran-90
language, by means of completely standard Fortran-90 features
(operator overloading and custom datatypes), to perform computations
with an arbitrarily high level of numeric precision (hundreds or
thousands of digits).  In most cases, only minor modifications need to
be made to an existing Fortran program -- e.g., change the type
statements of variables that you wish to be treated as multiprecision
variables, plus a few other minor details.  Three multiprecision
datatypes are supported: mp_real, mp_int and mp_complex.

The MPFUN90 package is described in detail in the technical paper
mpfun90.tex, which is included in the MPFUN90 distribution package.
It can also be downloaded from
http://crd.lbl.gov/~dhbailey/dhbpapers/mpf90.pdf

Note, however, that there are a few differences between the package as
it is now and the original MPFUN90 package, which was documented in
the above paper.  These changes include: (1) several new
transcendental functions have been added; (2) the new MPFUN90 package
no longer supports mixed-mode operations between real*4 data and
multiprecision data, nor between complex*8 data and multiprecision
data (hardly anyone used these); (3) the procedure to set or modify
the working precision level has been changed; and (4) the procedure to
set or modify system parameters has been changed.  Details are given
below.

B. Fortran programming techniques

Several example application programs using the F90 arprec software can
be found in the f90 and the toolkit directory.

To use the package, first set the global variable mpipl, which is the
maximum precision level in digits, in the file mp_mod.f.  As
delivered, mpipl is set to 2000 digits, which means that the package
can handle user programs specifying up to 2000 digit accuracy.  If
mpipl is changed, mp_mod.f must be recompiled.  With mpipl set to 2000
digits, programs can define a working precision level to up and
including 2000 digits.

Modifying an existing Fortran-90 program to use the ARPEC library is
generally quite easy, because of the translation facilities in
mp_mod.f.  A sample user program is:

  program main
    use mpmodule
    implicit none
    type (mp_real) a, b
    call mpinit (500)
    a = 1.d0
    b = cos(a)**2 + sin(a)**2 - 1.d0
    call mpwrite(6, b)
    stop
  end program

This verifies that cos^2(1) + sin^2(1) = 1 to 500 digit accuracy.  The
line "use mpmodule", as shown above, must be included at the beginning
of each subroutine or function subprogram that uses multiprecision
datatypes.  Multiprecision variables are declared using a Fortran-90
defined type statement such as the following.

   type (mp_real) a, b, c(10)
   type (mp_integer) k1, k2, k3
   type (mp_complex) z1, z2(5,5), z3

Most operators and generic function references, including many
mixed-mode type combinations, have been overloaded (extended) to work
with multiprecision 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 real*4 or real*8 accuracy.
For example, with the code

   real*8 d1
   type (mp_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 fullly accurate
multiprecision quotient is required, this should be written:

  real*8 d1
  type (mp_real) t1, t2
   ...
  t1 = cos (t2) + mpreal (d1) / 3.d0

which forces all operations to be performed with multiprecision
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 multiprecision accuracy is required,
one should write

   type (mp_real) t1
   ...
   t1 = '1.1'

The quotes enclosing 1.1 specify to the compiler that the constant is
to be converted to binary using multiprecision arithmetic, before
assignment to t1.  Quoted constants may only appear in assignment
statements such as this.

One difference between this version of MPFUN90 and the original
version of MPFUN90 is that the package no longer supports mixed-mode
operations between "real" (single precision real, ie real*4) data and
multiprecision data, nor between "complex" (single precision real, ie
complex*8) data and multiprecision data.  If your programs have such
datatypes, convert these to real*8 and complex*16, respectively,
before attempting to change the program to use MPFUN90.

C. Functions defined with multiprecision arguments

F90 functions defined with mp_int arguments:
  Arithmetic:  + - * / **
  Comparison tests:  == < > <= >= /=
  Others: abs, max, min

f90 functions defined with mp_real arguments:
  Arithmetic:  + - * / **
  Comparison tests:  == < > <= >= /=
  Others: abs, acos, aint, anint, asin, atan, atan2, cos, dble, erf,
  erfc, exp, int, log, log10, max, min, mod, mpcsshr (cosh and sinh),
  mpcssnf (cos and sin), mpranf (random number generator in (0,1)),
  mpnrtf (n-th root), sign, sin, sqr, sqrt, tan

Note that several of the above transcendental functions are new --
they were not included in the original MPFUN90 package.

D.  Input/output of multiprecision data

Input and output of multiprecision data is performed using the special
subroutines mpread and mpwrite.  The first argument of these
subroutines is the Fortran I/O unit number, while additional arguments
(up to 9 arguments) are scalar variables or array elements of the
appropriate type.  Example:

   type (mp_real) a, b, c(n)
   ...
   call mpread (6, a, b)
   do j = 1, n
     call mpwrite (6, c(j))
   enddo

When using mpread, each input numerical value must start on a separate
line, and end with a comma.  Here are three valid examples:

   1.1,
   3.14159 26535 89793,
   10^3 x 3.14159 26535 89793 23846
     26433 83279 50288,

When read using mpread, these constants will be converted using full
multiprecision accuracy.

One can also read and write multiprecision variables and arrays using
Fortran unformatted (binary) I/O, as in

   type (mp_real) t1, a(30)
   write (11, t1)
   write (12, a)

Data written to a file in this fashion can be read with a similar
unformatted read statement, but only on the same system that it was
written on.  Unformatted files written by MPFUN90 programs are not
compatible with unformatted files written by ARPREC.

E. Handling precision level

The initial working precision, and the maximum for this run, is set by

   call mpinit (idig)

where idig is the desired precision in digits.  The current working
precision, in digits (idig) and words (iwds), can be obtained by

   call mpgetprec (idig)     [or]
   call mpgetprecwords (iwds)

(one word contains 48 bits, or about 14.44 digits).  The working
precision level may be changed by calling

   call mpsetprec (idig)     [or]
   call mpsetprecwords (iwds)

NOTE: the above procedure for setting, changing and fetching the
precision level is new.  This change was made to be compatible with
the ARPREC convention.

The output precision, which by default is limited to 56 digits no
matter what the working precision level, is given in the global
variable mpoud, and may be changed or fetched simply by writing

   mpoud = idig     [or]
   idig = mpoud

This may optionally be done by using subroutine calls:

   call mpsetoutputprec (idig)    [or]
   call mpgetoutputprec (idig)

F. Multiprecision system variables

A number of global variables are defined in the Fortran-90 wrapper,
which may be accessed and in some cases changed by the user.  Here is
a listing and brief description of these variables:

Integer parameters [cannot be changed by the user during execution]:

mpipl       Max prec level, digits              User sets in mp_mod.f
mpldb       Logical unit for some output (=6)   Set in mp_mod.f
mpnbt       Bits per word (= 48)                Set in mp_mod.f
mpwds       Max prec level in words             Set in mp_mod.f, from mpipl

Integer variables [may be changed by user during exeuction]:

mpoud       Number of digits output with mpwrite; default = 56

Integer variables [accessed via mpgetpar/mpsetpar (see note below)]:

mpidb      Debug level for arprec routines; default = 0
mpier      Error flag and error number; default = 0
mpird      Rounding mode (0, 1 or 2); default = 1
mpmcr      Threshold for advanced routines; default = 7
mpndb      Number of debug words output; default = 22
mpker      Array of 72 error handling flags; default each entry = 2

Multiprecision real parameters [calculated during initialization]:

mpl02                Log(2)
mpl10                Log(10)
mppic                Pi
mplrg                Very large mp_real value = 2^(48*2^27)
mpsml                Very small mp_real value = 2^(-48*2^27)

Multiprecision real variables [may be changed by user during exeuction]:

mpeps                Epsilon for user program; default = 10^(-prec)

The system parameters mpidb, mpier, mpird, mpmcr, mpndb and mpker can
be stored or fetched as follows:

   integer mpidb
   ...
   call mpsetpar ('mpndb', 10)         [or]
   call mpgetpar ('mpidb', mpidbx)

In MPFUN90, one can also actually directly set and reference these,
since they are global variables, e.g.

   mpndb = 10
   mpidbx = mpidb        [or]

However, this usage is *not recommended* since the resulting code is
not compatible with the ARPREC package.

yejet

http://citeseer.ist.psu.edu/bailey95fortran.html
http://crd.lbl.gov/~dhbailey/dhbpapers/mpf90.pdf