1MP USER'S GUIDE (FOURTH EDITION)
0THE CHARACTER SET USED FOLLOWS, WITH SPECIAL CHARACTERS (QUOTE, LESS
THAN, GREATER THAN, UNDERLINE) AND LOWER CASE LETTERS UNDERNEATH THE
CORRESPONDING UPPER CASE ONES.
0ABCDEFGHIJKLMNOPQRSTUVWXYZ$0123456789+-*/=(),. QUOTE LT GT UNDERLINE
abcdefghijklmnopqrstuvwxyz ' < > _
0This file is optionally distributed in upper-case only, using the
standard Fortran character set plus QUOTE ('), LESS THAN (<) and GREATER
THAN (>).
0The standard Fortran carriage-control characters are in column 1. Only
blank, '1', '+' and '0' are used. There are at most 57 lines per page,
and at most 72 characters (excluding the carriage-control character and
trailing blanks) per line.
1
0
0
0
0
0
0 MP USER'S GUIDE (Fourth Edition)
+ _________________________________
0
0
0 Richard P. Brent
0
0
0 Department of Computer Science
0 Australian National University
0 Box 4, Canberra, ACT 2600
0 Australia
0
0
0 Technical Report TR-CS-81-08
0 June 1981
0
0
0 Copyright 1981 by R. P. Brent
1 MP User's Guide Page 2
+ _______________
0CONTENTS
+________
0
1. General description of MP page 3
0 1.1 Introduction 3
1.2 Restrictions and assumptions 4
1.3 Note on correctness 5
1.4 History and likely future development 5
1.5 Summary of useful MP routines 7
0
2. EXAMPLE program 8
0
3. Parameters in COMMON /MPCOM/ 10
0 3.1 Description of the parameters 10
3.2 Methods of setting the parameters 12
3.3 Efficiency considerations 13
0 4. The Augment interface 14
0 4.1 Introduction 14
4.2 The description deck 15
4.3 JACOBI - an example of the use of Augment 18
0 5. Summary for Augment users 20
0 5.1 Notes 20
5.2 Summary 21
5.3 Synonyms 25
5.4 Reserved words 25
0
6. Description of MP routines 26
0
7. Test programs 68
0 7.1 TEST 68
7.2 TEST2 68
7.3 Other test programs 68
0 8. Comparison with earlier versions of MP 69
0 8.1 New routines and capabilities 69
8.2 Incompatibilities with earlier versions of MP 70
0 9. Installation instructions 71
0 9.1 Essentials 71
9.2 Conversion notes 72
9.3 Desirable changes 73
9.4 Note on Fortran 77 73
1Sec. 1.1 MP User's Guide Page 3
+ _______________
01. GENERAL DESCRIPTION OF MP
+____________________________
01.1 Introduction
+________________
0MP is a multiple-precision floating-point arithmetic package. It is
almost completely machine-independent, and should run on any machine
with an ANSI Standard Fortran (ANS X3.9-1966) compiler, sufficient
memory, and a wordlength (for integer arithmetic) of at least 16 bits.
A precompiler (Augment) which facilitates the use of the MP package is
available - see Section 4 for details.
0MP has been tested on various machines including a Univac 1108, a Univac
1100/82, a DEC 10, an IBM 360/50, 360/91 and 370/168, and a PDP 11/45,
with various Fortran compilers. These machines have effective integer
wordlengths ranging from 16 to 48 bits. A wordlength of 12 bits would
probably be sufficient, but this has not been tested.
0MP works with normalized floating-point numbers. The base (B) and
number of digits (T) are arbitrary, subject to some restrictions given
below, and may be varied dynamically.
0T-digit floating-point numbers are stored in integer arrays of dimension
at least T+2, with the following conventions
0 word 1 = sign (0, -1 or +1),
word 2 = exponent (to base B),
words 3 to T+2 = normalized fraction (one base B digit per word).
0Note that words 2 to T+2 are undefined if the sign is zero. Integer
arrays used for storage of T-digit floating-point numbers in the above
format will be referred to as (unpacked) multiple-precision variables or
(unpacked) multiple-precision numbers. There is also a packed
representation, where two base B digits are stored in each word - for
details see the description of MPPACK in Section 6.
0The user may ask for rounded or truncated arithmetic to be used. With
rounded arithmetic the basic arithmetic operations (addition,
multiplication and division) give the correctly rounded T-digit, base B
result. With truncated arithmetic the relative error in the result is
at most B**(1-T). Directed roundings (useful for interval arithmetic)
are also available - for details see Section 3.
0Rounding options are implemented for some elementary and special
functions (e.g. MPSQRT) but not yet for others. For details see Section
6. Most routines give at worst a result
y = f(x) which could be obtained by making an O(B**(1-T)) relative
perturbation in the argument x, evaluating the function f exactly, then
making an O(B**(1-T)) relative perturbation in the result.
0Exponents can lie in the range 1-M, ... , +M inclusive, where M is set
by the user. Thus, numbers from B**(-M) to almost B**M are
representable. On underflow during an arithmetic operation, the result
is set to zero by subroutine MPUNFL. (This is the default action and
may be modified - see Section 6.) On overflow subroutine MPOVFL is
called and execution is terminated with an error message.
1Sec. 1.2 MP User's Guide Page 4
+ _______________
0Error messages are printed on logical unit LUN, where LUN is set by the
user, and then execution is terminated by a call to subroutine MPERR.
0The parameters B, T, M, LUN etc. are stored in a labelled COMMON block,
COMMON /MPCOM/. A working array of sufficient size must be provided in
blank COMMON. Details are given in Section 3.1.
01.2 Restrictions and assumptions
+________________________________
0B (the base) must be at least 2.
0T (the number of digits) must be at least 2.
0M (the exponent range) must be greater than T and less than 1/4 of the
largest machine-representable integer.
08*B**2-1 must not exceed the largest machine-representable integer.
0B**(T-1) should be at least 10**7 (this restriction does not apply to
the basic arithmetic operations).
0B, T etc. must be initialized by calling either MPINIT or MPSET2 before
calling any other MP routines. For details see Section 6.
0Blank COMMON rather than labelled COMMON is used for working storage
because a curious restriction in the ANSI (1966) Fortran standard
requires that a labelled COMMON block be declared with the same length
in each subprogram in which it is declared, but this restriction
fortunately does not apply to blank COMMON. Although space in blank
COMMON must be provided for the MP routines to use as working storage,
this space need not start at the first word of blank COMMON (see the
description of MPSET2 in Section 6). Thus, programs which already use
blank COMMON can make use of the MP package.
0MP variables used as arguments of subroutines need not be distinct.
Thus
CALL MPADD (X, X, X)
and
CALL MPEXP (X, X)
0are allowable. However, distinct arrays which overlap should not be
used. For example, CALL MPEXP (X(1), X(2)) is not allowed if X is an
integer array holding a multiple-precision variable.
0
It is assumed that the compiler passes addresses of arrays used as
arguments in subroutine calls (i.e. call by reference), and does not
check for array bounds violations (either for arguments of subroutines
or for arrays in COMMON). Apart from these violations, MP is written in
a portable subset of ANSI Standard Fortran (ANSI X3.9-1966). This has
been checked by the PFORT verifier. The only machine-dependent routine
is MPUPW (which unpacks characters stored several to a word). Other
routines which may require trivial changes are MPINIT, MPIS and MPLARG
(see comments in Sections 6 and 9 below). For comments on Fortran 77,
see Section 9.4.
1Sec. 1.3 MP User's Guide Page 5
+ _______________
01.3 Note on correctness
+_______________________
0It is clearly impossible to thoroughly test even one MP routine with all
possible combinations of B, T, M etc. It is also impracticable to
formally prove correctness of any nontrivial MP routines using present
theorem-proving techniques. Thus, no warranty can be given that the MP
package is correct. To gain confidence in it we rely on
0 1. Careful modular design, using the idea of 'information hiding'
and the avoidance of 'tricks' and 'side-effects' as far as
possible.
0 2. Error analysis (in the style of Wilkinson).
0 3. Limited testing including, where possible, comparison of
results obtained by at least two independent methods (e.g.
MPLN, MPLNGS and MPLNI compute logarithms by three independent
methods).
0Although no absolute assurance is possible, a user may be confident that
a result obtained using MP is correct if
0 1. He uses more digits than are strictly necessary.
0 2. He repeats the computation using a different base B.
0 3. He uses two independent mathematical formulations of the
problem.
0It is unlikely that a bug in MP, the user's Fortran compiler, or the
user's hardware, would remain undetected if steps 1-3 were followed.
0MP makes minimal assumptions - primarily that INTEGER arithmetic (on
numbers which are not too large) is exact.
0MP does not use Fortran library routines such as SIN, ALOG, EXP, SQRT,
and uses only 'A1' format for output of multiple-precision numbers.
0REAL and DOUBLE PRECISION arithmetic are used only in routines for
conversion of these types to and from multiple-precision etc.
0
1.4 History and likely future development
+_________________________________________
0The first working version of MP (version 731101) was written in November
1973. Between 1973 and 1978 a considerable number of improvements and
additions were made, and nonstandard Fortran constructs were eliminated.
In April 1978 the Augment interface routines (see Section 4.1) were
added and the first version of the Augment description deck was written.
The version of MP published in TOMS and formerly available from the ACM
Algorithms Distribution service (version 770217) did not include Augment
interface routines, but a version of MP which includes them (version
780802) is now available from the Algorithms Distribution service. The
most recent version of MP may be obtained from the author at the address
given on the title page of this User's Guide.
1Sec. 1.4 MP User's Guide Page 6
+ _______________
0In 1979 the storage allocation scheme was improved, rounding options
were implemented for many MP routines, dependence on Fortran library
routines such as EXP and ALOG and on REAL arithmetic was eliminated,
many routines were modified to allow packed as well as unpacked
multiple-precision arguments, and corresponding changes were made in the
Augment description deck. For more details, see Section 8.1. The MP
User's Guide (Second Edition) was extensively revised to produce the
third edition, and at the same time it was converted from upper-case to
duo-case. Only minor changes and corrections have been made since then.
0For an introduction to the design and philosophy of MP, see - R. P.
Brent, 'A Fortran multiple-precision arithmetic package', ACM Trans.
Math. Software 4 (March 1978), 57-70. Additional details are given in
'Algorithm 524 - MP, a Fortran multiple-precision arithmetic package',
ibid, 71-81, and in Section 6 below. The mathematical theory behind the
algorithms used in the MP package is described in - R. P. Brent,
'Unrestricted algorithms for elementary and special functions',
Information Processing 80, North Holland, 1980, 613-619.
0As mentioned above, a preprocessor (Augment) which facilitates the use
of the MP package is available. See - R. P. Brent, J. A. Hooper and
J. M. Yohe, 'An Augment interface for Brent's multiple-precision
arithmetic package', ACM Trans. Math. Software 6 (1980), 146-149, and
Sections 4 and 5 below.
0Motivation for certain details of our implementation of T-digit, base B
arithmetic is given in - C. H. Reinsch, 'Principles and preferences for
computer arithmetic', ACM Signum Newsletter 14, 1 (March 1979), 12-27.
In designing the Augment interface we have been influenced by the IFIP
Working Group 2.5 proposals - see B. Ford (ed.), 'Parameterization of
the environment for transportable numerical software', ACM Trans. Math.
Software 4 (1978), 100-103. See also - W. S. Brown and S. I. Feldman,
'Environment parameters and basic functions for floating-point
computations', ACM Signum Newsletter 14, 1 (March 1979), 42-45.
0In the future we hope to implement rounding options for more MP
routines, and write a multiple-precision interval arithmetic package
which uses MP and takes advantage of the directed rounding options. An
Augment interface to facilitate the use of the interval arithmetic
package would be provided. (A similar project has been attempted by
John P. Jeter - see his 'A variable precision interval data type
extension to Fortran', M. Sc. thesis, Univ. of South West Louisiana,
July 1979.)
0A never-ending project is to implement multiple-precision versions of
ever more special functions, and to improve the efficiency of those
multiple-precision routines already implemented.
0Correspondence concerning MP should be addressed to Prof. R. P. Brent at
the address given on the title page. Reports of bugs, information on
applications for which the MP package has proved useful, and suggestions
for the future development of the package are particularly welcome.
1Sec. 1.5 MP User's Guide Page 7
+ _______________
01.5 Summary of useful MP routines
+_________________________________
0We mention here the names of those MP routines which are likely to be of
interest to someone using the MP package without the aid of the Augment
interface. Routines which are called by other MP routines but are
unlikely to be called directly by a user of the package are omitted.
For the user of Augment, the summary given in Section 5.2 will be more
useful. Details of all MP routines may be found in Section 6.
0Basic arithmetic
0 MPADD, MPADDI, MPADDQ, MPDIV, MPDIVI, MPMUL, MPMULI, MPMULQ,
MPREC, MPSUB
0Powers and roots
0 MPPWR, MPPWR2, MPQPWR, MPROOT, MPSQRT
0Elementary functions
0 MPASIN, MPATAN, MPATN2, MPCIS, MPCOS, MPCOSH, MPEXP, MPEXP1,
MPLG10, MPLN, MPLNI, MPLNS, MPSIN, MPSINH, MPTAN, MPTANH
0Special functions
0 MPBESJ, MPDAW, MPEI, MPERF, MPERFC, MPGAM, MPGAMQ, MPLI, MPLNGM
0Constants
0 MPBERN, MPEPS, MPEUL, MPMAXR, MPMINR, MPPI, MPZETA
0Input and output
0 MPFIN, MPFOUT, MPIN, MPOUT, MPUNFR, MPUNFW
0Conversion
0 MPCAM, MPCDM, MPCIM, MPCMD, MPCMDE, MPCMI, MPCMR, MPCMRE, MPCQM,
MPCRM
0Comparison
0 MPCMPA, MPCMPD, MPCMPI, MPCMPQ, MPCMPR, MPCOMP, MPEQ, MPGE, MPGT,
MPLE, MPLT, MPNE
0General utility routines
0 MPABS, MPCEIL, MPCHEB, MPCHEV, MPCMF, MPCMIM, MPDIGS, MPDIM,
MPFLOR, MPGCDA, MPGCDB, MPINIT, MPMAX, MPMIN, MPMOD, MPNEG,
MPPACK, MPPARA, MPPARB, MPPOLY, MPSETR, MPSET2, MPSIGN, MPSTR,
MPUNPK
0Example and test programs
0 EXAMPLE, JACOBI, TEST, TEST2
1Sec. 2 MP User's Guide Page 8
+ _______________
02. EXAMPLE PROGRAM
+__________________
0The following program illustrates the straightforward use of the MP
package without the aid of Augment. An example using Augment is given
in Section 4.3.
0
0C $$ ****** EXAMPLE ******
C
C THIS PROGRAM COMPUTES PI AND EXP(PI*SQRT(163/9)) TO 100
C DECIMAL PLACES, AND EXP(PI*SQRT(163)) TO 90 DECIMAL PLACES,
C AND WRITES THEM ON LOGICAL UNIT 6.
C
C CORRECT OUTPUT (EXCLUDING HEADINGS) IS AS FOLLOWS
C
C 3.14159265358979323846264338327950288419716939937510
C 58209749445923078164062862089986280348253421170680
C 640320.00000000060486373504901603947174181881853947577148
C 57603665918194652218258286942536340815822646477590
C 262537412640768743.99999999999925007259719818568887935385633733699086
C 2707537410378210647910118607312951181346
C
C WORKING SPACE IN BLANK COMMON, PARAMETERS IN
C COMMON /MPCOM/ (WHICH HAS LENGTH 23)
COMMON R
COMMON /MPCOM/ B, T, M, LUN, MXR, SPTR, MXSPTR, DUMMY
INTEGER B, T, M, LUN, MXR, SPTR, MXSPTR, DUMMY(16), R(500)
C
C WE HAVE T .LE. 62 IF WORDLENGTH AT LEAST 16 BITS, AND WORKING SPACE
C IS AT MOST 500 WORDS (LESS IF WORDLENGTH IS GREATER THAN 16 BITS).
C
C VARIABLES NEED T+2 .LE. 64 WORDS AND ALLOW 110 CHARACTERS FOR
C DECIMAL OUTPUT
INTEGER PI(64), X(64), C(110)
C
C CALL MPSET2 TO SET OUTPUT LOGICAL UNIT = 6 AND EQUIVALENT
C NUMBER OF DECIMAL PLACES TO AT LEAST 110. THE LAST THREE
C PARAMETERS ARE THE DIMENSIONS OF PI (OR X) AND THE LOWER
C AND UPPER INDICES OF BLANK COMMON AVAILABLE TO MP.
CALL MPSET2 (6, 110, 64, 1, 500)
C
C COMPUTE MULTIPLE-PRECISION PI
CALL MPPI(PI)
C
C CONVERT TO PRINTABLE FORMAT (F110.100) AND WRITE
CALL MPOUT (PI, C, 110, 100)
WRITE (LUN, 10) B, T, C
10 FORMAT (32H1EXAMPLE OF MP PACKAGE, BASE =, I9,
$ 12H, DIGITS =, I4 /// 11H PI TO 100D //
$ 11X, 60A1 / 21X, 50A1)
C
C SET X = SQRT(163/9), THEN MULTIPLY BY PI
CALL MPQPWR (163, 9, 1, 2, X)
CALL MPMUL (X, PI, X)
1Sec. 2 MP User's Guide Page 9
+ _______________
0C SET X = EXP(X)
CALL MPEXP (X, X)
C
C CONVERT TO PRINTABLE FORMAT AND WRITE
CALL MPOUT (X, C, 110, 100)
WRITE (LUN, 20) C
20 FORMAT (/ 28H EXP(PI*SQRT(163/9)) TO 100D //
$ 11X, 60A1 / 21X, 50A1)
C
C SET X = X**3 = EXP(PI*SQRT(163))
CALL MPPWR (X, 3, X)
C
C WRITE IN FORMAT F110.90
CALL MPOUT (X, C, 110, 90)
WRITE (LUN, 30) C
$ 1X, 70A1 / 21X, 40A1)
WRITE (LUN, 40) MXSPTR
40 FORMAT (/ 21H END OF EXAMPLE, USED, I4,
$ 23H WORDS OF WORKING SPACE //)
STOP
END
1Sec. 3.1 MP User's Guide Page 10
+ _______________
03. PARAMETERS IN COMMON /MPCOM/
+_______________________________
0
3.1 Description of the parameters
+_________________________________
0The labelled COMMON block, COMMON /MPCOM/, contains the following 23
parameters which are used by various routines in the MP package and
referred to in Section 6. In the description below, 'default value'
means the value set by a call to MPSET2 (or MPINIT or MPSET, since they
call MPSET2). All the parameters are of type INTEGER.
0
0(1) BASE or B The base of the multiple-precision number
representation. B must be at least 2, and 8*B**2-1
must be representable as a single-precision integer.
The default value is the largest power of two
satisfying this condition.
0(2) NUMDIG or T The number of multiple-precision digits (must be at
least 2). The default value is sufficient to give
accuracy equivalent to at least DECPL decimal places.
0(3) MAXEXP or M The maximum exponent of multiple-precision numbers.
Exponents in the range 1-M, ... , +M are allowed, so
numbers from B**(-M) to almost B**M are representable.
M must be greater than T. The default setting and
maximum allowable value is int (MXINT/4).
0(4) LUN The logical unit for output by routines MPFOUT, MPOUTF,
MPDUMP, MPERR, MPERRM etc. The default value is the
first argument of the last call to MPSET2 (or 6 if
MPINIT as distributed with the MP package is called).
0(5) MXR Words MNSPTR to MXR of blank COMMON are used as working
space by the MP package. MXR-MNSPTR must be
sufficiently large (see Section 6). The default values
of MNSPTR and MXR are the 4-th and 5-th arguments of
the last call to MPSET2 (or 1 and 500 if MPINIT is used
as distributed with the MP package).
0(6) SPTR The current pointer to the top of the working storage
stack (plus one). SPTR should be in the range
MNSPTR <= SPTR <= MXR + 1. The default initial setting
is MNSPTR.
0(7) MXSPTR The maximum value of SPTR during the current run. The
default initial setting is MNSPTR.
0(8) MNSPTR The minimum value of SPTR during the current run. See
the description of MXR above.
0(9) MXEXPN The maximum exponent of any multiple-precision number
during the current run. The default initial setting is
-M. Set to M+1 if multiple-precision overflow occurs.
1Sec. 3.1 MP User's Guide Page 11
+ _______________
0(10) MNEXPN The minimum exponent of any multiple-precision number
during the current run. The default initial setting is
M+1. Set to -M if a multiple-precision underflow
occurs.
0(11) RNDRL Integer in the range 0, ... , 3 indicating the type of
rounding to be performed. Generally,
0 0 means truncated (chopped) arithmetic,
1 means rounded arithmetic,
2 means round down (towards -infinity),
3 means round up (towards +infinity).
0 However, these interpretations do not always apply.
For details of the effect of RNDRL, see Section 6. The
default setting is 0 (truncated arithmetic).
0(12) KTUNFL Counter for the number of multiple-precision underflows
which have occurred. The default initial setting is
zero.
0(13) MXUNFL If positive, execution is terminated after MXUNFL
multiple-precision underflows have occurred. If zero,
any number of multiple-precision underflows are allowed
(this is the default).
0(14) DECPL A conservative estimate of the 'equivalent' number of
floating-point decimal places. The default setting is
the second argument in the last call to MPSET2 (or 43
if MPINIT is used as distributed). The default
settings of B and T satisfy
0 B**(T - 1) >= 10**(DECPL - 1)
0(15) MT2 The dimension of arrays used for multiple-precision
variables. The default setting is the third argument
in the last call to MPSET2 (or 27 if MPINIT is used as
distributed).
0(16) MXINT The largest single-precision integer of the form
2**k - 1 such that INTEGER arithmetic with operands and
results in the range -MXINT, ... , +MXINT is performed
exactly. Usually k = (wordlength in bits) - 1. See
also the description of MPLARG in Section 6.
0(17) EXWID The number of characters in the exponent field for
output of multiple-precision numbers using MPFOUT. The
default value is 6. (Since one position is taken by the
exponent character EXPCH and one by the sign, this
means that exponents in the range -9999 to +9999 are
allowable with the default setting.)
0(18) INRECL The input record length assumed by MPFIN. Must not
exceed 80 (the default value) unless MPFIN and MPCHK
are modified.
1Sec. 3.2 MP User's Guide Page 12
+ _______________
0(19) INBASE The base of 'decimal' numbers read by MPIN, MPFIN etc.
(e.g. INBASE = 8 means octal input, INBASE = 16 means
hexadecimal input). Must be in the range 2, ... , 16.
The default setting is ten, for decimal input.
0(20) OUTBAS The base for 'decimal' output by MPOUT, MPFOUT etc.
Default setting ten, must be in the range 2, ... , 16.
0(21) EXPCH The exponent character used by MPFOUT. Must not be a
digit, sign or blank. The default setting is 'E'.
0(22) CHWORD The number of characters per single-precision integer
word. The default setting is determined by MPUPW.
0(23) ONESCP 1 or 0 if the machine uses one's complement integer
arithmetic, 0 otherwise. The default setting is
determined by MPSET2 and MPUPW so that MPUPW works.
03.2 Methods of setting the parameters
+_____________________________________
0It would, of course, be possible to declare the parameters B, T, M, LUN,
... , ONESCP in COMMON /MPCOM/ and initialize them by assignment
statements or by a BLOCK DATA subprogram. However, this method is not
recommended. Instead, we recommend calling MPINIT or MPSET2 to set the
default values (see above and Section 6), and then using MPPARB or PARAM
to change the setting of any parameters for which the default values are
not desired.
0For example, to use the MP package with 20 place rounded floating
decimal arithmetic, we could do the following -
0 Without using Augment Using Augment
+ _____________________ _____________
0 CALL MPINIT (MP) INITIALIZE MP
CALL MPPARB (20, 'NUMDIG') PARAM ('NUMDIG') = 20
CALL MPPARB (10, 'BASE') PARAM ('BASE') = 10
CALL MPPARB (1, 'RNDRL') PARAM ('RNDRL') = 1
0MPPARC may be used instead of MPPARB. It is less mnemonic, but slightly
faster, so might be preferred inside a loop. For details see Section 6.
0The current setting of the parameters may be determined by MPPARA or
PARAM. For example, MPPARA ('NUMDIG') will return the current number of
multiple-precision digits. For users of Augment, PARAM ('NUMDIG') on
the right-hand side of an assignment statement will return the same
information.
0Augment users may prefer to set the first three parameters in COMMON
/MPCOM/ with the field functions BASE, NUMDIG and MAXEXP instead of
PARAM. For example, if X is a variable of type MULTIPLE, the number of
multiple-precision digits may be set to 20 by the statement
0 NUMDIG (X) = 20
0For further details, see Section 5.2.
1Sec. 3.3 MP User's Guide Page 13
+ _______________
03.3 Efficiency considerations
+_____________________________
0For efficiency choose B fairly large, subject to the restrictions given
in Section 1.2. For example, if the wordlength is
0 48 bits, could use B = 4194304 or 1000000,
36 bits, could use B = 65536 or 100000,
32 bits, could use B = 16384 or 10000,
24 bits, could use B = 1024 or 1000,
18 bits, could use B = 128 or 100,
16 bits, could use B = 64 or 10,
12 bits, could use B = 16 or 10.
0Of course, for special purposes such as simulation of decimal or binary
arithmetic, B may be set as small as 2. If a significant amount of
decimal input/output is to be performed, it may be best to choose B a
power of ten (MPIN and MPOUT will take advantage of this).
0Avoid multiplication and division by multiple-precision numbers, as
these take O(T**2) operations, whereas multiplication and division by
(single-precision) integers take O(T) operations. See the descriptions
of MPDIV, MPDIVI, MPMUL and MPMULI in Section 6.
0The default setting of RNDRL is 0 (i.e. truncated arithmetic). The use
of other values of RNDRL is more expensive in both time and space.
1Sec. 4.1 MP User's Guide Page 14
+ _______________
04. THE AUGMENT INTERFACE
+________________________
04.1 Introduction
+________________
0Augment is a preprocessor which allows the introduction of nonstandard
types (e.g. multiple-precision variables) into Fortran programs. For
details, see - F.D. Crary, 'A versatile precompiler for nonstandard
arithmetics', ACM Trans. Math. Software 5 (June 1979), 204-217, and the
references given there.
0A 'description deck' has been written to enable Augment to be used in
conjunction with the MP package. This greatly simplifies the task of
writing a program for a multiple-precision computation, or converting a
single (or double) precision routine to multiple precision.
0For example, if Augment is used we can declare X, Y and Z to be
multiple-precision variables by a declaration such as
0 MULTIPLE X, Y, Z
or
IMPLICIT MULTIPLE (A-H, O-Z)
0and subsequently write expressions such as
0 X = Y + Z*EXP(X+1)/Y
0instead of the equivalent
0 CALL MPADDI (X, 1, TEMP)
CALL MPEXP (TEMP, TEMP)
CALL MPMUL (Z, TEMP, TEMP)
CALL MPDIV (TEMP, Y, TEMP)
CALL MPADD (Y, TEMP, X)
0The Augment interface can be used with MP version 780420, or later
versions. For more details, see - R. P. Brent, J. A. Hooper and J. M.
Yohe, 'An Augment interface for Brent's multiple-precision arithmetic
package', ACM Trans. Math. Software 6 (1980), 146-149.
0A convenient summary of the operations implemented in the MP package and
the methods of invoking them via Augment is given in Section 5.2.
1Sec. 4.2 MP User's Guide Page 15
+ _______________
04.2 The description deck
+________________________
0The 'description deck' which describes the MP package is as follows.
0COMMENT AUGMENT DESCRIPTION DECK FOR THE MULTIPLE-PRECISION ARITHMETIC
PACKAGE OF R. P. BRENT.
0 THREE TYPES OF VARIABLE ARE DEFINED HERE -
0 MULTIPLE (STANDARD MULTIPLE-PRECISION NUMBERS),
MULTIPAK (PACKED MULTIPLE-PRECISION NUMBERS), AND
INITIALIZE (USED ONLY AS A DEVICE TO PERSUADE
AUGMENT TO INITIALIZE THE MP PACKAGE).
0 WORKING SPACE SHOULD BE ALLOCATED AND THE MP PACKAGE
INITIALIZED BY THE DECLARATION
0 INITIALIZE MP
0 IN THE MAIN PROGRAM.
0 THIS DESCRIPTION DECK ASSUMES THAT EACH MP NUMBER REQUIRES NO
MORE THAN 27 WORDS (14 IN PACKED FORMAT). THIS IS SUFFICIENT
FOR ABOUT 43 DECIMAL PLACE ACCURACY ON A 16-BIT MACHINE, AND
HIGHER ACCURACY ON A MACHINE WITH A LONGER WORDLENGTH.
0 SEE COMMENTS IN ROUTINE MPINIT FOR THE METHOD OF CHANGING THE
PRECISION, ALSO REGARDING COMMON DECLARATIONS.
0*DESCRIBE MULTIPAK
0DECLARE INTEGER (14), KIND SAFE SUBROUTINE, PREFIX MPK
0CONVERSION CTP (CIM, INTEGER, $, UPWARD), CTP (CAM, HOLLERITH)
0SERVICE COPY (=MPSTR)
0*DESCRIBE MULTIPLE
0DECLARE INTEGER (27), KIND SAFE SUBROUTINE, PREFIX MP
0OPERATOR + (, NULL UNARY, PRV, $), - (NEG, UNARY),
+ (ADD, BINARY3, PRV, $, $, $, COMM), * (MUL),
- (SUB,,,,,, NONCOMM), / (DIV), ** (PWR2),
0 + (ADDI,,,, INTEGER), * (MULI), / (DIVI), ** (PWR),
* (IMUL,,, INTEGER, $, $, NONCOMM),
0 .EQ. (EQ, BINARY2, PRV, $, LOGICAL, COMM), .NE. (NE),
.GE. (GE,,,,, NONCOMM), .GT. (GT), .LE. (LE), .LT. (LT),
0 + (, NULL UNARY, PRV, MULTIPAK),
- (NEG, UNARY, PRV, MULTIPAK),
+ (KADD, BINARY3, PRV, MULTIPAK, MULTIPAK, $, COMM),
* (KMUL), - (KSUB,,,,,, NONCOMM), / (KDIV), ** (PWR2),
1Sec. 4.2 MP User's Guide Page 16
+ _______________
0 * (KMLI,,,, INTEGER), / (KDVI), ** (PWR),
* (KIML,,, INTEGER, MULTIPAK, $, NONCOMM),
0 .EQ. (EQ, BINARY2, PRV, MULTIPAK, LOGICAL, COMM), .NE. (NE),
.GE. (GE,,,,, NONCOMM), .GT. (GT), .LE. (LE), .LT. (LT)
0TEST MPSIGA (SIGA, INTEGER)
0
FIELD SGN (SIGA, SIGB, ($), INTEGER), EXPON (EXPA, EXPB),
BASE (BASA, BASB), NUMDIG (DIGA, DIGB), MAXEXP (MEXA, MEXB),
DIGIT (DGA, DGB, ($, INTEGER)),
PARAM (PARA, PARB, (HOLLERITH)), FIO (FIN, FOUT, (INTEGER),
$), UNFIO (UNFR, UNFW),
EXPON (EXPA, EXPB, (MULTIPAK), INTEGER), NUMDIG (DIGA, DIGB),
SGN (SIGA, SIGB)
0
FUNCTION ABS (ABS, ($), $), DABS (ABS), ASIN (ASIN), DASIN (ASIN),
ARSIN (ASIN), DARSIN (ASIN), ARCSIN (ASIN), ATAN (ATAN),
DATAN (ATAN), ARTAN (ATAN), DARTAN (ATAN), ARCTAN (ATAN),
CMF (CMF), CEIL (CEIL), FLOOR (FLOR), CMIM (CMIM),
COS (COS), DCOS (COS), COSH (COSH), DCOSH (COSH),
DAW (DAW), DDAW (DAW), EI (EI), DEI (EI), ERF (ERF),
DERF (ERF), ERFC (ERFC), DERFC (ERFC), EXP (EXP),
DEXP (EXP), EXP1 (EXP1), FRAC (CMF), GAM (GAM),
DGAM (GAM), GAMMA (GAM), DGAMMA (GAM), AINT (CMIM),
DINT (CMIM), LI (LI), DLI (LI), LN (LN), LOG (LN),
ALOG (LN), DLOG (LN), LOG10 (LG10), ALOG10 (LG10),
DLOG10 (LG10), LNGM (LNGM), ALNGM (LNGM), DLNGM (LNGM),
LNGS (LNGS), LNS (LNS), REC (REC), SIN (SIN), DSIN (SIN),
SINH (SINH), DSINH (SINH), SQRT (SQRT), DSQRT (SQRT),
TAN (TAN), DTAN (TAN), TANH (TANH), DTANH (TANH),
0 SNGL (STR), DBLE (STR),
0 ART1 (ART1, (INTEGER)), LN (LNI), LNI (LNI), LOG (LNI),
ALOG (LNI), DLOG (LNI), ZETA (ZETA),
0 CAM (CAM), CAM (CAM, (HOLLERITH)),
0 MAX (MAX, ($, $)), AMAX1 (MAX), DMAX1 (MAX), MIN (MIN),
AMIN1 (MIN), DMIN1 (MIN), DIM (DIM), DDIM (DIM),
GCD (GCDA), ATAN2 (ATN2), DATAN2 (ATN2), ARTAN2 (ATN2),
MOD (MOD), AMOD (MOD), DMOD (MOD),
SIGN (SIGN), DSIGN (SIGN),
0 BESJ (BESJ, ($, INTEGER)), ROOT (ROOT),
0 MPINF (INF (SUBROUTINE), ($, INTEGER, INTEGER, HOLLERITH),
LOGICAL), MPOUTF (OUTF (SUBROUTINE)),
MPINF (INF (SUBROUTINE), ($, INTEGER, INTEGER, INTEGER)),
MPOUTF (OUTF (SUBROUTINE)),
0 INTEGR (INTG, ($), LOGICAL),
1Sec. 4.2 MP User's Guide Page 17
+ _______________
0 COMP (COMP, ($, $), INTEGER), CMPA (CMPA),
COMP (CMPI, ($, INTEGER)), COMP (CMPR, ($, REAL)),
COMP (CMPD, ($, DOUBLE PRECISION)),
COMP (CMPQ, ($, INTEGER, INTEGER)),
0 ADDQ (ADDQ, ($, INTEGER, INTEGER), $), MULQ (MULQ),
QPWR (QPWR, (INTEGER, INTEGER, INTEGER, INTEGER)),
CQM (CQM, (INTEGER, INTEGER)), CTM (CQM),
0 GAM (GAMQ), DGAM (GAMQ), GAMQ (GAMQ),
0 BERN (BERN, (INTEGER, INTEGER), MULTIPAK),
0 INT (CMI (SUBROUTINE), ($), INTEGER),
IDINT (CMI (SUBROUTINE)), IFIX (CMI (SUBROUTINE)),
0 ABS (ABS, (MULTIPAK), MULTIPAK), ATAN (ATAN,, $), COS (COS),
COSH (COSH), EXP (EXP), LN (LN), LOG (LN), ALOG (LN),
LOG10 (LG10), ALOG10 (LG10), SIN (SIN), SINH (SINH),
SQRT (SQRT), TAN (TAN), TANH (TANH),
EXP1 (EXP1), LNS (LNS),
0 MAX (MAX, (MULTIPAK, MULTIPAK)), MIN (MIN), DIM (DIM),
ATAN2 (ATN2), MOD (MOD), SIGN (SIGN),
ROOT (ROOT, (MULTIPAK, INTEGER))
0CONVERSION CTM (CDM, DOUBLE PRECISION, $, UPWARD), CTM (CIM, INTEGER),
CTM (CRM, REAL), CTM (UNPK, MULTIPAK), CTM (CAM, HOLLERITH),
CTD (CMD (SUBROUTINE), $, DOUBLE PRECISION, DOWNWARD),
CTI (CMI (SUBROUTINE),, INTEGER), CTR (CMR (SUBROUTINE),,
REAL), CTP (PACK,, MULTIPAK)
0SERVICE COPY (STR)
0*DESCRIBE INITIALIZE
0DECLARE INTEGER (1), KIND SAFE SUBROUTINE, PREFIX MPI
SERVICE COPY (STR), INITIAL (NIT)
0COMMENT END OF AUGMENT DESCRIPTION DECK FOR MP PACKAGE
1Sec. 4.3 MP User's Guide Page 18
+ _______________
04.3 JACOBI - an example of the use of Augment
+_____________________________________________
0The program which follows illustrates the use of the MP package with the
Augment interface. It is quite straightforward, and does not by any
means illustrate all the possibilities provided by the Augment
interface.
0
< Machine-dependent statement(s) to execute Augment >
PRINT SUPPRESS
< Machine-dependent statement(s) to insert the description deck >
0*BEGIN
*DISABLE PRINT, OUTPUT
*ENABLE SOURCE
0< Machine-dependent statement(s) to (later) execute Fortran compiler >
0C
C PROGRAM TO VERIFY AN IDENTITY OF JACOBI USING THE MP
C PACKAGE AND AUGMENT.
C
C THE PROGRAM READS A NUMBER X IN FREE-FIELD FORMAT ACCEPTABLE TO
C MPIN. IF X IS NON-POSITIVE IT HALTS. OTHERWISE IT COMPUTES
C AND PRINTS FN(X), FN(1/X) AND (FN(X)-FN(1/X))/FN(X),
C WHERE FN(X) IS THE SUM FROM N = -INFINITY TO +INFINITY OF
C SQRT(X)*EXP(-PI*(N*X)**2).
C THE IDENTITY VERIFIED IS
C FN(X) = FN(1/X)
C
C DECLARE VARIABLES AND INITIALIZE MP PACKAGE. ON SOME SYSTEMS BLANK
C COMMON MUST BE DECLARED HERE - SEE COMMENTS IN ROUTINE MPINIT.
C
MULTIPLE X, F1, F2, FN
INITIALIZE MP
C
C SET INPUT RECORD LENGTH TO 72 (TO AVOID READING
C SEQUENCE NUMBERS), DEFAULT VALUE IS 80.
C
PARAM (6HINRECL) = 72
C
C READ MP X FROM UNIT 5 IN FREE FORMAT, STOP IF X NOT POSITIVE.
C
10 X = FIO (5)
IF (X.LE.0) STOP
C
C WRITE HEADING, X, FN(X), AND FN(1/X) TO 40S IN STANDARD FORMAT
C
WRITE (6, 20)
20 FORMAT (//41H X, FN(X), FN(1/X), (FN(X)-FN(1/X))/FN(X)/)
FIO (40) = X
F1 = FN(X)
FIO (40) = F1
F2 = FN(1/X)
FIO (40) = F2
1Sec. 4.3 MP User's Guide Page 19
+ _______________
0C WRITE (F1-F2)/F1 TO 6 SIGNIFICANT FIGURES.
C NOTE THAT AN MP EXPRESSION CAN BE WRITTEN.
C
FIO (6) = (F1 - F2)/F1
GO TO 10
END
C
< Machine-dependent statement(s) to (later) execute Fortran compiler >
C
FUNCTION FN(X)
C
C RETURNS FN(X) = THE SUM FROM N = -INFINITY TO +INFINITY OF
C SQRT(X)*EXP(-PI*(N*X)**2), ASSUMING X POSITIVE.
C USES THE OBVIOUS METHOD, SO SLOW IF X SMALL.
C NOTE THAT X AND FN ARE BOTH TYPE MULTIPLE.
C TO ILLUSTRATE MIXED-MODE ARITHMETIC, SOME LOCAL
C VARIABLES ARE DECLARED TO BE PACKED.
C
MULTIPLE FN, X
MULTIPAK TM, FAC, PR
IF (X.LE.0) CALL MPERR
FN = 0
C
C AUGMENT CAN DEAL WITH THE FOLLOWING EXPRESSION AS IT KNOWS THAT X
C IS TYPE MULTIPLE, SO CALLS MPCAM TO CONVERT 'PI' TO MULTIPLE.
C
TM = EXP(-2HPI*X*X)
PR = TM
FAC = TM**2
C
C LOOP TO SUM INFINITE SERIES
C WARNING - NUMBER OF ITERATIONS IS PROPORTIONAL TO 1/X
C
10 FN = FN + TM
PR = PR*FAC
TM = TM*PR
C
C TEST FOR CONVERGENCE BY COMPARING EXPONENTS OF FN AND TM.
C WE COULD ALSO HAVE SAVED THE OLD VALUE OF FN AND SEEN IF
C STATEMENT 10 CHANGED IT.
C
IF (EXPON(FN)-EXPON(TM).LT.NUMDIG(X)) GO TO 10
FN = SQRT(X)*(2*FN+1)
RETURN
END
*END
0
.5
.3
1.+1
1.234567890123456789012345678901234567890123456789
0
1Sec. 5.1 MP User's Guide Page 20
+ _______________
05. SUMMARY FOR AUGMENT USERS
+____________________________
05.1 Notes
+_________
0We summarize as briefly as possible how MP routines may be invoked using
the Augment interface. Under 'Invocation' we give the typical way of
invoking the operation described under 'Description'. For further
details of calling sequences etc., the column 'Ref' gives a reference to
the relevant MP routine, and descriptions of these routines may be found
in Section 6. The column headed 'Res' gives the type of the result of
each operation.
0Types are indicated by one-letter abbreviations,
0 D = double-precision real (DOUBLE PRECISION),
G = packed or unpacked multiple (MULTIPAK or MULTIPLE)
(indicated by 'MP*' in the column headed 'Description'),
H = packed Hollerith (as in 'ABCDEF'),
I = single-precision integer (INTEGER),
L = Boolean (LOGICAL),
M = unpacked multiple (MULTIPLE),
P = packed multiple (MULTIPAK),
R = single-precision real (REAL).
0Arguments are written as XA, XB, XC, ... , where X indicates the type,
as above. XDUMMY means a dummy argument of the type indicated by X.
Such arguments serve no useful purpose, but Augment expects them.
0The conversion routines may be invoked implicitly by assignment
statements, for example -
0 MR = IA
0will be translated by Augment into a call to MPCIM.
0Packed (i.e. type MULTIPAK) and unpacked (i.e. type MULTIPLE)
multiple-precision variables may be mixed in simple arithmetic
expressions. For example, in the statement -
0 A = B + C/(D - E*F)
0any combination of packed and unpacked multiple-precision variables is
permissible. Augment will generate calls on the appropriate conversion
routines. We do not recommend using packed multiple-precision variables
in expressions involving function calls, without explicit type
conversion (using CTM and/or CTP), or unless it is mentioned as
permissible in the summary below. Augment will not always generate
calls on the appropriate conversion routines.
0Field functions may occur on either side of an assignment statement. For
example -
0 PARAM ('MXUNFL') = 10
0would set the maximum allowable number of MP underflows to 10.
1Sec. 5.2 MP User's Guide Page 21
+ _______________
05.2 Summary
+___________
0 Description Res Invocation Ref
+ ___________ ___ __________ ___
0 Addition
+ ________
0 Sum of two MP numbers M MA + MB MPADD
Sum of two packed MP numbers M PA + PB MPKADD
Sum of MP number and an integer M MA + IB MPADDI
Sum of MP number and rational IB/IC M ADDQ (MA, IB, IC) MPADDQ
0 Division
+ ________
0 Quotient of two MP numbers M MA / MB MPDIV
Quotient of two packed MP numbers M PA / PB MPKDIV
Quotient of MP number and an integer M MA / IB MPDIVI
Quotient of packed MP and integer M PA / IB MPKDVI
0 Multiplication
+ ______________
0 Product of two MP numbers M MA * MB MPMUL
Product of two packed MP numbers M PA * PB MPKMUL
Product of MP number and integer M MA * IB MPMULI
Product of packed MP and integer M PA * IB MPKMLI
Product of integer and MP number M IA * MB MPIMUL
Product of integer and packed MP M IA * PB MPKIML
Product of MP and rational IB/IC M MULQ (MA, IB, IC) MPMULQ
0 Reciprocal
+ __________
0 Reciprocal of MP number M REC (MA) MPREC
0 Subtraction
+ ___________
0 Difference of two MP numbers M MA - MB MPSUB
Difference of two packed MP numbers M PA - PB MPKSUB
0 Powers and roots
+ ________________
0 Raise MP* number to integer power M GA ** IB MPPWR
Raise MP number to MP power M MA ** MB MPPWR2
Raise packed MP to packed MP power M PA ** PB MPPWR2
Raise rat. IA/IB to rat. power IC/ID M QPWR(IA,IB,IC,ID) MPQPWR
IB-th root of MP* number M ROOT (GA, IB) MPROOT
0
Elementary functions
+ ____________________
0 Arcsin of MP number (in radians) M ASIN (MA) MPASIN
Arctan of MP* number (in radians) M ATAN (GA) MPATAN
Arctan of 1/IA, IA > 1 M ART1 (IA) MPART1
Arctan of MA/MB (both MP numbers) M ATAN2 (MA, MB) MPATN2
Arctan of PA/PB (both packed MP) M ATAN2 (PA, PB) MPATN2
Cosine of MP* number (radians) M COS (GA) MPCOS
Hyperbolic cosine of MP* number M COSH (GA) MPCOSH
1Sec. 5.2 MP User's Guide Page 22
+ _______________
0 Description Res Invocation Ref
+ ___________ ___ __________ ___
0 Elementary functions continued
+ ______________________________
0 Exponential of MP* number M EXP (GA) MPEXP
(Exponential - 1) of MP* number M EXP1 (GA) MPEXP1
Natural logarithm of MP* number M LOG (GA) MPLN
Logarithm (base 10) of MP* number M LOG10 (GA) MPLG10
Natural logarithm of (1 + MP* number) M LNS (GA) MPLNS
Natural logarithm of integer > 0 M LOG (IA) MPLNI
Sine of MP* number (radians) M SIN (GA) MPSIN
Hyperbolic sine of MP* number M SINH (GA) MPSINH
Square root of MP* number M SQRT (GA) MPSQRT
Tangent of MP* number (radians) M TAN (GA) MPTAN
Hyperbolic tangent of MP* number M TANH (GA) MPTANH
0
Special functions
+ _________________
0 Bessel function of first kind,
MP argument MA, integer order IB M BESJ (MA, IB) MPBESJ
Dawson's integral of MP argument M DAW (MA) MPDAW
Exponential integral of MP argument M EI (MA) MPEI
Error function of MP number M ERF (MA) MPERF
Complementary error function M ERFC (MA) MPERFC
Gamma function of MP number M GAM (MA) MPGAM
Gamma function of rational IA/IB M GAM (IA, IB) MPGAMQ
GCD of MP numbers which are integers M GCD (MA, MB) MPGCDA
Logarithmic integral of MP number M LI (MA) MPLI
Logarithm of Gamma function M LNGM (MA) MPLNGM
Riemann zeta function of integer > 1 M ZETA (IA) MPZETA
0
Constants
+ _________
0 Bernoulli numbers (result is array) P BERN (IA, IB) MPBERN
MP machine precision M CTM ('EPS') MPEPS
Euler's constant M CTM ('EUL') MPEUL
Largest positive MP number M CTM ('MAXR') MPMAXR
Smallest positive MP number M CTM ('MINR') MPMINR
Pi M CTM ('PI') MPPI
Decimal constants, e.g. 6.3E5 M CTM ('6.3E5$') MPCAM
0
Input and output
+ ________________
0 Read MP number from unit IA in
free format M FIO (IA) MPFIN
0 Write MP number MB with IA decimal
places on default output unit FIO (IA) = MB MPFOUT
0 Read IB words from unit IC under
format HD and convert to MP L MPINF(MA,IB,IC,HD)
number MA, result error flag MPINF
1Sec. 5.2 MP User's Guide Page 23
+ _______________
0 Description Res Invocation Ref
+ ___________ ___ __________ ___
0 Input and output continued
+ __________________________
0 Write IB words representing MP
number MA on default unit, with
IC places after decimal point, L MPOUTF(MA,IB,IC,HD)
under format HD, result error flag MPOUTF
0 Read MP number from unit IA,
unformatted M UNFIO (IA) MPUNFR
0 Write MP number MB on unit IA,
unformatted UNFIO (IA) = MB MPUNFW
0
0 Conversion
+ __________
0 Double-precision to multiple M CTM (DA) MPCDM
Integer to multiple M CTM (IA) MPCIM
Real to multiple M CTM (RA) MPCRM
Rational IA/IB to multiple M CTM (IA, IB) MPCQM
Packed multiple to unpacked multiple M CTM (PA) MPUNPK
Packed Hollerith to multiple M CTM (HA) MPCAM
Packed Hollerith stored in integer
or integer array to multiple M CAM (IA) MPCAM
Multiple to double-precision D CTD (MA) MPCMD
Multiple to integer (truncated) I CTI (MA) MPCMI
Multiple to real R CTR (MA) MPCMR
Unpacked multiple to packed multiple P CTP (MA) MPPACK
0
0 Comparison (result is +1, 0 or -1)
+ __________
0 Compare absolute value of MP numbers I CMPA (MA, MB) MPCMPA
Compare MP number with integer I COMP (MA, IB) MPCMPI
Compare MP number with real I COMP (MA, RB) MPCMPR
Compare MP number with double-prec. I COMP (MA, DB) MPCMPD
Compare MP number with rat. IB/IC I COMP (MA, IB, IC) MPCMPQ
Compare two MP numbers I COMP (MA, MB) MPCOMP
0
0 Relational
+ __________
0 Compare MP* numbers for equality L GA .EQ. GB MPEQ
(Similarly for .GE., L GA .GE. GB MPGE
.GT., .LE., .LT. and .NE.) etc. etc.
0
Test
+ ____
0 Three-way branch (n1, n2, n3 labels) IF (MA) n1,n2,n3 MPSIGA
1Sec. 5.2 MP User's Guide Page 24
+ _______________
0 Description Res Invocation Ref
+ ___________ ___ __________ ___
0 Field functions
+ _______________
0 Sign of MP* number I SGN (GA) MPSIGA
MPSIGB
Exponent of MP* number I EXPON (GA) MPEXPA
MPEXPB
IB-th digit of MP number I DIGIT (MA, IB) MPDGA
MPDGB
Number of MP digits (T) I NUMDIG (GDUMMY) MPDIGA
MPDIGB
Maximum exponent of MP numbers (M) I MAXEXP (MDUMMY) MPMEXA
MPMEXB
Multiple-precision base (B) I BASE (MDUMMY) MPBASA
MPBASB
Parameter in COMMON /MPCOM/, I PARAM (HA) MPPARA
argument HA may be any of 'BASE', MPPARB
'NUMDIG', 'MAXEXP', 'LUN', 'SPTR',
'MXSPTR', 'MNSPTR', 'MXEXPN',
'MNEXPN', 'RNDRL', 'KTUNFL',
'MXUNFL', 'DECPL', 'MT2', 'MXINT',
'EXWID', 'INRECL', 'INBASE',
'OUTBAS', 'EXPCH', 'CHWORD',
or 'ONESCP'
0 See under 'Input and output' M FIO (IA) MPFIN
MPFOUT
See under 'Input and output' M UNFIO (IA) MPUNFR
MPUNFW
0 Miscellaneous
+ _____________
0 Unary minus of MP number M -MA MPNEG
Unary minus of packed MP number P -PA MPNEG
Assignment of MP number (result = MA) M MA MPSTR
Assignment of packed MP number P PA MPSTR
Absolute value of MP number M ABS (MA) MPABS
Absolute value of packed MP number P ABS (PA) MPABS
Fractional part of MP number M FRAC (MA) MPCMF
Integer part of MP number M INT (MA) MPCMIM
Ceiling function of MP number M CEIL (MA) MPCEIL
Floor function of MP number M FLOOR (MA) MPFLOR
Returns true if MA is exact integer L INTEGR (MA) MPINTG
Maximum of two MP numbers M MAX (MA, MB) MPMAX
Maximum of two packed MP numbers M MAX (PA, PB) MPMAX
Minimum of two MP numbers M MIN (MA, MB) MPMIN
Minimum of two packed MP numbers M MIN (PA, PB) MPMIN
Max (0, MA - MB) for MP numbers M DIM (MA, MB) MPDIM
Max (0, PA - PB) for packed numbers M DIM (PA, PB) MPDIM
MA - int(MA/MB)*MB for MP numbers M MOD (MA, MB) MPMOD
PA - int(PA/PB)*PB for packed numbers M MOD (PA, PB) MPMOD
abs(MA)*sign(MB) for MP numbers M SIGN (MA, MB) MPSIGN
abs(PA)*sign(PB) for packed numbers M SIGN (PA, PB) MPSIGN
Initialization of MP package INITIALIZE MP MPINIT
1Sec. 5.3 MP User's Guide Page 25
+ _______________
05.3 Synonyms
+____________
0In some cases Augment will recognize synonyms for the names given in the
column headed 'Invocation' in the summary above. For example, AMAX1 may
be used in place of MAX. A list of synonyms follows. In most cases
these are usable only for unpacked multiple-precision arguments. It is
easy to modify or expand the list by making appropriate changes to the
Augment description deck.
0 Generic name Synonyms
+ ____________ ________
0 ABS DABS
AINT DINT
ASIN ARCSIN, ARSIN, DARSIN, DASIN
ATAN ARCTAN, ARTAN, DARTAN, DATAN
ATAN2 ARTAN2, DATAN2
COS DCOS
COSH DCOSH
DAW DDAW
DIM DDIM
EI DEI
ERF DERF
ERFC DERFC
EXP DEXP
GAM DGAM, DGAMMA, GAMMA,
GAMQ (rational argument only)
INT IDINT, IFIX
LI DLI
LOG ALOG, DLOG, LN, LNI (integer argument only)
LNGM ALNGM, DLNGM
LOG10 ALOG10, DLOG10
MAX AMAX1, DMAX1
MIN AMIN1, DMIN1
SIGN DSIGN
MOD AMOD, DMOD
SIN DSIN
SINH DSINH
SQRT DSQRT
TAN DTAN
TANH DTANH
05.4 Reserved words
+__________________
0When writing programs which use MP via the Augment interface, it is best
to avoid using the following identifiers except with their reserved
meanings as indicated above -
0 INITIALIZE,
MPxxxx (for any letters or digits xxxx),
MULTIPAK,
MULTIPLE.
0Care should also be taken to avoid using the names given in the column
headed 'Invocation' in the summary above, or in the list of synonyms
above, in a manner which would confuse Augment.
1Sec. 6 MP User's Guide Page 26
+ _______________
06. DESCRIPTION OF MP SUBROUTINES
+________________________________
0We give first the method of calling each multiple-precision routine
directly, second (third, ...) alternative methods (if any) using the
Augment interface described in Sections 4 and 5.
0Unless otherwise noted, X, Y and Z denote integer arrays representing
multiple-precision variables (often called (unpacked) multiple-precision
numbers, or (unpacked) multiple-precision variables, or variables of
type MULTIPLE), ERR and LV denote LOGICAL variables, I, J, K, L, IX etc.
denote single-precision INTEGER variables, RX, RY etc. denote
single-precision REAL variables, DX, DY etc. denote DOUBLE PRECISION
variables, and A denotes a packed Hollerith string.
0For definitions of B, T, M, LUN, MXR, RNDRL, INBASE, OUTBAS, ... see
Section 3.1. 'Space' means the size of the work area in COMMON, i.e.
MXR+1-MNSPTR (see Section 3.1). MPGAM, MPLNGM and MPZETA require
O(T**2) space, other MP routines require O(T) space.
0Time and space bounds such as O(T**2) are as T tends to infinity with
everything else fixed. For the definition of the function m(T)
appearing in some time bounds, see the description of MPMUL below.
0MPABS
+_____
0CALL MPABS (X, Y) or Y = ABS (X)
0Sets Y = abs(X) for multiple-precision variables X and Y. X and Y may
be both packed or both unpacked. This is an exception to the general
rule that an unpacked result is returned even if the argument(s) are
packed. The result is exact.
0MPADD
+_____
0CALL MPADD (X, Y, Z) or Z = X + Y
0Adds X and Y, forming result in Z, where X, Y and Z are
multiple-precision numbers.
0Rounding is defined by the parameter RNDRL in COMMON /MPCOM/ (see
Section 3.1) as follows -
0 RNDRL = 0 - Truncate towards zero if X*Y >= 0, away from zero if
X*Y < 0, in both cases using one guard digit, so the
result is exact if severe cancellation occurs.
0 RNDRL = 1, 2 or 3 - See description of MPCQM. Sufficient guard
digits are used to give the best possible result.
0MPADDI
+______
0CALL MPADDI (X, IY, Z) or Z = X + IY
0Adds multiple-precision X to integer IY giving multiple-precision Z.
Rounding options are implemented as for MPADD.
1Sec. 6 MP User's Guide Page 27
+ _______________
0MPADDQ
+______
0CALL MPADDQ (X, I, J, Y) or Y = ADDQ (X, I, J)
0Adds the rational number I/J to the multiple-precision number X, giving
a multiple-precision result in Y. For the effect of rounding options
see the descriptions of MPCQM and MPADD.
0MPADD3
+______
0Called by MPADD, does inner loops of addition. Not recommended for
independent use.
0MPART1
+______
0CALL MPART1 (N, Y) or Y = ART1 (N)
0Computes multiple-precision Y = arctan(1/N), for integer N > 1. Uses
the Taylor series arctan(x) = x - x**3/3 + x**5/5 - ... , where x = 1/N.
Rounding options are implemented as for MPEXP.
0
MPASIN
+______
0CALL MPASIN (X, Y) or Y = ASIN (X)
0Returns Y = arcsin(X), for multiple-precision numbers X and Y, where
abs(X) <= 1. The result is in the range -pi/2 to +pi/2. Uses MPATAN,
so time = O(m(T)T). Rounding options are not yet implemented, no guard
digits used.
0
MPATAN
+______
0CALL MPATAN (X, Y) or Y = ATAN (X)
0Returns Y = arctan(X) for (packed or unpacked) multiple-precision X,
unpacked multiple-precision Y, using an O(m(T)T) method. An
asymptotically faster method would be to combine Newton's method with
MPTAN. For another asymptotically faster method, see - R. P. Brent,
'Fast multiple-precision evaluation of elementary functions', J. ACM 23
(1976), 242-251, and the description of MPPIGL. Result is in the range
-pi/2 to +pi/2. Rounding options are implemented as for MPEXP.
0
MPATN2
+______
0CALL MPATN2 (X, Y, Z) or Z = ATAN2 (X, Y)
0Returns Z = arctan (X/Y) if Y nonzero,
Z = pi/2 if Y = 0,
0where X and Y are (both packed or both unpacked) multiple-precision
variables, Z is an unpacked multiple-precision variable. Rounding
options are implemented as for MPEXP.
1Sec. 6 MP User's Guide Page 28
+ _______________
0MPBASA
+______
0I = MPBASA (X) or I = BASE (X)
0Returns the multiple-precision base (first word of COMMON /MPCOM/ - see
Section 3.1). X is a dummy multiple-precision argument (Augment expects
one), I is an integer.
0MPBASB
+______
0CALL MPBASB (I, X) or BASE (X) = I
0Sets the multiple-precision base (first word of COMMON /MPCOM/) to I.
I should be an integer such that I > 1 and (8*I*I-1) is representable as
a single-precision integer (see Section 3.1). X is a dummy
multiple-precision argument (Augment expects one).
0Setting the base does not automatically convert multiple-precision
numbers to the new base. This can be done by converting to decimal
using MPFOUT, changing the base, and converting back using MPFIN.
0MPBERN
+______
0CALL MPBERN (N, P, X) or X(1) = BERN (N, (PARAM ('MT2') + 1)/2)
0Computes the Bernoulli numbers B(2) = 1/6, B(4) = -1/30, B(6) = 1/42,
B(8)= -1/30, ... , B(2N), which are defined by the generating function
y/(exp(y)-1). N and P are single-precision integers, with 2*P > T+1.
0For use without Augment, X should be a one-dimensional integer array of
dimension at least P*N. The Bernoulli numbers B(2), B(4), ... , B(2N)
are returned in packed format in X, with B(2J) in locations
X((J-1)*P+1), ... , X(J*P). Thus, to get B(2J) in unpacked
multiple-precision format in Y, one should CALL MPUNPK (X(IX), Y) after
calling MPBERN, where IX = (J-1)*P+1.
0Alternatively (simpler but nonstandard) - X may be a two-dimensional
integer array declared with dimension (P, N1), where N1 >= N and
2*P > T+1. Then B(2), B(4), ... , B(2N) are returned in packed format
in X, with B(2J) in X(1,J), ... , X(P,J). Thus, to get B(2J) in the
usual multiple-precision format in Y one should CALL MPUNPK (X(1, J), Y)
after calling MPBERN.
0For use with Augment, declare
0 MULTIPAK X(N1)
0where N1 >= N, and use (PARAM('MT2') + 1)/2 as the second argument.
0The well-known recurrence is unstable (losing about 2J bits of relative
accuracy in the computed B(2J)), so we use a different recurrence
derived by multiplying the generating function y/(exp(y)-1) by sinh(y/2)
and equating coefficients. (This method was suggested by Christian
Reinsch, and is faster than the method used in earlier versions of
MPBERN.)
1Sec. 6 MP User's Guide Page 29
+ _______________
0The relation
0 B(2J) = -2((-1)**J)factorial(2J)zeta(2J)/((2pi)**(2J))
0is used if zeta(2J) equals 1 to working accuracy. The relative error in
B(2J) is O((J**2)*(B**(1-T))). For further details, see - R. P. Brent,
'Unrestricted algorithms for elementary and special functions',
Information Processing 80, North Holland, 1980, 613-619.
0Time = O(T*min(N,T)**2 + N*m(T)).
0Rounding options are not implemented, no guard digits used.
0The above remarks assume that N is positive. If N is negative, they
apply with N replaced by abs(N), except that the precision decreases
linearly, from full precision for B(2) to low precision for B(2N). This
is faster, and sufficiently accurate for most applications where the
Bernoulli numbers are required to generate coefficients in an
Euler-Maclaurin sum formula.
0
0MPBESJ
+______
0CALL MPBESJ (X, NU, Y) or Y = BESJ (X, NU)
0Returns Y = J (NU, X), the first-kind Bessel function of order NU, for
small integer NU, multiple-precision X and Y. The method used is
Hankel's asymptotic expansion if abs(X) is large, the power series if
abs(X) is small, and the backward recurrence method otherwise. Results
for negative arguments are defined by
0 J (-NU, X) = J (NU, -X) = ((-1)**NU) * J (NU, X).
0Time = O(m(T)T) for fixed X and NU, increases as X and NU increase,
unless X is large enough for asymptotic series to be used. Space
required is large (compared to most multiple-precision routines), but
O(T).
0Rounding options are not yet implemented. Uses truncated arithmetic
with some guard digits.
0MPBES2
+______
0Uses the backward recurrence method to evaluate J (NU, X), where X is a
multiple-precision variable, NU (the index) is an integer, and J is the
Bessel function of the first kind. Assumes that NU >= 0 and X > 0. For
normalization the identity
0 J(0,X) + 2*J(2,X) + 2*J(4,X) + ... = 1
0is used. Called by MPBESJ and not recommended for independent use.
Rounding options are not implemented, no guard digits used.
1Sec. 6 MP User's Guide Page 30
+ _______________
0MPCAM
+_____
0CALL MPCAM (A, X) or X = CTM (A) or X = CAM (A)
0Converts the (packed) Hollerith string A to a multiple-precision number
X. A can be a string of digits acceptable to routine MPIN and
terminated by '$', e.g. '-5.367$' or '1E-99$', or one of the following
special strings
0 'EPS' (multiple-precision machine-precision, see MPEPS),
'EUL' (Euler's constant 0.5772..., see MPEUL),
'MAXR' (largest valid multiple-precision number, see MPMAXR),
'MINR' (smallest positive multiple-precision number, see MPMINR),
'PI' (pi = 3.14..., see MPPI).
0Actually, only the first two characters of these special strings are
significant. For example, 'MA', 'MAX' and 'MAXR' are equivalent.
0The error message '*** MXR TOO SMALL ... ***' is usually caused by
omission of the sentinel '$'.
0For the effect or rounding options, see the descriptions of MPIN, MPEPS,
MPMAXR, MPMINR and MPPI.
0Warning to Augment users - use CAM(A) and not CTM(A) if A is declared as
an integer array. (Otherwise Augment will generate a call to MPCIM
instead of MPCAM.)
0MPCDM
+_____
0CALL MPCDM (DX, Z) or Z = DX or Z = CTM (DX)
0Converts double-precision DX to multiple-precision Z. Some numbers will
not convert exactly on machines with base other than two, four or
sixteen, or if B is not a power of two, or if T is too small. Thus,
MPCDM should be used only to obtain starting approximations etc., and
not where full multiple-precision accuracy is required. Rounding
options are not implemented. For accurate initialization of
multiple-precision variables, use MPCAM, MPCIM, MPCQM or MPQPWR.
0
MPCDM, MPCMD, MPCMDE and MPCMPD are not called by any other routines in
the MP package, so they may be omitted if double-precision is not
available.
0
MPCEIL
+______
0CALL MPCEIL (X, Y) or Y = CEIL (X)
0Sets Y = ceiling (X), i.e. the smallest integer not less than X, where X
and Y are unpacked multiple-precision variables. Rounding is defined by
RNDRL in COMMON /MPCOM/, as for MPEXP (this is only relevant if X is
large and positive, for otherwise ceiling (X) is exactly representable
as a multiple-precision number).
1Sec. 6 MP User's Guide Page 31
+ _______________
0MPCHEB
+______
0CALL MPCHEB (C, NC, N, IND) or CALL MPCHEB (C, PARAM ('MT2'), N, IND)
0Converts the power series coefficients C(1), ... , C(N) to Chebyshev
series coefficients (which overwrite C(1), ... , C(N)). It is assumed
that on summation of the Chebyshev series the constant term is halved.
0IND may have the value 0, -1 or +1, with the following meanings -
0 0 - C(j) is the coefficient of x**(j-1) on input, and of T(j-1, x) on
output, where T(0, x), T(1, x), ... are the Chebyshev polynomials,
e.g.
T(3, x) = 4*x**3 - 3*x.
0 -1 - C(j) is the coefficient of x**(2j-1) on input, and of T(2j-1,
x) on output (useful for approximation of odd functions).
0 +1 - C(j) is the coefficient of x**(2j-2) on input, and of T(2j-2,
x) on output (useful for approximation of even functions).
0Conceptually, C is an array of multiple-precision variables. If MPCHEB
is called directly, C is a two-dimensional INTEGER array (with first
dimension NC at least T+2, second dimension at least N). If Augment is
used, C should be of type MULTIPLE, with dimension at least N, and NC
should be PARAM ('MT2').
0Rounding options are not implemented. Uses truncated arithmetic and no
guard digits, time = O(N*N*T), space = O(T).
0MPCHEV
+______
0CALL MPCHEV (C, NC, N, IND, X, Y)
0 or CALL MPCHEV (C, PARAM ('MT2'), N, IND, X, Y)
0Returns Y = Chebyshev series evaluated at X, where X and Y are unpacked
multiple-precision variables, C is an array containing the Chebyshev
coefficients. For a description of the arguments C, NC, N and IND, see
the notes on MPCHEB above.
0The algorithm used is described in Chapter 8 of - Modern Computing
Methods (Notes on Applied Science No. 16, HMSO, London, 1961).
TIme = O(N*m(T)), space = O(T).
0Rounding options are not implemented. Uses truncated arithmetic and no
guard digits.
0MPCHGB
+______
0J = MPCHGB (B1, B2, N)
0Returns J such that B1**abs(J) >= B2**abs(N), i.e.
0 abs(J) >= abs(N)*log(B2)/log(B1),
1Sec. 6 MP User's Guide Page 32
+ _______________
0and sign(J) = sign(N), where B1, B2 and J are integers, B1 > 1, B2 > 0,
B1*B2 <= MXINT. The value of J returned is close to minimal. Uses only
integer arithmetic.
0MPCHK
+_____
0CALL MPCHK
0Checks legality of parameters in COMMON /MPCOM/ (see Section 3.1), and
updates MXSPTR. If an illegal parameter in COMMON /MPCOM/ is detected,
an error message is written on unit LUN, and execution is terminated.
0MPCIM
+_____
0CALL MPCIM (I, Z) or Z = I or Z = CTM (I)
0Converts single-precision integer I to multiple-precision Z, where
abs(I) <= B**T.
0MPCIS
+_____
0CALL MPCIS (X, C, S, BOTH)
0Returns C = cos (X) and S = sin (X) if BOTH = .TRUE.,
or just C = cos (X) if BOTH = .FALSE.
0X, C and S are unpacked multiple-precision variables, BOTH is a logical
variable. Faster than calling MPCOS and MPSIN with the same first
argument. The algorithm is similar to that used in MPEXP (with
imaginary argument), time = O(sqrt(T)m(T)). Rounding options are
implemented as follows -
0 RNDRL = 0 or 1 - Absolute error < 0.6*B**(-T)
(but relative error may be large if abs(X) > 1).
RNDRL = 2 - Lower bound on the true result.
RNDRL = 3 - Upper bound on the true result.
0MPCMD
+_____
0CALL MPCMD (X, DZ) or DZ = X or DZ = CTD (X)
0Converts multiple-precision X to double-precision DZ. Assumes X is in
allowable range for double-precision numbers (otherwise double-precision
floating-point overflow or underflow may occur). There may be some loss
of accuracy if B is not a power of the base used for double-precision
floating-point arithmetic, or if T is too small. Rounding options are
not implemented. See also the description of MPCDM.
0MPCMDE
+______
0CALL MPCMDE (X, N, DX)
0Returns integer N and double-precision DX such that (multiple-precision)
0 X = DX*OUTBAS**N (approximately),
1Sec. 6 MP User's Guide Page 33
+ _______________
0where 1 <= abs(DX) < OUTBAS (unless X = 0). The default value of OUTBAS
is ten (see Section 3.1). It is assumed that X is not so large or small
that N overflows. X may be packed or unpacked. Rounding options are
not implemented. See also the description of MPCDM.
0
MPCMEF
+______
0CALL MPCMEF (X, N, Y)
0Given multiple-precision X, returns integer N and multiple-precision Y
such that
0 X = (OUTBAS**N)*Y and 1 <= abs(Y) < OUTBAS
0(unless X = 0, when N = Y = 0). The default value of OUTBAS is ten (see
Section 3.1). It is assumed that X is not so large or small that N
overflows.
0Rounding options are not fully implemented, but the directed rounding
options (RNDRL = 2 or 3) give correct (though not best possible) lower
and upper bounds on the true result.
0
MPCMF
+_____
0CALL MPCMF (X, Y) or Y = CMF (X) or Y = FRAC (X)
0For multiple-precision X and Y, returns fractional part of X in Y, i.e.
0 Y = X - integer part of X (truncated towards 0).
0Rounding options are irrelevant as the result is exact.
0
MPCMI
+_____
0CALL MPCMI (X, I) or I = X or I = CTI (X)
0Converts multiple-precision X to integer I, assuming that X is not too
large (else use MPCMIM). X is truncated towards zero (irrespective of
RNDRL). I is returned as zero if abs(int(X)) > MXINT. The user may
check for this possibility by testing if
0 ((X(1).NE.0) .AND. (X(2).GT.0) .AND. (I.EQ.0))
0is true on return from MPCMI.
0MPCMIM
+______
0CALL MPCMIM (X, Y) or Y = CMIM (X) or Y = INT (X)
0Returns Y = integer part of X (truncated towards 0), for
multiple-precision X and Y. Use if abs(Y) > MXINT (otherwise MPCMI is
preferable.) X may be packed or unpacked if MPCMIM is called directly.
1Sec. 6 MP User's Guide Page 34
+ _______________
0MPCMP
+_____
0J = MPCMP (X, Y)
0Compares the unpacked multiple-precision numbers X and Y, returning
0 +1 if X > Y,
0 if X = Y,
-1 if X < Y.
0X and Y must be unpacked. If they are packed, MPCOMP should be used
instead. The result is exact.
0
MPCMPA
+______
0J = MPCMPA (X, Y) or J = CMPA (X, Y)
0Compares abs(X) with abs(Y) for multiple-precision X and Y, returning
0 +1 if abs(X) > abs(Y),
0 if abs(X) = abs(Y),
-1 if abs(X) < abs(Y).
0X and/or Y may be packed or unpacked if MPCMPA is called directly, but
must be unpacked if called via Augment. The result is exact.
0
MPCMPD
+______
0J = MPCMPD (X, DI) or J = COMP (X, DI)
0Compares multiple-precision X with double-precision DI, returning
0 +1 if X > DI,
0 if X = DI,
-1 if X < DI.
0X may be packed or unpacked if MPCMPD is called directly. Uses MPCDM to
convert DI to multiple-precision, so the result may be incorrect if
MPCDM is inaccurate. See the description of MPCDM.
0
MPCMPI
+______
0J = MPCMPI (X, I) or J = COMP (X, I)
0Compares multiple-precision number X with integer I, returning
0 +1 if X > I,
0 if X = I,
-1 if X < I.
0X may be packed or unpacked if MPCMPI is called directly, but must be
unpacked if called via Augment. The result is exact.
1Sec. 6 MP User's Guide Page 35
+ _______________
0MPCMPQ
+______
0K = MPCMPQ (X, I, J) or K = COMP (X, I, J)
0Compares multiple-precision X with the rational number I/J, where I and
J are single-precision integers, J nonzero. Returns
0 +1 if X > I/J,
0 if X = I/J,
-1 if X < I/J.
0Rounding options are irrelevant as the result is exact (the rational
number I/J is never computed).
0
0MPCMPR
+______
0J = MPCMPR (X, RI) or J = COMP (X, RI)
0Compares multiple-precision number X with real number RI, returning
0 +1 if X > RI,
0 if X = RI,
-1 if X < RI.
0X may be packed or unpacked if MPCMPR is called directly. Uses MPCRM to
convert RI to multiple-precision, so the result may be incorrect if
MPCRM is inaccurate. See the description of MPCRM.
0
0MPCMR
+_____
0CALL MPCMR (X, RZ) or RZ = X or RZ = CTR (X)
0Converts multiple-precision X to single-precision real RZ. Assumes that
X is in the allowable range for single-precision real variables
(otherwise floating-point overflow or underflow may occur). There may
be some loss of accuracy if B is not a power of the base used for
single-precision real arithmetic. Rounding options are not implemented.
0
MPCMRE
+______
0CALL MPCMRE (X, N, RX)
0Returns integer N and single-precision real RX such that
multiple-precision
0 X = RX*OUTBAS**N (approximately),
0where 1 <= abs(RX) < OUTBAS (unless X = 0). It is assumed that X not so
large or small that N overflows. The default value of OUTBAS is ten
(see Section 3.1). Rounding options are not implemented.
1Sec. 6 MP User's Guide Page 36
+ _______________
0MPCMUL
+______
0CALL MPCMUL (XR, XI, YR, YI, ZR, ZI)
0Computes the complex product
0 ZR + i*ZI = (XR + i*XI) * (YR + i*YI),
0where i**2 = -1. XR, XI, YR, YI, ZR, ZI are unpacked multiple-precision
variables. Rounding options are not implemented. Uses truncated
arithmetic with no guard digits. Time = O(m(T)).
0MPCOMP
+______
0J = MPCOMP (X, Y) or J = COMP (X, Y)
0Compares the multiple-precision numbers X and Y, returning
0 +1 if X > Y,
0 if X = Y,
-1 if X < Y.
0X and/or Y may be packed or unpacked if MPCOMP is called directly, but
must be unpacked if called via Augment. The result is exact.
0MPCOS
+_____
0CALL MPCOS (X, Y) or Y = COS (X)
0Returns Y = cos(X) for multiple-precision X and Y, using MPCIS. X may
be packed or unpacked, Y is unpacked. Time = O(sqrt(T)m(T)). Rounding
options are implemented as for MPCIS, so the absolute error is small but
the relative error may be large.
0MPCOSH
+______
0CALL MPCOSH (X, Y) or Y = COSH (X)
0Returns multiple-precision Y = COSH(X) for multiple-precision X (not too
large). X may be packed or unpacked, Y is unpacked. Rounding options
are not yet implemented, uses no guard digits.
Time = O(sqrt(T)m(T)).
0MPCQM
+_____
0CALL MPCQM (I, J, Q) or Q = CQM (I, J) or Q = CTM (I, J)
0Converts the rational number I/J to multiple-precision Q. Rounding
options are implemented. Rounding is controlled by the parameter RNDRL
in COMMON /MPCOM/ (see Section 3.1), as follows -
0 0 - Round towards zero (i.e. chop or truncate). The resulting error is
less than 1 ulp (unit in the last place). This is the default
option if MPINIT or MPSET2 is used. It is the fastest option, and
is satisfactory for most purposes.
1Sec. 6 MP User's Guide Page 37
+ _______________
0
1 - Round to nearest, preferring even last digit if a tie occurs.
Error is less than or equal to 0.5 ulp. This option gives
worst-case error bounds half as large as for RNDRL = 0, and
average-case error considerably better than for RNDRL = 0.
However, it is expensive in time and space. If RNDRL = 0 does not
give sufficient accuracy, it is usually better to increase T than
to use RNDRL = 1.
0 2 - Round down (towards -infinity). Error is less than 1 ulp. This
option is useful for interval arithmetic.
0 3 - Round up (towards +infinity). Error is less than 1 ulp. This
option is useful for interval arithmetic.
0MPCRM
+_____
0CALL MPCRM (RX, Z) or Z = RX or Z = CTM (RX)
0Converts single-precision real RX to multiple-precision Z. Some numbers
will not convert exactly on machines with base other than two, four or
sixteen, or if B is not a power of 2, or if T is too small. Thus, MPCRM
should only be used to generate starting approximations etc., and not
where full multiple-precision accuracy is required. For accurate
initialization of multiple-precision variables, use MPCAM, MPCIM, MPCQM
or MPQPWR. Rounding options are not implemented.
0MPDAW
+_____
0CALL MPDAW (X, Y) or Y = DAW (X)
0Returns
0 Y = Dawson's integral of multiple-precision argument X
0 = exp(-X**2)*(integral from 0 to X of exp(u**2)du),
0where X and Y are multiple-precision variables. X may be packed or
unpacked if MPDAW is called directly. Rounding options are implemented
as for MPEXP.
0MPDGA
+_____
0I = MPDGA (X, N) or I = DIGIT (X, N)
0Returns the N-th digit of the multiple-precision number X, 0 < N <= T.
Returns zero if X is zero.
0MPDGB
+_____
0CALL MPDGB (I, X, N) or DIGIT (X, N) = I
0Sets the N-th digit of the multiple-precision number X to I. N must be
in the range 0 < N <= T, I must be in the range 0 <= I < B (and I > 0 if
N = 1). The sign and exponent of X are not changed.
1Sec. 6 MP User's Guide Page 38
+ _______________
0MPDIGA
+______
0I = MPDIGA (X) or I = NUMDIG (X)
0Returns T, the number of multiple-precision digits (the second word of
COMMON /MPCOM/). X is a dummy multiple-precision argument (Augment
expects one).
0MPDIGB
+______
0CALL MPDIGB (I, X) or NUMDIG (X) = I
0Sets T, the number of multiple-precision digits (second word of COMMON
/MPCOM/), to I, which should be an integer in the range 1 < I < MT2-1.
X is a dummy multiple-precision argument (Augment expects one).
0MPDIGS
+______
0J = MPDIGS (N)
0Returns the number of multiple-precision (base B) digits required for
the equivalent of at least N floating 'decimal' places. The inequality
0 B**(MPDIGS - 1) >= OUTBAS**(N - 1)
0will be satisfied, usually with the minimal such value of MPDIGS. If
OUTBAS has its default value (ten), the working precision may be set to
the equivalent of at least N floating decimal places by the statement
0 PARAM ('NUMDIG') = MPDIGS (N) (using Augment)
or
CALL MPPARB (MPDIGS (N), 'NUMDIG') (not using Augment).
0MPDIGV
+______
0J = MPDIGV (C)
0Returns the integer value
0 0 , ... , 9 , 10, ... , 15 if C is the corresponding character
'0', ... , '9', 'A', ... , 'F'
0and the value returned is less than INBASE (default value ten, see
Section 3.1). Otherwise returns -1. The character value should be
represented in an integer variable as if read under A1 format.
0MPDIGW
+______
0C = MPDIGW (J)
0Returns the character '0', ... , 'F' if J is the corresponding integer,
0 , ... , 15,
0returns '*' otherwise. The character is stored in an integer as if read
under A1 format.
1Sec. 6 MP User's Guide Page 39
+ _______________
0MPDIM
+_____
0CALL MPDIM (X, Y, Z) or Z = DIM (X, Y)
0Sets Z = X - min (X, Y) = max (0, X-Y) for multiple-precision variables
X, Y and Z. X and Y may be both packed or both unpacked, Z is unpacked.
Rounding options are implemented as for MPSUB.
0
MPDIV
+_____
0CALL MPDIV (X, Y, Z) or Z = X/Y
0Sets Z = X/Y, for multiple-precision X, Y and Z. Y must not be zero.
Uses MPDIVL for small T or if RNDRL > 0, MPREC and MPMUL otherwise.
Thus time = O(m(T) + RNDRL*T**2). Rounding options are implemented as
for MPCQM.
0
MPDIVI
+______
0CALL MPDIVI (X, IY, Z) or Z = X/IY
0Divides unpacked multiple-precision X by the nonzero single-precision
integer IY, giving multiple-precision Z in time O(T). This is much
faster than division by a multiple-precision number. Rounding options
are implemented as for MPCQM.
0MPDIVL
+______
0CALL MPDIVL (X, Y, Z)
0Sets Z = X/Y for unpacked multiple-precision variables X, Y and Z
(Y nonzero). Uses the 'schoolboy' or 'long division' algorithm, so
time = O(T**2). An alternative method is to use MPREC and MPMUL. This
would be faster for large T, but MPDIVL is faster for small T. Rounding
options are implemented as for MPCQM.
0MPDIV2
+______
0Routine called by MPDIVI, not recommended for independent use.
0MPDIV3
+______
0Routine called by MPDIVL, not recommended for independent use.
0MPDUMP
+______
0CALL MPDUMP (X)
0Dumps the (packed or unpacked) multiple-precision number X on unit LUN.
Writes the sign, exponent and digits in suitable I format (or just the
sign if X is zero). Embedded blanks should be interpreted as zeros.
MPDUMP is sometimes useful for debugging, but is not recommended in
general for output of MP numbers.
1Sec. 6 MP User's Guide Page 40
+ _______________
0MPEI
+____
0CALL MPEI (X, Y) or Y = EI (X)
0Returns
0 Y = ei (X) = -E1 (X)
= (principal value integral from -infinity to X of exp(u)/u du),
0for multiple-precision numbers X and Y, using the power series for small
abs(X), the asymptotic series for large abs (X), and the continued
fraction for medium negative X. The relative error in Y is small unless
X is very close to the zero 0.37250741078136663446... of ei(X), and then
the absolute error in Y is O(B**(1-T)). In any case the error in Y
could be induced by an O(B**(1-T)) relative perturbation in X. Rounding
options are not yet implemented. Time = O(m(T)T).
0
0MPEPS
+_____
0CALL MPEPS (X) or X = CTM ('EPS')
0Sets multiple-precision X to the (multiple-precision) machine precision,
that is an upper bound on the smallest representable positive number X
such that the relative error in the basic multiple-precision operations
of addition, subtraction, multiplication and division does not exceed X
unless the result underflows. In fact
0 X = 1.01*B**(1-T) (rounded up) if RNDRL = 0,
= 0.5*B**(1-T) (rounded up) if RNDRL = 1,
= B**(1-T) if RNDRL = 2 or 3.
0
0MPEQ
+____
0LV = MPEQ (X, Y) or LV = (X .EQ. Y) or IF (X .EQ. Y) ...
0Returns logical value LV of (X .EQ. Y) for multiple-precision X and Y.
MPEQ must be declared LOGICAL unless the Augment interface is used.
X and/or Y may be packed or unpacked.
0
0MPERF
+_____
0CALL MPERF (X, Y) or Y = ERF (X)
0Returns Y = erf(X)
= sqrt(4/pi)*(integral from 0 to X of exp(-u**2) du)
0for multiple-precision X and Y. X may be packed or unpacked if MPERF is
called directly, must be unpacked if MPERF is called via Augment.
Rounding options are implemented as for MPEXP.
1Sec. 6 MP User's Guide Page 41
+ _______________
0MPERFC
+______
0CALL MPERFC (X, Y) or Y = ERFC (X)
0Returns Y = erfc(X) = 1 - erf(X) for multiple-precision numbers X and Y.
X may be packed or unpacked if MPERFC is called directly, must be
unpacked if MPERFC is called via Augment. Rounding as for MPEXP.
0MPERF2
+______
0Computes exp(X**2)*(integral from 0 to X of exp(-u*u) du)
0for multiple-precision X, using the power series for small X, and MPEXP
for large X. Called by MPERF, not recommended for independent use.
Rounding options are not implemented, no guard digits are used.
0MPERF3
+______
0Trys to return
0 exp(X**2)*(integral from X to infinity of exp(-u**2)du),
or
exp(-X**2)*(integral from 0 to X of exp(u**2) du),
0using the asymptotic expansion, where X is a multiple-precision
variable. The condition for success is approximately that
0 X > sqrt(T*ln(B)).
0Called by MPERF, MPERFC and MPDAW, not recommended for independent use.
Rounding options are not implemented, uses no guard digits.
0MPERR
+_____
0CALL MPERR
0This routine is called when a fatal error condition is encountered, and
after a message has been written on logical unit LUN. As supplied in
the multiple-precision package, MPERR writes a message which includes
the version number in the format YYMMDD, then stops. On many systems it
can easily be modified to give a trace-back (e.g. with Univac 1100
Fortran V this can be obtained by replacing STOP by RETURN 0).
0Since MPERR is only called when a fatal error condition has been
detected, it does not return to the calling routine.
0MPERRM
+______
0CALL MPERRM (A)
0The argument A is a packed Hollerith string terminated by '$'. MPERRM
writes the string on unit LUN, then calls MPERR. Only the first 71
characters of A are significant. They will be preceeded and followed by
'***'. MPERRM is called by MP routines when a fatal error is
encountered, and does not return to the calling routine.
1Sec. 6 MP User's Guide Page 42
+ _______________
0MPEUL
+_____
0CALL MPEUL (G) or G = CTM ('EUL')
0Returns multiple-precision G = Euler's constant (gamma = 0.57721566...)
to almost full multiple-precision accuracy. The method was discovered
by Edwin McMillan and Richard Brent, and is faster than the method of
Sweeney (used in earlier versions of MPEUL). See - R. P. Brent and E.
M. McMillan, 'Some new algorithms for high-precision computation of
Euler's constant', Math. Comp. 34(1980), 305-312. Time = O(T**2).
Rounding options are implemented as for MPEXP.
0MPEXP
+_____
0CALL MPEXP (X, Y) or Y = EXP (X)
0Returns Y = exp(X) for multiple-precision X and Y. X may be packed or
unpacked if MPEXP is called directly. The method used is described in -
R. P. Brent, 'The complexity of multiple-precision arithmetic', in The
Complexity of Computational Problem Solving, Queensland University
Press, Brisbane, 1976, 126-165. See also - R. P. Brent, 'Unrestricted
algorithms for elementary and special functions', Information Processing
80, North-Holland, 1980, 613-619. Time = O(sqrt(T)m(T)).
0Rounding is controlled by the parameter RNDRL in COMMON /MPCOM/ (see
Section 3.1), as follows -
0 0 or 1 - error less than 0.6 ulp (units in the last place),
0 2 - computed result is lower bound on the true result,
0 3 - computed result is upper bound on the true result.
0
MPEXPA
+______
0I = MPEXPA (X) or I = EXPON (X)
0Returns the exponent of the multiple-precision number X (or large
negative exponent -M if X is zero). X may be packed or unpacked if
MPEXPA is called directly.
0MPEXPB
+______
0CALL MPEXPB (I, X) or EXPON (X) = I
0Sets the exponent of the multiple-precision number X to I, unless X is
zero (when its exponent is unchanged). X must be a valid
multiple-precision number (either zero or normalized). X may be packed
or unpacked if MPEXPB is called directly. Note that multiple-precision
underflow/overflow will occur if I is too small/large.
0MPEXP1
+______
0CALL MPEXP1 (X, Y) or Y = EXP1 (X)
1Sec. 6 MP User's Guide Page 43
+ _______________
0Returns
0 Y = exp(X) - 1
0where X and Y are multiple-precision numbers and -1 < X < 1. Uses an
O(sqrt(T)m(T) algorithm described in - R. P. Brent, 'The complexity of
multiple-precision arithmetic', in Complexity of Computational Problem
Solving, Univ. of Queensland Press, Brisbane, 1976, 126-165. Rounding
options are implemented as for MPEXP.
0MPFIN
+_____
0CALL MPFIN (LUNIT, X) or X = FIO (LUNIT)
0Reads (unpacked) multiple-precision number X from Fortran logical unit
LUNIT (LUNIT >= 0). It is assumed that a record of length INRECL <= 80
may be read in '(80A1)' format from unit LUNIT, and then converted to
multiple-precision using MPIN. Thus, MPFIN may be used to read
free-format fixed or floating-point numbers, one per record, from unit
LUNIT. The input base INBASE has default value ten (i.e. decimal input
- see Section 3.1). For further details see under MPIN.
0MPFLOR
+______
0CALL MPFLOR (X, Y) or Y = FLOOR (X)
0Sets Y = floor (X), i.e. the largest integer not exceeding X, where X
and Y are unpacked multiple-precision variables. Rounding is defined by
RNDRL in COMMON /MPCOM/, as for MPEXP (this is only relevant if X is
large and negative, for otherwise floor (X) is exactly representable as
a multiple-precision number).
0
MPFOUT
+______
0CALL MPFOUT (X, N) or FIO (N) = X
0Writes the multiple-precision number X on unit LUN in floating-point
format, with N significant 'decimal' digits (N > 1). The exponent field
width is EXWID (default is 6, including exponent character EXPCH and
sign - see Section 3.1). The default exponent character EXPCH is 'E'
(or '$' if OUTBAS = 15 or 16), and the default output base OUTBAS is ten
(i.e. decimal output). Rounding options are implemented as for MPOUTE.
X may be packed or unpacked if MPFOUT is called directly, but must be
unpacked if called via Augment.
0MPGAM
+_____
0CALL MPGAM (X, Y) or Y = GAM (X)
0Computes multiple-precision Y = Gamma(X) for multiple-precision argument
X, using MPGAMQ if 240*X is a small integer, otherwise using MPLNGM.
Space required is about the same as for MPLNGM, i.e. O(T**2),
time = O(T**3). Rounding options are not yet implemented and no guard
digits are used.
1Sec. 6 MP User's Guide Page 44
+ _______________
0MPGAMQ
+______
0CALL MPGAMQ (I, J, X) or Y = GAMQ (I, J) or Y = GAM (I, J)
0Returns
0 X = Gamma (I/J),
0where X is multiple-precision, I and J are small integers. The method
used is reduction of the argument to (0, 1) and then a direct expansion
of the defining integral truncated at a sufficiently high limit, using
2T digits to compensate for cancellation. Time = O(T**2) if I/J is not
too large. If I/J > 100 (approximately) it is faster to use MPGAM.
(MPGAMQ is very slow if I/J is is very large, because the relation
Gamma(X+1) = X*Gamma(X) is used repeatedly.) Rounding options are not
yet implemented.
0MPGCD
+_____
0CALL MPGCD (K, L)
0Returns K = K/GCD and L = L/GCD, where GCD is the greatest common
divisor of the initial K and L (single-precision integers). Called by
various MP routines.
0
MPGCDA
+______
0CALL MPGCDA (X, Y, Z) or Z = GCD (X, Y)
0Returns Z = greatest common divisor of X and Y. (By definition,
GCD (X, 0) = GCD (0, X) = abs(X), GCD (X, Y) >= 0.) X, Y and Z are
integers represented as multiple-precision numbers, and must satisfy
abs(X) < B**T, abs(Y) < B**T. The method is a straight-forward
implementation of the Euclidean algorithm, and Lehmer's trick is not
used, although it might be faster (see Amer. Math. Monthly 45 (1938),
227-233). Time = O(T**2). X and/or Y may be packed or unpacked if
MPGCDA is called directly. Rounding options are irrelevant as the
result is exact.
0MPGCDB
+______
0CALL MPGCDB (X, Y)
0Returns (X, Y) as (X/Z, Y/Z) where Z is the GCD of initial X and Y,
which are integers represented as multiple-precision numbers, and must
satisfy abs(X) < B**T, abs(Y) < B**T. Time = O(T**2). Rounding options
are irrelevant as the result is exact.
0
MPGD
+____
0J = MPGD (N)
0Returns ceiling (ln (max (1, abs(N))) / ln(B)),
1Sec. 6 MP User's Guide Page 45
+ _______________
0i.e. the minimum J >= 0 such that
0 B**J >= abs(N).
0This function is useful for computing the number of guard digits
required for various (multiple-precision) base B calculations.
0MPGD3
+_____
0CALL MPGD3 (N, TG)
0Sets T = T + 1 + MPGD (100*N) if it is safe to call MPGD, and the same
or slightly more otherwise. Also sets TG to the new value of T. N and
TG are integer arguments. Called by various MP routines when it is
necessary to increase the working precision.
0MPGE
+____
0LV = MPGE (X, Y) or LV = (X .GE. Y) or IF (X .GE. Y) ...
0Returns the logical value of (X > Y) for multiple-precision X and Y.
MPGE must be declared LOGICAL unless the Augment interface is used.
X and/or Y may be packed or unpacked.
0MPGET
+_____
0J = MPGET (X, N)
0Returns X(N), where X is an integer array of dimension at least N.
Necessary to avoid some compiler diagnostics.
0MPGT
+____
0LV = MPGT (X, Y) or LV = (X .GT. Y) or IF (X .GT. Y) ...
0Returns the logical value of (X > Y) for multiple-precision X and Y.
MPGT must be declared LOGICAL unless the Augment interface is used. X
and/or Y may be packed or unpacked.
0MPHANK
+______
0Tries to compute the Bessel function J (NU, X) using Hankel's asymptotic
expansion, for nonnegative integer NU and multiple-precision X. Time =
O(T**3). Called by MPBESJ, not recommended for independent use.
Rounding options are not implemented, uses no guard digits.
0MPIMUL
+______
0CALL MPIMUL (IY, X, Z) or Z = IY*X
0Equivalent to
0 CALL MPMULI (X, IY, Z) or Z = X*IY.
0See the description of MPMULI for further details.
1Sec. 6 MP User's Guide Page 46
+ _______________
0MPIN
+____
0CALL MPIN (C, X, N, ERROR)
0Converts the characters (assumed to have been read under NA1 format) in
C(1), ... , C(N) to a multiple-precision number in X. If C represents a
valid number, ERROR is returned as .FALSE. If C does not represent a
valid number, ERROR is returned as .TRUE. and X as zero. Leading and
trailing blanks are allowed, embedded blanks (except between the number
and its sign) are forbidden. If there is no decimal point one is
assumed to lie just to the right of the last decimal digit. X is a
multiple-precision number, C an integer array, N an integer, and ERROR
logical.
0The input base is determined by the parameter INBASE in COMMON /MPCOM/
(see Section 3.1). INBASE has default value ten, but may be set to any
value in two, ... , sixteen.
0 Examples of valid numbers Examples of invalid numbers
+ _________________________ ___________________________
0 - 123456789 12 345
3.14159 123.456E -67
-44. 1.2.3
.0001234 E123
123.456D789 64.4E+
-.1234566-789 ++12.3
+999+88 E3.
0Rounding is determined by the parameter RNDRL in COMMON /MPCOM/ as
follows -
0 0 or 1 - Round to (approximately) the nearest representable base
B, T-digit number. The error is less than 0.6 units in
the last (base B) place.
0 2 - Round down (towards -infinity) to a base B, T-digit
number,
0 3 - Round up (towards +infinity) to a base B, T-digit number.
0MPINE
+_____
0CALL MPINE (C, X, N, J, ERROR)
0MPINE is the same as MPIN except that the result (X) is multiplied by
INBASE**J, where J is a single-precision integer. For details of the
other arguments, see MPIN. Useful for floating-point input of
multiple-precision numbers. The user can read the exponent into J
(using any suitable format) and the fraction into C (using NA1 format),
then call MPINE to convert to multiple-precision.
0Note - in early versions of the MP package, MPIN was unable to deal with
numbers with exponents (e.g. 5.3E-26). Now that MPIN can deal with such
numbers there is no real need for MPINE, so it is included only for
compatibility with early versions of the package.
1Sec. 6 MP User's Guide Page 47
+ _______________
0MPINF
+_____
0CALL MPINF (X, N, UNIT, IFORM, ERR) or ERR = MPINF (X, N, UNIT, IFORM)
0 or IF (MPINF (X, N, UNIT, IFORM)) ...
0Reads N words from logical unit abs(UNIT) using format IFORM, then
converts to multiple-precision number X using routine MPIN. IFORM
should contain a format which allows for reading N words in A1 format,
e.g. '(80A1)'.
0ERR is returned as .TRUE. if MPIN could not interpret input as a
multiple-precision number or if N not positive, otherwise .FALSE. In
the former case X is returned as zero.
0Note - for a more convenient, though less flexible, method of reading
multiple-precision numbers, see MPFIN.
0MPINIT
+______
0CALL MPINIT (MP) or INITIALIZE MP
0Declares blank COMMON (see Section 3.1) and calls MPSET2 to initialize
parameters. MP is a dummy integer argument. The Augment declaration
0 INITIALIZE MP
0causes the statement
0 CALL MPINIT (MP)
0to be generated.
0Warning - as distributed MPINIT assumes output unit LUN = 6, at most 25
multiple-precision digits (i.e. MT2 = 27), MNSPTR = 1, and MXR = 500.
If the Augment description deck is changed, MPINIT should be changed
accordingly (see Section 9.2).
0MPINTG
+______
0LV = MPINTG (X) or LV = INTEGR (X) or IF (INTEGR (X)) ...
0Returns .TRUE. if the (unpacked) multiple-precision number X is an exact
integer, .FALSE. otherwise. LV is a LOGICAL variable. MPINTG must be
declared LOGICAL if it is called directly.
0MPIO
+____
0CALL MPIO (C, N, UNIT, IFORM, ERR)
0Writes C(1), ... , C(N) in format IFORM if UNIT > 0,
Reads C(1), ... , C(N) in format IFORM if UNIT <= 0,
0in both cases using logical unit abs(UNIT). Unformatted I/O is
performed if IFORM = 'U'. C is an integer array of dimension at least N.
1Sec. 6 MP User's Guide Page 48
+ _______________
0N and UNIT are integers, ERR is a logical variable. ERR is returned
as .TRUE. iff N <= 0. (It would be desirable to return ERR = .TRUE. if
an I/O error occurred. This is easily done on most systems, but can not
be done with ANSI Standard Fortran. See Section 9.3.)
0MPIS
+____
0LV = MPIS (J, K)
0Assumes that J and K are INTEGER words each containing one character in
the leftmost character position (read under A1 format or set by a 1Hx
data statement or unpacked by MPUPK). MPIS returns .TRUE. if the two
characters are equal, .FALSE. otherwise. LV is a LOGICAL variable, and
MPIS must be declared to be LOGICAL.
0On most machines, MPIS (J, K) is equivalent to (J .EQ. K), but on the
Burroughs B6700 it is equivalent to (J .IS. K). For conversion notes,
see Section 9.2.
0MPKADD
+______
0CALL MPKADD (X, Y, Z) or Z = X + Y
0The same as CALL MPADD (X, Y, Z) except that X and Y are packed
multiple-precision variables. (They may be packed or unpacked if MPKADD
is called directly - similarly for MPKDIV, MPKDVI, MPKIML, MPKMLI,
MPKMUL and MPKSUB.) For further details, see under MPADD.
0MPKDIV
+______
0CALL MPKDIV (X, Y, Z) or Z = X/Y
0The same as CALL MPDIV (X, Y, Z) except that X and Y are packed
multiple-precision variables. For further details see under MPDIV.
0MPKDVI
+______
0CALL MPKDVI (X, IY, Z) or Z = X/IY
0The same as CALL MPDIVI (X, IY, Z) except that X is a packed
multiple-precision variable. For further details, see under MPDIV.
0MPKIML
+______
0CALL MPKIML (IY, X, Z) or Z = IY*X
0The same as CALL MPMULI (X, IY, Z) except that X is a packed
multiple-precision variable. For further details, see under MPMULI.
0MPKMLI
+______
0CALL MPKMLI (X, IY, Z) or Z = X*IY
0The same as CALL MPMULI (X, IY, Z) except that X is a packed
multiple-precision variable. For further details, see under MPMULI.
1Sec. 6 MP User's Guide Page 49
+ _______________
0MPKMUL
+______
0CALL MPKMUL (X, Y, Z) or Z = X*Y
0The same as CALL MPMUL (X, Y, Z) except that X and Y are packed
multiple-precision variables. For further details, see under MPMUL.
0MPKSUB
+______
0CALL MPKSUB (X, Y, Z) or Z = X - Y
0The same as CALL MPSUB (X, Y, Z) except that X and Y are packed
multiple-precision variables. For further details, see under MPSUB.
0MPK3V
+_____
0CALL MPK3V (MPXX, X, Y, Z)
0Calls MPXX (XU, YU, Z) after unpacking X and Y to give XU and YU. X and
Y are packed or unpacked multiple-precision variables, Z is an unpacked
multiple-precision variable, MPXX is a subroutine taking three unpacked
multiple-precision arguments and returning a result as the third
argument (e.g. MPMUL). Called by MPKADD, MPKMUL, MPKSUB etc.
0MPK3V2
+______
0CALL MPK3V2 (MPXX, X, N, Y)
0Calls MPXX (XU, N, Z) after unpacking X to give XU. X is a packed or
unpacked multiple-precision variable, N is an integer, Z is an unpacked
multiple-precision variable, MPXX is a subroutine taking three arguments
of the correct types and returning a result as the third argument (e.g.
MPMULI). Called by MPKDVI, MPKMLI etc.
0
MPLARG
+______
0CALL MPLARG (MXINT)
0Returns MXINT <= maximum representable INTEGER of the form 2**k - 1.
Integer arithmetic must be performed exactly on integers in the range
-MXINT, ... , +MXINT, so on some machines (e.g. Burroughs B6700 and
Cyber 76) MXINT must be smaller than the largest representable integer.
We say that the 'effective' wordlength is k+1 bits.
0Note - The version of MPLARG supplied with the MP package does not work
on all machines. For conversion notes, see Section 9.2.
0MPLE
+____
0LV = MPLE (X, Y) or LV = (X .LE. Y) or IF (X .LE. Y) ...
0Returns logical value of (X < Y) for multiple-precision X and Y. MPLE
must be declared type logical unless Augment interface used. X and/or Y
may be packed or unpacked.
1Sec. 6 MP User's Guide Page 50
+ _______________
0MPLG10
+______
0CALL MPLG10 (X, Y) or Y = LOG10 (X)
0Returns Y = ln(X)/ln(10) = log(X) to base 10, for multiple-precision X
and Y. X may be packed or unpacked, Y is unpacked. Uses MPLN and
MPLNI. For restrictions, time bounds, and the effect of rounding
options, see the description of MPLN.
0
MPLI
+____
0CALL MPLI (X, Y) or Y = LI (X)
0Returns
0 Y = Li(X) = logarithmic integral of X
= principal value integral from 0 to X of du/ln(u),
0using MPEI. X and Y are multiple-precision numbers, X > 0, X .NE. 1.
X may be packed or unpacked if MPLI is called directly. Rounding
options are not yet implemented. Error in Y could be induced by an
O(B**(1-T)) relative perturbation in X followed by similar perturbation
in Y. Thus relative error in Y is small unless X is close to 1 or to
the zero 1.45136923488338105028... of Li(X). Time = O(m(T)T).
0
MPLN
+____
0CALL MPLN (X, Y) or Y = LN (X) or Y = LOG (X) or Y = ALOG (X)
0Returns Y = ln(X), for multiple-precision X and Y, using MPLNS. X may
be packed or unpacked, Y is unpacked. The integer part of ln(X) must be
representable as a single-precision integer. Time = O(sqrt(T)m(T)).
For small integer X, MPLNI is faster. Asymptotically faster methods
exist (e.g. the Gauss-Salamin method), but are not useful unless T is
large (see the descriptions of MPATAN, MPEXP1, MPLNGS and MPPIGL).
Rounding options are implemented as for MPEXP.
0
MPLNGM
+______
0CALL MPLNGM (X, Y) or Y = LNGM (X)
0Returns multiple-precision Y = ln(Gamma(X)) for positive
multiple-precision X, using Stirling's asymptotic expansion. (X may be
packed or unpacked if MPLNGM is called directly.) Slower than MPGAMQ
(unless X large) and uses more space, so use MPGAMQ and MPLN if X is
rational and not too large, say X less than 100. Time = O(T**3).
0Space = nl*((T+3)/2) + O(T) words, where nl is the number of terms used
in the asymptotic expansion, nl < 2+0.125*T*ln(B). Thus space = O(T**2)
words.
0Rounding options are not yet implemented, uses no guard digits.
1Sec. 6 MP User's Guide Page 51
+ _______________
0MPLNGS
+______
0CALL MPLNGS (X, Y)
0Returns Y = ln(X) for multiple-precision X and Y, using the
Gauss-Salamin algorithm based on the arithmetic-geometric mean
iteration. This is described in - R. P. Brent, 'Multiple-precision
zero-finding methods and the complexity of elementary function
evaluation', in Analytic Computational Complexity (ed. by J.F. Traub),
Academic Press, 1976, 151-176. (If abs(X-1) < 1/B, MPLNS is used.)
Time = O(ln(T)m(T)) + O(T**2) if abs(X-1) > 1/B, as for MPLNS otherwise.
0MPLNGS is slower than MPLN unless T is large (greater than about 500),
so it is recommended for testing purposes only. Rounding options are
implemented as for MPEXP.
0
MPLNI
+_____
0CALL MPLNI (N, X) or X = LNI (N) or X = ALOG (N)
or X = LN (N) or X = LOG (N)
0Returns multiple-precision X = ln(N) for small positive integer N, using
a rapidly converging series and MPL235. Time = O(T**2), faster than
MPLN. Rounding options are implemented as for MPEXP.
0
MPLNS
+_____
0CALL MPLNS (X, Y) or Y = LNS (X)
0Returns multiple-precision Y = ln(1+X) if X is a multiple-precision
number satisfying the condition abs(X) < 1/B, error otherwise. X may
be packed or unpacked if MPLNS is called directly. Uses Newton's method
to solve the equation exp1(-z) = X, then sets Y = -z. (Here exp1(z) =
exp(z) - 1 is computed using MPEXP1.) Time = O(sqrt(T)m(T)). Rounding
options are implemented as for MPEXP.
0
MPLT
+____
0LV = MPLT (X, Y) or LV = (X .LT. Y) or IF (X .LT. Y) ...
0Returns logical value of (X < Y) for multiple-precision X and Y. MPLT
must be declared LOGICAL unless the Augment interface used. X and/or Y
may be packed or unpacked.
0
MPL235
+______
0Computes ln ((2**i)*(3**j)*(5**k)) = i*ln(2) + j*ln(3) + k*ln(5) for
small integer i, j and k. ln(81/80), ln(25/24) and ln(16/15) are
calculated first, using rapidly converging series. Time = O(T**2).
Called by MPLNI, not recommended for independent use. Rounding options
are not implemented, uses no guard digits.
1Sec. 6 MP User's Guide Page 52
+ _______________
0MPMAX
+_____
0CALL MPMAX (X, Y, Z) or Z = MAX (X, Y)
0Sets Z = max (X, Y) where X, Y and Z are multiple-precision variables.
X and/or Y may be packed or unpacked if MPMAX is called directly.
Rounding options are irrelevant as the result is exact.
0MPMAXR
+______
0CALL MPMAXR (X) or X = CTM ('MAXR')
0Sets X to the largest possible positive multiple-precision number, i.e.
B**M - B**(M-T). Rounding options are irrelevant as the result is
exact.
0MPMEXA
+______
0I = MPMEXA (X) or I = MAXEXP (X)
0Returns the maximum allowable exponent of multiple-precision numbers
(the third word of COMMON /MPCOM/ - see Section 3.1). X is a dummy
multiple-precision argument (Augment expects one).
0MPMEXB
+______
0CALL MPMEXB (I, X) or MAXEXP (X) = I
0Sets the maximum allowable exponent of multiple-precision numbers to I.
I should be greater than T, and 4*I should be representable as a
single-precision integer. X is a dummy multiple-precision argument
(Augment expects one).
0MPMIN
+_____
0CALL MPMIN (X, Y, Z) or Z = MIN (X, Y)
0Sets Z = min (X, Y) where X, Y and Z are multiple-precision variables.
X and/or Y may be packed or unpacked if MPMIN is called directly.
Rounding options are irrelevant as the result is exact.
0MPMINR
+______
0CALL MPMINR (X) or X = CTM ('MINR')
0Sets X to the smallest positive normalized multiple-precision number,
i.e. X = B**(-M), where M is the third word of COMMON /MPCOM/ (see
Section 3.1). Rounding options are irrelevant as the result is exact.
0MPMLP
+_____
0Performs inner multiplication loop for MPMUL. Carries are not
propagated in inner loop, which saves time at the expense of space.
This routine is usually called more often than any other MP routine, so
it is worthwhile to optimize it.
1Sec. 6 MP User's Guide Page 53
+ _______________
0MPMOD
+_____
0CALL MPMOD (X, Y, Z) or Z = MOD (X, Y)
0Returns Z = X - Y*int(X/Y) for multiple-precision X, Y and Z. X and Y
may be both packed or both unpacked, Z is unpacked. Here 'int' means
the integer part, truncated towards zero.
0Notes - 1. MPMOD returns Z as zero if Y is zero.
2. Time = O (max (1, ln(abs(X/Y))) * T**2)),
which is large if abs(X) >> abs(Y).
3. Rounding options are irrelevant as the result is exact.
0
MPMOVE
+______
0CALL MPMOVE (X, TX, Y, TY)
0Assumes that X and Y are multiple-precision numbers with TX and TY
digits respectively. X may be packed or unpacked, Y is unpacked (but
not necessarily initialized). MPMOVE moves X to Y, padding with
trailing zeros if TX < TY, or rounding as specified by RNDRL (see the
description of MPCQM) if TX > TY. Note that TX and TY are integer
arguments, and the value of T (second word of COMMON /MPCOM/) is
irrelevant.
0MPSTR (X, Y) is equivalent and faster if TX = TY = T and X is unpacked.
0
0MPMUL
+_____
0CALL MPMUL (X, Y, Z) or Z = X*Y
0Multiplies X and Y, returning result in Z, for multiple-precision
variables X, Y and Z. The simple O(T**2) algorithm is used, with at
least two guard digits. Advantage is taken of any zero digits in X, but
not in Y. Asymptotically faster algorithms are known (see for example
Knuth, Vol. 2), but are difficult to implement in Fortran in an
efficient and machine-independent manner.
0In description of other multiple-precision routines, m(T) is the time to
perform T-digit multiple-precision multiplication. Thus m(T) = O(T**2)
with the present version of MPMUL, but m(T) = O(T.log(T)log(log(T))) is
theoretically possible (see Knuth, Vol. 2).
0Rounding options are implemented. They are controlled by the parameter
RNDRL in COMMON /MPCOM/ (see Section 3.1) as follows -
0 RNDRL = 0 - truncate towards zero, error less than 1.01 units in the
last place (ulp),
0 RNDRL = 1 - round to nearest representable multiple-precision number,
preferring last digit even in case of a tie, error at most
0.5 ulp (unless the result underflows),
1Sec. 6 MP User's Guide Page 54
+ _______________
0 RNDRL = 2 - round down (towards -infinity), error less than 1 ulp,
0 RNDRL = 3 - round up (towards +infinity), error less than 1 ulp.
0Note that RNDRL = 0 is the fastest (approximately twice as fast as the
other options) because less guard digits are required. Thus, RNDRL = 0
is recommended for general use and is the default option.
0
MPMULI
+______
0CALL MPMULI (X, IY, Z) or Z = X*IY
0Multiplies multiple-precision X by single-precision integer IY giving
multiple-precision Z. Time = O(T), which is faster than MPMUL.
Multiplication by 1 may be used to normalize a number even if the last
digit is B. Rounding options are implemented as for MPCQM.
0
MPMULQ
+______
0CALL MPMULQ (X, I, J, Y) or Y = MULQ (X, I, J)
0Multiplies multiple-precision X by rational number I/J, giving
multiple-precision Y. Here I and J are single-precision integers, J
nonzero. Time = O(T).
0Rounding options are determined by RNDRL in COMMON /MPCOM/ (see Section
3.1) as follows -
0 RNDRL = 0 - as for MPMULI followed by MPDIVI,
0 RNDRL = 1, 2, or 3 - as for MPCQM.
0Augment users should note that Y = MULQ (X, I, J) is usually faster than
Y = X * CTM (I,J).
0
MPMULS
+______
0CALL MPMULS (X, I, J, K, L)
0Sets X = X*I*J/(K*L) for unpacked multiple-precision X, integer I, J,
K and L (K*L nonzero). Calls MPMULQ once if I*J and K*L are not too
large, otherwise calls MPMULQ twice. Rounding is not best possible, but
the directed rounding options give correct upper and lower bounds.
0
MPNE
+____
0LV = MPNE (X, Y) or LV = (X .NE. Y) or IF (X .NE. Y) ...
0Returns logical value of (X .NE. Y) for multiple-precision X and Y.
MPNE must be declared LOGICAL unless the Augment interface used. X
and/or Y may be packed or unpacked.
1Sec. 6 MP User's Guide Page 55
+ _______________
0MPNEG
+_____
0CALL MPNEG (X, Y) or Y = -X
0Sets Y = -X for multiple-precision numbers X and Y. X and Y may be both
packed or both unpacked. (This violates the usual rule of returning an
unpacked result regardless of whether the argument is packed or
unpacked. The reason is that Augment does not allow unary operators to
return a type different from that of the variable which they operate
on.)
0
0MPNEW
+_____
0CALL MPNEW (I)
0Returns integer index I such that words I, ... , I+T+1 of blank COMMON
are available for use. Sets SPTR = I+T+2 and updates MXSPTR (see
Section 3.1). MPSTOV is called if MXR is too small.
0Note - CALL MPNEW (I) is equivalent to CALL MPNEW2 (I, T+2).
0
MPNEW2
+______
0CALL MPNEW2 (I, J)
0Returns integer index I such that words I, ... , I+abs(J)-1 of blank
COMMON are available for use. Sets SPTR = I+abs(J) and updates MXSPTR
(see Section 3.1). MPSTOV is called if MXR is too small.
0
MPNZR
+_____
0Normalizes and rounds a multiple-precision number (represented in a
nonstandard format, with separate sign and exponent, and possibly some
guard digits). Rounding depends on RNDRL in COMMON /MPCOM/ - see the
description of MPCQM. Called by MPADD, MPDIVI, MPMUL, MPMULI, etc., and
not recommended for independent use.
0
MPOUT
+_____
0CALL MPOUT (X, C, P, N)
0Converts multiple-precision X to FP.N format in C, which may be printed
under PA1 format. (Here FP.N means sign, P-N-2 'decimal' places before
the 'decimal' point, and N places after it.) Note that N = -1 is
allowed, and effectively gives IP format (i.e., sign and P-1 'decimal'
digits). Digits after the 'decimal' point are blanked out if they could
not be significant. The output base is OUTBAS in COMMON /MPCOM/ (see
Section 3.1), in the range two to sixteen, with default value ten. P
and N are integers, C is an integer array of dimension at least P. X
may be packed or unpacked.
1Sec. 6 MP User's Guide Page 56
+ _______________
0Rounding is defined by the parameter RNDRL in COMMON /MPCOM/, as
follows -
0 0 or 1 - Round to approximately the nearest decimal number
representable in the specified format. The error in the
conversion to decimal is less than 0.6 units in the last
decimal place.
0 2 - Round down (towards -infinity) to a decimal number
representable in the specified format.
0 3 - Round up (towards +infinity) to a decimal number
representable in the specified format.
0Efficiency is higher if B is a power of 10 than if not. Time = O(T**2).
0Note - MPFOUT and MP40D are more convenient to use, although less
flexible than MPOUT.
0
MPOUTE
+______
0CALL MPOUTE (X, C, J, P)
0Assumes X is a multiple-precision number and C an integer array of
dimension at least P > 3. J and P are integers. On return J is the
exponent (to base OUTBAS) of X and the fraction is in C, ready to be
printed in A1 format. For example, we could print J and C in format
'(I10, 1X, PA1)'. The fraction has one place before the 'decimal' point
and P-3 after. The output base is OUTBAS in COMMON /MPCOM/ (see Section
3.1), with default value ten. Rounding is not the best possible, but
the directed rounding options (RNDRL = 2 and 3) give correct lower and
upper bounds on the true result (see the description of MPOUT). X may
be packed or unpacked.
0Note - MPFOUT is a more convenient (though less flexible) output
routine.
0
MPOUTF
+______
0CALL MPOUTF (X, P, N, IFORM, ERR) or ERR = MPOUTF (X, P, N, IFORM)
or IF (MPOUTF (X, P, N, IFORM)) ...
0Writes multiple-precision number X on logical unit LUN (fourth word of
COMMON /MPCOM/) in format IFORM after converting to FP.N 'decimal'
representation using MPOUT. The output base is OUTBAS, with default
value ten. For further details see description of MPOUT. IFORM should
contain format which allows for output of P words in A1 format, plus any
desired spacing etc. ERR is returned as .TRUE. iff P <= 0 (see
description of MPIO). X may be packed or unpacked if MPFOUT is called
directly. For rounding options see MPOUT. Time = O(T**2).
0Note - a more convenient (though less flexible) output routine is
MPFOUT.
1Sec. 6 MP User's Guide Page 57
+ _______________
0MPOUT2
+______
0CALL MPOUT2 (X, C, P, N, NB)
0Same as MPOUT except that output representation is in base NB, where two
< NB < sixteen, e.g. NB = eight gives octal output, NB = sixteen gives
hexadecimal. Output digits are 0123456789ABCDEF. X is a
multiple-precision number, P, N and NB are integers, C is an integer
array of dimension at least P. Time = O(T**2). For rounding options
see the description of MPOUT.
0Note - MPOUT2 is superfluous now that OUTBAS determines the output base
for MPOUT. It is included for compatibility with earlier versions of
the MP package.
0
MPOVFL
+______
0CALL MPOVFL (X)
0Called on multiple-precision overflow, that is when the exponent of the
multiple-precision number X would exceed M. Execution is terminated
with an error message after calling MPMAXR(X) and setting MXEXPN = M+1
(see Section 3.1).
0
MPPACK
+______
0CALL MPPACK (X, Y) or Y = X or Y = CTP (X)
0Assumes that X is a multiple-precision number represented as usual in an
integer array of dimension at least T+2, and Y is an integer array of
dimension at least int((T+3)/2). X is stored in a compact format in Y,
and may be retrieved by calling MPUNPK (Y, X).
0MPPACK and MPUNPK are useful if space is critical, for example when
working with large arrays of multiple-precision numbers. X may be
packed or unpacked if MPPACK is called directly. (If X is packed the
effect of MPPACK is the same as that of MPSTR.)
0Augment users - X is type MULTIPLE, Y is type MULTIPAK.
0
The packed format is as follows -
0 word 2 = exponent (as for unpacked multiple-precision numbers),
0 words 1 and 3, ... , int((T+3)/2) = base B digits packed two per
word (digits pq are represented as p*B+q), with sign
attached to word 1.
0Thus, zero is represented in both packed and unpacked formats by a zero
first word, with the following words undefined. The first word has
absolute value greater than 1 if and only if the multiple-precision
number is nonzero and is represented in the packed format.
1Sec. 6 MP User's Guide Page 58
+ _______________
0MPPARA
+______
0I = MPPARA (A) or I = PARAM (A)
0Returns the integer value of the word of COMMON /MPCOM/ corresponding to
the Hollerith string A. For the allowable values of A, see the
description of MPPARM.
0
MPPARB
+______
0CALL MPPARB (I, A) or PARAM (A) = I
0Sets the word of COMMON /MPCOM/ corresponding to the Hollerith string A
to the integer value I. For the allowable values of A, see the
description of MPPARM. For the allowable values of I, see Section 3.1.
0
MPPARC
+______
0CALL MPPARC (N, J)
0Sets the N-th word of COMMON /MPCOM/ to the integer value J, for
0 < N < 24. Faster than MPPARB, but less mnemonic.
0
0MPPARM
+______
0CALL MPPARM (HOLL, SET, J)
0HOLL is a Hollerith string, SET is a Boolean variable, and J is an
integer. The first three characters of HOLL should agree with the first
three characters of one of the following -
0 (1) 'BASE' (2) 'NUMDIG' (3) 'MAXEXP' (4) 'LUN'
(5) 'MXR' (6) 'SPTR' (7) 'MXSPTR' (8) 'MNSPTR'
(9) 'MXEXPN' (10) 'MNEXPN' (11) 'RNDRL' (12) 'KTUNFL'
(13) 'MXUNFL' (14) 'DECPL' (15) 'MT2' (16) 'MXINT'
(17) 'EXWID' (18) 'INRECL' (19) 'INBASE' (20) 'OUTBAS'
(21) 'EXPCH' (22) 'CHWORD' (23) 'ONESCP'
0The corresponding word in COMMON /MPCOM/ will be set to the integer J if
SET = .TRUE., or its value will be returned in J if SET = .FALSE. For
the meanings and default settings of these parameters in COMMON /MPCOM/,
see Section 3.1.
0
0MPPARN
+______
0J = MPPARN (N)
0Returns the value of the N-th word of COMMON /MPCOM/, for 0 < N < 24,
zero otherwise. Faster than MPPARA, but less mnemonic.
1Sec. 6 MP User's Guide Page 59
+ _______________
0MPPI
+____
0
CALL MPPI (X) or X = CTM ('PI')
0Sets multiple-precision X = pi = 3.14159... to the available precision.
Uses the formula
0 pi = 16*arctan(1/5) - 4*arctan(1/239)
0of Machin. Time = O(T**2). (For asymptotically faster methods, see the
description of MPPIGL.) Rounding options are implemented as for MPEXP.
0
0MPPIGL
+______
0CALL MPPIGL (X)
0Sets multiple-precision X = pi = 3.14159... to almost the available
precision. Uses the Gauss-Legendre algorithm. This method requires
time O(ln(T)m(T)), so it is slower than MPPI if m(T) = O(T**2), but
would be faster for large T if a faster multiplication algorithm were
used (see the description of MPMUL). For a description of the
Gauss-Legendre algorithm see - R. P. Brent, 'Multiple-precision
zero-finding and the complexity of elementary function evaluation', in
Analytic Computational Complexity (edited by J.F. Traub), Academic
Press, l976, 151-176.
0Rounding options are not implemented, uses no guard digits. Recommended
for testing purposes only.
0
0MPPOLY
+______
0
CALL MPPOLY (X, Y, IC, N)
0Sets Y = IC(1) + IC(2)*X + ... + IC(N)*X**(N-1), where X and Y are
multiple-precision numbers, and IC is an integer array of dimension at
least N > 0. Rounding is not best possible, but the directed rounding
options (RNDRL = 2 and 3) give correct lower and upper bounds on the
true result.
0
0MPPWR
+_____
0
CALL MPPWR (X, N, Y) or Y = X**N
0Returns Y = X**N, for multiple-precision X and Y, integer N, with
0**0 = 1. X may be packed or unpacked, Y is unpacked. Rounding options
are implemented as for MPEXP.
1Sec. 6 MP User's Guide Page 60
+ _______________
0MPPWRA
+______
0CALL MPPWRA (X, N, Y)
0Returns Y = X**abs(N) for (unpacked) multiple-precision X and Y,
integer N, with 0**0 = 1. Uses truncated rather than rounded
arithmetic, and no guard digits, but directed rounding options give
correct upper and lower bounds. Called by MPEXP, MPPWR etc., not
recommended for independent use (use MPPWR instead).
0
MPPWR2
+______
0CALL MPPWR2 (X, Y, Z) or Z = X**Y
0Returns Z = X**Y for multiple-precision numbers X, Y and Z, where X is
positive (X = 0 is allowed if Y > 0). X and Y may be both packed or
both unpacked, Z is unpacked. Uses MPLN and MPEXP, so slower than MPPWR
and MPQPWR. Rounding options are implemented as for MPEXP.
0
MPQPWR
+______
0CALL MPQPWR (I, J, K, L, X) or X = QPWR (I, J, K, L)
0Sets multiple-precision X = (I/J)**(K/L) for integers I, J, K and L
(with the obvious conditions to ensure that a real result is defined).
Uses MPROOT if abs(L) is small, otherwise uses MPLNI and MPEXP.
Rounding options are implemented as for MPEXP.
0Augment users - X = QPWR (I, J, K, L) is usually faster than
X = CTM (I, J) ** CTM (K, L).
0
0MPREC
+_____
0CALL MPREC (X, Y) or Y = REC (X)
0Returns Y = 1/X, for multiple-precision X and Y. Uses MPDIVL if T is
small or RNDRL > 0, MPROOT otherwise. Time = O(m(T)) if RNDRL = 0,
O(T**2) if RNDRL > 0.
0Augment users - Y = REC (X) is faster than Y = 1/X.
0
0MPRESN
+______
0CALL MPRESN (N)
0Restores T, M and RNDRL which are assumed to have been saved on the
stack in blank COMMON by a call to MPSAVN (N). The integer N should be
as returned by MPSAVN. SPTR is restored to its value before the call to
MPSAVN (i.e. N). For further details, see the description of MPSAVN.
1Sec. 6 MP User's Guide Page 61
+ _______________
0MPRES2
+______
0CALL MPRES2 (N)
0Has the same effect as CALL MPRESN (N) except that SPTR is not
changed.
0
MPREVR
+______
0CALL MPREVR (I)
0If I > 0 and RNDRL = 2 or 3, RNDRL is set to 3 or 2 respectively (this
reverses the direction of rounding - see the description of MPCQM).
0MPRND
+_____
0CALL MPRND (X, TX, Y, TY, K)
0Moves X + S*K*abs(X)*B**(-TY) appropriately rounded to Y, where X is
an unpacked multiple-precision number with TX digits, Y is an unpacked
multiple-precision number with TY digits, K is an integer, and
0 S = 0 if RNDRL = 0 or 1,
-1 if RNDRL = 2 (rounding down),
+1 if RNDRL = 3 (rounding up).
0Note that RNDRL = 0 has the same effect as RNDRL = 1, and that the value
of T is irrelevant.
0MPROOT
+______
0CALL MPROOT (X, N, Y) or Y = ROOT (X, N)
0Returns Y = X**(1/N) for integer N such that a real result is defined,
(packed or unpacked) multiple-precision X and (unpacked)
multiple-precision Y. Uses Newton's method without divisions, so
time = O(m(t)). Rounding options are implemented as for MPEXP.
0Augment users - Y = ROOT (X, N) is faster than Y = X ** CTM (1, N)
(and Y = X**(1/N) is incorrect as 1/N is an
integer expression, hence usually zero).
0MPSAVN
+______
0CALL MPSAVN (N)
0Saves T, M and RNDRL on the stack in blank COMMON (they may be restored
by calling MPRESN). Returns N = old value of SPTR.
0 T is saved in word N of blank COMMON,
M is saved in word N+1, and
RNDRL is saved in word N+2.
0Note - N and SPTR must not have the same address.
1Sec. 6 MP User's Guide Page 62
+ _______________
0MPSCAL
+______
0CALL MPSCAL (X, N, J)
0Sets X = X*(N**J) for unpacked multiple-precision X, integer J, and
small positive integer N. Intended for use when 1 < N <= 100. Rounding
is not best possible, but the directed rounding options give correct
upper and lower bounds.
0
MPSET
+_____
0CALL MPSET (LUN, DECPL, MT2, MXR)
0Has the same effect as
0 CALL MPSET2 (LUN, DECPL, MT2, 1, MXR)
0(see the description of MPSET2). MPSET is redundant, but is included
for compatibility with earlier versions of the MP package.
0
MPSETR
+______
0CALL MPSETR (N)
0Sets the parameter RNDRL in COMMON /MPCOM/ to N. N should be an integer
with value 0, 1, 2 or 3 (see Section 3.1).
0
0MPSET2
+______
0CALL MPSET2 (LUN, DECPL, MT2, MNSPTR, MXR)
0MPSET2 may be called to initialize the parameters in COMMON /MPCOM/.
All five arguments are of type INTEGER. LUN is the Fortran logical unit
to be used for error messages. The base (B) and number of digits (T)
are set to give the equivalent of at least DECPL > 0 floating decimal
digits. MT2 should be the dimension of the arrays used for MP
variables, so an error occurs if the computed value of T exceeds MT2-2.
MXR is the size of blank COMMON declared in the calling program, and
MNSPTR is the first location of blank COMMON available for use by MP
routines. (They will use words MNSPTR, MNSPTR+1, ... , MXR of blank
COMMON (counting from 1) for working storage.)
0Let MXINT be the largest integer of the form 2**k - 1 representable as a
signed integer on the machine (for a more precise definition, see the
description of MPLARG). The computed B and T satisfy the conditions
0 B**(T-1) >= 10**(DECPL-1) and 8*B*B-1 <= MXINT. Approximately
minimal T and maximal B satisfying these conditions are chosen.
0MPSET2 also sets the following default values. For the meanings of the
various parameters, see Section 3.1.
1Sec. 6 MP User's Guide Page 63
+ _______________
0 M = int (MXINT/4), SPTR = MNSPTR, MXSPTR = MNSPTR, RNDRL = 0,
KTUNFL = 0, MXUNFL = 0, MNEXPN = M+1, MXEXPN = -M, MXINT = as
above, EXWID = 6, INRECL = 80, INBASE = 10, OUTBAS = 10,
EXPCH = 'E', CHWORD = value returned by MPUPW, ONESCP = 0 or
1 (see Section 3.1).
0If these are not all as desired, change after the call to MPSET2. For
the methods of doing this, see Section 3.2.
0
MPSIGA
+______
0I = MPSIGA (X) or I = SGN (X)
0Returns the sign of a packed or unpacked multiple-precision number X,
i.e.
0 +1 if X > 0,
0 if X = 0,
-1 if X < 0.
0
MPSIGB
+______
0CALL MPSIGB (I, X) or SGN (X) = I
0Sets the sign of a (packed or unpacked) multiple-precision number X to
the sign of I. (X is set to zero if I = 0.) The exponent and digits of
X are unchanged, but the result must be a valid (packed or unpacked)
multiple-precision number.
0
MPSIGN
+______
0CALL MPSIGN (X, Y, Z) or Z = SIGN (X, Y)
0Sets Z = abs(X)*sign(Y) for (both packed or both unpacked)
multiple-precision variables X and Y, unpacked multiple-precision Z,
where
0 sign (Y) = +1 if Y > 0,
= 0 if Y = 0,
= -1 if Y < 0.
0
MPSIN
+_____
0CALL MPSIN (X, Y) or Y = SIN (X)
0Returns Y = sin(X) for multiple-precision X and Y. X may be packed or
unpacked, Y is unpacked. Uses MPCIS, so time = O(sqrt(T)m(T)).
0Rounding options are implemented as for MPCIS. The absolute error in
the computed result is small, i.e. O(B**(-T)), but the relative error
may be large if X is close to a nonzero multiple of pi.
1Sec. 6 MP User's Guide Page 64
+ _______________
0MPSINH
+______
0CALL MPSINH (X, Y) or Y = SINH (X)
0Returns Y = sinh(X) for multiple-precision numbers X and Y, X not too
large. X may be packed or unpacked, Y is unpacked. The method is to
use MPEXP or MPEXP1, so time = O(sqrt(T)m(T)). Rounding options are not
yet implemented, uses no guard digits.
0
0MPSQRT
+______
0CALL MPSQRT (X, Y) or Y = SQRT (X)
0Returns Y = sqrt(X) for packed or unpacked multiple-precision X > 0,
unpacked multiple-precision Y. Uses MPROOT, so time = O(m(T)).
Rounding options are implemented as for MPEXP. (RNDRL = 0 or 1 gives
the exact result if it is exactly representable.)
0
0MPSTOV
+______
0CALL MPSTOV (N)
0MPSTOV is called if the working space in blank COMMON is too small. If
possible, the space should be expanded by at least N words, and MXR
should be increased by the number of words expanded. The new space must
be contiguous with the old (i.e. it must include words MXR+1, ... ,
MXR+N of blank COMMON).
0As distributed MPSTOV does nothing, because the method of expanding the
working space is machine and operating system dependent (see Section
9.3).
0
MPSTR
+_____
0CALL MPSTR (X, Y) or Y = X
0Sets Y = X for multiple-precision X and Y. X may be packed or unpacked.
Y will be packed if X is packed, unpacked if X is unpacked. (This is an
exception to the general rule that an unpacked result is returned even
if the argument(s) are packed.)
0
0MPSUB
+_____
0CALL MPSUB (X, Y, Z) or Z = X - Y
0Computes Z = X - Y for multiple-precision X, Y and Z. Rounding
options are implemented - see the description of MPADD with Y replaced
by -Y.
1Sec. 6 MP User's Guide Page 65
+ _______________
0MPSUBA
+______
0CALL MPSUBA (X, Y, I)
0Returns I = max (1, X(1) - Y(1)) for integer arrays X and Y. Called
by MPLARG and MPUPW with Hollerith arguments X and Y.
0MPTAN
+_____
0CALL MPTAN (X, Y) or Y = TAN (X)
0Sets Y = tan(X) for multiple-precision X and Y. X may be packed or
unpacked, Y is unpacked. Uses MPCIS, so time = O(sqrt(T)m(T)).
Rounding options are not implemented, but some guard digits are used.
0MPTANH
+______
0CALL MPTANH (X, Y) or Y = TANH (X)
0Returns Y = tanh(X) for multiple-precision X and Y. X may be packed or
unpacked, Y is unpacked. Time = O(sqrt(T)m(T)). Rounding options are
implemented as for MPEXP.
0MPTLB
+_____
0I = MPTLB (J)
0Returns an upper bound on abs(J)*ln(B)/ln(2), using only integer
arithmetic. Called by various MP routines.
0MPUNFL
+______
0CALL MPUNFL (X)
0Called on multiple-precision underflow, that is when the exponent of the
multiple-precision number X would be less than 1-M. The underflow
counter KTUNFL is incremented. Other action depends on the parameters
RNDRL and MXUNFL in COMMON /MPCOM/, as follows -
0 If RNDRL < 2 then
0 if MXUNFL = 0 (the default option), X is set to zero and execution
continues,
0 if MXUNFL > 0, X is set to zero and execution continues unless
MXUNFL underflows have occurred, when execution is
terminated by a call to MPERR.
0 If RNDRL = 2, action is as above except that X is set to the smallest
negative representable number, i.e. -(B**(-M)).
0 If RNDRL = 3, action is as above except that X is set to the smallest
positive representable number, i.e. B**(-M). (See
also under MPMINR.)
1Sec. 6 MP User's Guide Page 66
+ _______________
0MPUNFR
+______
0CALL MPUNFR (LUNIT, X) or X = UNFIO (LUNIT)
0Reads a multiple-precision number X from Fortran logical unit LUNIT,
assuming it has been written using MPUNFW (with the same B and T). X is
unpacked.
0MPUNFR and MPUNFW are faster than MPFIN and MPFOUT, so should be used if
temporary storage of multiple-precision numbers is required.
0MPUNFW
+______
0CALL MPUNFW (X, LUNIT) or UNFIO (LUNIT) = X
0Writes an multiple-precision number X without formatting on Fortran
logical unit LUNIT. X may be packed or unpacked if MPUNFW is called
directly, but must be unpacked if called via Augment. To save space on
LUNIT, X is packed before being written, so an unformatted record of
length int((T+3)/2) words is is actually written. Multiple-precision
numbers written with MPUNFW may be read with MPUNFR (so long as the
values of B and T have not been changed). See also the description of
MPUNFR above.
0MPUNPK
+______
0CALL MPUNPK (Y, X) or X = Y or X = CTM (Y)
0Unpacks the multiple-precision number which is stored in packed format
in Y, and stores the result in X. This reverses the effect of
0 CALL MPPACK (X, Y).
0If called directly, Y may be packed or unpacked, and in the latter case
the effect is the same as that of MPSTR. For further details see the
description of MPPACK.
0Augment interface users - X is type MULTIPLE,
Y is type MULTIPAK.
0MPUPDT
+______
0CALL MPUPDT (J)
0Sets MXEXPN = max (J, MXEXPN) and MNEXPN = min (J, MNEXPN).
0MPUPK
+_____
0CALL MPUPK (SOURCE, DEST, LDEST, LFIELD)
0This subroutine unpacks a packed Hollerith string SOURCE, placing one
character per word in the integer array DEST (as if read in A1 format).
It continues unpacking until it finds a sentinel ($) or until it has
filled LDEST words of the array DEST. The length of the unpacked string
is returned in LFIELD. Thus 0 <= LFIELD <= LDEST.
1Sec. 6 MP User's Guide Page 67
+ _______________
0MPUPW
+_____
0CALL MPUPW (W, C, N)
0W and N are INTEGER variables, C is an INTEGER array. When called with
a packed character string in W, MPUPW returns
0 N = number of characters per word
and
C(1), ... , C(N) = the characters in W (unpacked,
left justified, blank filled).
0Important note
+______________
0MPUPW is machine-dependent, and the version supplied with the MP package
does not work on all machines (e.g. machines which do not perform
integer arithmetic on full words, machines in which characters are
stored in an unusual order, and machines with decimal or
sign-and-magnitude arithmetic). For conversion hints, see Section 9.2.
0
0MPZETA
+______
0CALL MPZETA (N, X) or X = ZETA (N)
0Returns multiple-precision X = zeta(N) for integer N > 1, where zeta(N)
is the Riemann zeta function (the sum from j = 1 to infinity of
j**(-N)). Uses the Euler-Maclaurin series unless N = 2, 3, 4, 6 or 8.
In the worst case space = k*T/2 + O(T), where k is the number of terms
used in the Euler-Maclaurin series, k = O(T). Time = O(T**3) in
general, O(T**2) if N = 2, 3, 4, 6 or 8. Rounding options are
implemented as for MPEXP.
0
MP40D
+_____
0CALL MP40D (N, X)
0Output routine called by TEST program. Writes multiple-precision X to N
'decimal' places on unit LUN. The output base is OUTBAS (default value
ten - see Section 3.1). It is assumed that -OUTBAS < X < OUTBAS**2.
Rounding options are implemented as for MPOUT.
0
MP40F
+_____
0CALL MP40F (N, X)
0Equivalent to
0 CALL MPFOUT (X, N)
0MP40F is redundant, but is included for compatibility with earlier
versions of the MP package.
1Sec. 7.1 MP User's Guide Page 68
+ _______________
07. TEST PROGRAMS
+________________
0
07.1 TEST
+________
0
This main program computes the constants given in Appendix A of Knuth,
The Art of Computer Programming, Vol. 3. The constants are printed in
the same order as they are given in Knuth. Their correct values
(rounded to 40 decimal places) are also given in the comments in the
source of TEST.
0TEST computes the constants to 40 decimal places, but to increase the
accuracy it is only necessary to change the statement PLACES = 40, and
possibly the parameters of the call to MPSET2 and the dimensions of the
arrays, as described in the comments in the source of TEST. The
constants are given to 1000 decimal places in - R. P. Brent, 'Knuth's
constants to 1000 decimal and 1100 octal places', Tech. Report 47,
Computer Centre, Australian National University, Canberra (Sept. 1975).
0Execution time (CPU SUPs) on a Univac 1100/82 computer using FTN (level
9R1) is about 2.464 seconds.
0
07.2 TEST2
+_________
0
This main program tests various multiple-precision routines, especially
those not called by program TEST. It computes the constants given (in a
different order) in Computer Approximations (by Hart, Cheney, Lawson,
Maehly, Mesztenyi, Rice, Thacher and Witzgall), John Wiley, 1968,
Appendix C, pp. 182-183, and various other constants which are described
in the comments in the source of TEST2.
0Most of the constants are computed to 40 significant figures, with
working precision equivalent to at least 42 significant figures. To
increase the precision, it is only necessary to make a few changes
described in the comments. The correct output on a Univac 1100/82 is
given in the comments. The output may be slightly different on other
machines, especially if they have integer wordlength other than 36 bits.
0Execution time is about 25 times longer than that of TEST.
0
07.3 Other test programs
+_______________________
0
The EXAMPLE program (given in Section 2) and the JACOBI program (given
in Section 4.3) test various routines in the MP package, although not as
thoroughly as TEST and TEST2 do. The JACOBI program also provides a
partial test of the Augment precompiler and the MP description deck (see
Section 4).
1Sec. 8.1 MP User's Guide Page 69
+ _______________
08. COMPARISON WITH EARLIER VERSIONS OF MP
+_________________________________________
08.1 New routines and capabilities
+_________________________________
0For a brief description of the differences between MP version 781124 and
earlier versions, see Section 1.3. The main differences between version
800207 (or later) of the MP package and version 781124 are -
0 1. Introduction of the labelled common block COMMON /MPCOM/ (see
Section 3.1).
0 2. Improved dynamic storage allocation. A stack is now implemented,
with stack pointer SPTR, minimum and maximum pointers MNSPTR and
MXSPTR (see Section 3.1). The environment may be saved with
MPSAVN, working space allocated with MPNEW or MPNEW2, and the
environment later restored with MPRESN.
0 Earlier versions of the MP package did not have a genuine stack,
and it was necessary for the high-level routines to know how much
space the low-level routines required, which made it very difficult
to modify any low-level routines. Now it is much easier to add
high-level routines or to modify low-level routines, as it is not
necessary to know the space requirements of routines called by them
or calling them.
0 3. Although the stack is still in blank COMMON (for a reason given in
Section 1.2), it no longer must start at the first word of blank
COMMON. Thus, routines which use blank COMMON for other purposes
can call MP routines. The facility to start the stack anywhere in
blank COMMON is also useful on systems where working storage can be
obtained at runtime.
0 4. Rounding options have been implemented, both for the basic
arithmetic operations, and for most of the elementary and special
functions implemented in the MP package (see Section 6). Thus,
implementation of a multiple-precision interval-arithmetic package
should now be straightforward.
0 5. Many routines in the MP package now accept both packed and unpacked
arguments, and the Augment description deck has been modified
accordingly. To achieve this the packed representation had to be
changed (see below).
0 6. The MP routines are now independent of the Fortran routines ALOG,
SQRT etc., which makes MP useful for testing them.
0 7. The MP routines now depend very little on single-precision REAL
arithmetic. Only the conversion and comparison routines MPCDM,
MPCMD, MPCMDE, MPCMPD, MPCMPR, MPCMR, MPCMRE and MPCRM use REAL or
DOUBLE PRECISION arithmetic. This makes the MP package useful for
testing the correctness of REAL arithmetic, and practically
eliminates machine-dependence of the results obtained using the
package. It also makes implementation of most of the routines in
the MP package possible on small machines which do not have
floating-point (i.e. REAL) arithmetic.
1Sec. 8.2 MP User's Guide Page 70
+ _______________
0The following routines have been added to the MP package since the
release of version 781124. For descriptions of these routines, see
Section 6.
0MPATN2, MPCEIL, MPCHEB, MPCHEV, MPCHGB, MPCIS, MPCMP, MPCMPD, MPCMPQ,
MPCMUL, MPDIGS, MPDIGV, MPDIGW, MPDIM, MPDIVL, MPDIV2, MPDIV3, MPERRM,
MPFIN, MPFLOR, MPFOUT, MPGD, MPGD3, MPGET, MPIMUL, MPINTG, MPIS, MPKADD,
MPKDIV, MPKDVI, MPKIML, MPKMLI, MPKMUL, MPKSUB, MPK3V, MPK3V2, MPLARG,
MPLG10, MPMOD, MPMOVE, MPMULS, MPNEW, MPNEW2, MPPARA, MPPARB, MPPARC,
MPPARM, MPPARN, MPPWRA, MPRESN, MPRES2, MPREVR, MPRND, MPSAVN, MPSCAL,
MPSETR, MPSET2, MPSIGN, MPSTOV, MPSUBA, MPTLB, MPUNFR, MPUNFW, MPUPDT,
MPUPW.
0
8.2 Incompatibilities with earlier versions of MP
+_________________________________________________
0The parameters B, T, M, LUN and MXR, which were formerly the first five
words of blank COMMON, are now the first five words of the labelled
common block COMMON /MPCOM/ (see Section 3.1). Thus, any routines which
manipulate these parameters directly will have to declare them in
COMMON /MPCOM/.
0The format for packed multiple-precision numbers has been changed - the
exponent is now stored in the second word, not the first. The reason
for this change is that with the new format it is possible to tell, by
inspecting the first word of a multiple-precision number, whether it is
packed or unpacked (unless it is zero, when only the first word is
defined anyway). For more details, see the description of MPPACK in
Section 6.
0Some routines used to return an INTEGER error flag. The type of all
error flags is now LOGICAL. In all cases the integer value 0 has been
replaced by the logical value .FALSE., and the integer value 1 has been
replaced by the logical value .TRUE.. Routines affected are MPERF3,
MPHANK, MPIN, and MPINE.
0The calling sequence of MPADD3 has been changed. This routine should
not be called directly by users of the MP package.
0The routines MPINE, MPOUT2, MPSET and MP40F are now redundant, and may
be deleted from later versions of the MP package. It is recommended
that MPIN, MPOUT (with the correct setting of OUTBAS), MPSET2 and MPFOUT
(respectively) should be used instead.
0The following routines have been removed from the MP package since the
release of version 781124 -
0MPADD2, MPCLR, MPEXT, MPKSTR, MPMUL2, MPSIN1, MP40E, MP40G, TESTV and
TIMEMP.
0MPKSTR is no longer needed, as MPSTR now works for both packed and
unpacked arguments. The test program TESTV may be reconstructed by
making trivial changes to the program TEST. The other routines are
low-level routines whose disappearance should not affect a user of the
MP package.
1Sec. 9.1 MP User's Guide Page 71
+ _______________
09. INSTALLATION INSTRUCTIONS
+____________________________
09.1 Essentials
+______________
0MP is normally distributed on an unlabelled, 9 track, Ebcdic or Ascii,
800 or 1600 fpi magnetic tape, with logical record length 80 characters,
blocking factor 12, and the following 6 files -
0 File 1 - Comments and EXAMPLE program.
File 2 - MP subroutines.
File 3 - Test programs TEST and TEST2.
File 4 - This User's Guide (duo-case version).
File 5 - Augment description deck and JACOBI program using it.
File 6 - This User's Guide (upper-case version).
0To install MP, read these six files. Print file 4 or 6 (the User's
Guide) using the first character as standard Fortran printer control.
0Check the source of routines MPINIT, MPIS, MPLARG, and MPUPW (in file
2), make any necessary changes, then compile the routines in files 1, 2
and 3 (preferably with compiler optimization options turned off).
0Convert compiled routines from file 2 into a relocatable library.
0Execute the programs EXAMPLE, TEST and TEST2 from files 1 and 3 (after
linking to required routines from file 2) and check that output is
correct (see the comments in the source of the programs). Output is on
unit 6 unless the first argument in each call to MPSET2 is changed.
0If all has gone well, try recompiling file 2 with compiler optimization
options turned on, and rerun the test programs. Any problems which
appear are probably due to bugs in your compiler, not in the MP
routines. (Such problems exist with some releases of the Univac FTN and
DEC 10 (F10) compilers, for example.) The routines whose optimization
is most worthwhile are MPADD3, MPDIVI, MPDIV2, MPDIV3, MPMLP and MPNZR.
0The Augment description deck is supplied with the MP package, but the
Augment precompiler is not. If you want to use the Augment interface,
obtain Augment from the Programming Services Supervisor, Mathematics
Research Center, 610 Walnut Street, Madison, Wisconsin 53706, and get it
running. This should not be too hard if you have a Univac 1100, IBM
360/370, CDC 6000/7000, DEC 10, or Honeywell 600/6000, as Augment has
already been implemented on these machines. (Augment is written mainly
in portable Fortran, but there are a few machine-dependent routines.)
0Next, use the description deck supplied in file 5 and run the JACOBI
test program which follows the description deck in file 5 (after
inserting a few machine-dependent control cards). It will be necessary
to change the dimension statements in the description deck and modify
routine MPINIT if your machine has wordlength less than 16 bits, and
desirable to do likewise if the wordlength is greater than 16 bits
(because arrays used for multiple-precision variables will be declared
larger than necessary). See the comments in MPINIT regarding
declaration of COMMON if your system insists that all common blocks be
declared in the main program.
1Sec. 9.2 MP User's Guide Page 72
+ _______________
09.2 Conversion notes
+____________________
0
To convert MPINIT and the description deck, let
0 i = dimension of arrays for MP variables (see the description deck
and MPINIT),
0 j = dimension of arrays for packed MP variables (see the
description deck),
0 k = size of work area in blank COMMON
= MXR + 1 - MNSPTR (see Section 3.1).
0
Suppose precision equivalent to at least d decimal places is required on
a machine with effective wordlength w bits (MXINT = 2**(w-1)-1 must be
representable as a signed integer - see also the description of MPLARG
in Section 6).
0Then
i = T + 2,
0 j = int ((T + 3)/2),
and
k <= T*T + 15*T + 200 (usually much less),
where
T = ceiling (1 + (d-1)*log(10)/log(B))
0is the number of base B digits to be used, B = 2**int(w/2 - 2).
0
0To convert MPIS on a Burroughs 6000/7000 machine, change 'EQ' to 'IS' in
the first executable statement of MPIS. On some other machines it may
be necessary to mask off some characters or shift right before the
comparison.
0
To convert MPLARG on a Cyber 76, Burrough 6000/7000, or PDP 11 machine,
set the appropriate logical variable in the DATA statements of MPLARG to
.TRUE.. On other machines it may be necessary to set WDLEN to the
effective wordlength (for integer arithmetic) in MPLARG.
0
To convert MPUPW on most machines, replace the body of MPUPW by
0 INTEGER W, C(chword)
N = chword
DECODE (N, 10, W) C
10 FORMAT (80A1)
0where 'chword' is the number of characters per word. On some CDC
machines the arguments of DECODE must be given in a different order.
(If chword > 10, trivial changes are required in MPSET2 and MPUPW.)
1Sec. 9.3 MP User's Guide Page 73
+ _______________
09.3 Desirable changes
+_____________________
0The following machine-dependent changes are desirable but not essential,
and may not be possible on some systems.
0MPIO could be modified to return ERR = .TRUE. on a read/write error.
0MPUPK could be modified to replace a compiler-generated end-of-string
sentinel (if any) by '$', so that string arguments of MPCAM need not be
terminated by '$'. See the DATA statement for SENTNL in MPUPK.
0MPSTOV could be modified to expand the working storage when necessary.
0MPERR could be modified to give a trace-back when an error is detected.
0MPDIGV could be made much more efficient. This would be worthwhile if
MPIN or MPFIN were to be used extensively.
0
09.4 Note on Fortran 77
+______________________
0The MP routines do not conform strictly to the new 'Fortran 77' Standard
(ANS X3.9-1978). There are two main violations of the standard. One is
that arrays which are actual parameters of subprograms, and are declared
with dimension (1), should instead be declared with dimension (*), e.g.
0 SUBROUTINE MPABC (X, Y)
INTEGER X(1), Y(1)
...
0should be replaced by
0 SUBROUTINE MPABC (X, Y)
INTEGER X(*), Y(*)
...
0The other violation is the use of the Hollerith data type, which is not
included in the Fortran 77 Standard.
0Most Fortran 77 compilers accept these violations as 'extensions' to the
new standard, in order to be compatible with the old Fortran Standard
(ANS X3.9-1966). This is true, for example, of the Univac 1100 FTN
compiler (level 9R1). Thus, we do not anticipate any major problems
with Fortran 77 compilers.
1End of MP User's Guide.