Richard P. Brent. 1978. Algorithm 524: MP, A Fortran Multiple-Precision Arithmetic Package [A1]. ACM Trans. Math. Softw. 4, 1 (March 1978), 71–81. https://doi.org/10.1145/355769.355776
- for: multiple precision floating point arithmetic and evaluating elementary and special functions
- calgo note: You almost surely want the newer version (June 1981) in netlib/bmp.
In general, we recommend the use of a more modern package, for example:
- Paul Zimmermann's MPFR package (written in C using GMP);
- David Bailey's MPFUN90 package (written in Fortran 90) or ARPREC package (written in C++/Fortran 90).
files:
524.gz524.f77: all in one file
524.zipDoc/Drivers/Src/
See: Doc/userg.as
0 MP USER'S GUIDE
**************
0
0 RICHARD P. BRENT,
COMPUTER CENTRE,
AUSTRALIAN NATIONAL UNIVERSITY,
BOX 4, CANBERRA, ACT 2600,
AUSTRALIA
0
0 TECHNICAL REPORT NO. 54,
SEPTEMBER 1976 (LAST REVISED AUGUST 1978)
0 CONTENTS PAGE
******** ****
0 1. GENERAL DESCRIPTION OF MP 1.1
RESTRICTIONS AND EFFICIENCY CONSIDERATIONS 1.2
HISTORY AND REFERENCES 1.3
0 2. SUMMARY OF MP ROUTINES 2.1
0 3. EXAMPLE PROGRAM 3.1
0 4. THE AUGMENT INTERFACE 4.1
THE DESCRIPTION DECK 4.1
EXAMPLE PROGRAM (JACOBI) USING AUGMENT 4.3
RESERVED WORDS 4.5
0 5. INSTALLATION INSTRUCTIONS 5.1
CONVERSION NOTES 5.2
0 6. DESCRIPTION OF MP ROUTINES AND TEST PROGRAMS 6.1
MP SUBROUTINES 6.1
MP TEST PROGRAMS 6.23
0 7. INDEX OF LINE NUMBERS 7.1
1 MP USERS GUIDE (AUGUST 1978) PAGE 1.1
01 GENERAL DESCRIPTION OF MP
***************************
0 MP 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 COMPILER, SUFFICIENT MEMORY, AND A WORDLENGTH OF AT LEAST
16 BITS. SOME MODIFICATIONS WOULD BE NECESSARY FOR A WORDLENGTH OF LESS
THAN 16 BITS.
0 MP HAS BEEN TESTED ON VARIOUS MACHINES INCLUDING A UNIVAC 1108 (E LEVEL
FORTRAN V), A UNIVAC 1100/42 (E AND T LEVEL FORTRAN V, ASCII FORTRAN, AND
RALPH), A DEC 10 (FORTRAN 10 (/NOOPT), AND FORTRAN 40), AN IBM 360/50
(FORTRAN G AND FORTRAN H, OPT = 2), AN IBM 360/91 AND 370/168 (FORTRAN H
EXTENDED, OPT = 2), A CYBER 76 (FTN 4.2, OPT = 1) AND A PDP 11/45 (DOS).
THESE MACHINES HAVE EFFECTIVE INTEGER WORDLENGTHS RANGING FROM 16 TO 48
BITS.
0 MP 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.
0 T-DIGIT FLOATING-POINT NUMBERS ARE STORED IN INTEGER ARRAYS OF DIMENSION
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).
0 NOTE THAT WORDS 2 TO T+2 ARE UNDEFINED IF SIGN = 0.
0 ARITHMETIC IS ROUNDED, AND FOUR GUARD DIGITS ARE USED FOR ADDITION AND
MULTIPLICATION, SO THE CORRECTLY ROUNDED RESULT IS USUALLY PRODUCED.
DIVISION, SQRT ETC ARE DONE BY NEWTONS METHOD, BUT GIVE THE EXACT RESULT IF
IT CAN BE REPRESENTED WITH T-2 DIGITS. OTHER ROUTINES (MPSIN, MPLN ETC)
USUALLY GIVE A RESULT Y = F(X) WHICH COULD BE OBTAINED BY MAKING AN
O(B**(1-T)) PERTURBATION IN X, EVALUATING F EXACTLY, THEN MAKING AN
O(B**(1-T)) PERTURBATION IN Y.
0 EXPONENTS CAN LIE IN THE RANGE -M, ... , +M INCLUSIVE, WHERE M IS SET BY
THE USER. ON UNDERFLOW DURING AN ARITHMETIC OPERATION, THE RESULT IS SET
TO ZERO BY SUBROUTINE MPUNFL. ON OVERFLOW SUBROUTINE MPOVFL IS CALLED AND
EXECUTION IS TERMINATED WITH AN ERROR MESSAGE.
0 ERROR 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. IT
IS ASSUMED THAT LOGICAL RECORDS OF UP TO 80 CHARACTERS MAY BE WRITTEN ON
UNIT LUN. A WORKING ARRAY OF SIZE MXR (SEE BELOW) MUST BE PROVIDED IN
COMMON.
1 MP USERS GUIDE (AUGUST 1978) PAGE 1.2
0 THE PARAMETERS B, T, M, LUN AND MXR ARE PASSED TO THE UTILITY ROUTINES IN
COMMON, TOGETHER WITH A WORKING ARRAY R WHICH MUST BE SUFFICIENTLY LARGE
(SEE BELOW). MOST ROUTINES USE THE STATEMENTS
0 COMMON B, T, M, LUN, MXR, R
INTEGER B, T, R(1)
0 AND IT IS ASSUMED THAT MXR IS SET TO THE DIMENSION OF R IN THE CALLING
PROGRAM, AND THAT MXR IS SUFFICIENTLY LARGE (SEE SECTION 6).
0 RESTRICTIONS AND EFFICIENCY CONSIDERATIONS
******************************************
0 B (THE BASE) MUST BE AT LEAST 2.
0 T (NUMBER OF DIGITS) MUST BE AT LEAST 2.
0 M (EXPONENT RANGE) MUST BE GREATER THAN T AND LESS THAN 1/4 THE LARGEST
MACHINE-REPRESENTABLE INTEGER.
0 8*B**2-1 MUST BE NO GREATER THAN THE LARGEST MACHINE-REPRESENTABLE
INTEGER, AND THE INTEGERS 0, 1, ... , B MUST BE EXACTLY REPRESENTABLE AS
SINGLE-PRECISION FLOATING-POINT NUMBERS.
0 B**(T-1) SHOULD BE AT LEAST 10**7.
0 B AND T MAY BE SET TO GIVE THE EQUIVALENT OF A SPECIFIED NUMBER OF DECIMAL
PLACES BY CALLING MPSET (SEE SECTION 6), OR MAY BE SET DIRECTLY BY THE
USER. IF MPSET IS NOT CALLED, THE USER MUST REMEMBER TO INITIALIZE M, LUN
AMD MXR (SEE ABOVE) AS WELL AS B AND T BEFORE CALLING ANY MP ROUTINES.
0 IT WOULD BE POSSIBLE TO USE LABELLED COMMON INSTEAD OF BLANK COMMON
THROUGHOUT, AND SET DEFAULT INITIALIZATIONS IN A DATA STATEMENT. BLANK
COMMON RATHER THAN LABELLED COMMON IS USED FOR WORKING STORAGE ONLY 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.
0 FOR EFFICIENCY CHOOSE B FAIRLY LARGE, SUBJECT TO THE RESTRICTIONS GIVEN
ABOVE. FOR EXAMPLE, IF THE WORDLENGTH IS
0 48 BITS, COULD USE B = 4194304 OR 1000000,
36 BITS, COULD USE B = 65536 OR 10000,
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.
0 AVOID MULTIPLICATION AND DIVISION BY MP NUMBERS, AS THESE TAKE O(T**2)
OPERATIONS, WHEREAS MULTIPLICATION AND DIVISION BY (SINGLE-PRECISION)
INTEGERS TAKE O(T) OPERATIONS. SEE THE COMMENTS ON MPDIV, MPDIVI, MPMUL
AND MPMULI IN SECTION 6.
1 MP USERS GUIDE (AUGUST 1978) PAGE 1.3
0 MP NUMBERS USED AS ARGUMENTS OF SUBROUTINES NEED NOT BE DISTINCT. FOR
EXAMPLE,
0 CALL MPADD (X, Y, Y)
AND
CALL MPEXP (X, X)
0 ARE ALLOWABLE. HOWEVER, DISTINCT ARRAYS WHICH OVERLAP SHOULD NOT BE USED.
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 OR FOR ARRAYS IN COMMON).
APART FROM THESE VIOLATIONS, MP IS WRITTEN IN ANSI STANDARD FORTRAN (ANSI
X3.9-1966). THIS HAS BEEN CHECKED BY THE PFORT VERIFIER. THE ONLY
MACHINE-DEPENDENT ROUTINE IS MPUPK (WHICH UNPACKS CHARACTERS STORED SEVERAL
TO A WORD). OTHER ROUTINES WHICH MAY REQUIRE TRIVIAL CHANGES ARE MPSET,
MPINIT AND TIMEMP (SEE COMMENTS IN SECTIONS 4 AND 6 BELOW).
0 HISTORY AND REFERENCES
**********************
0 THE FIRST WORKING VERSION OF MP (VERSION 731101) WAS WRITTEN BY R. P. BRENT
IN NOVEMBER 1973. BETWEEN 1973 AND 1978 A CONSIDERABLE NUMBER OF
IMPROVEMENTS AND ADDITIONS WERE MADE. THE LAST MAJOR REVISION WAS IN APRIL
1978, WHEN THE AUGMENT INTERFACE ROUTINES (SEE SECTION 4 BELOW) WERE ADDED.
0 FOR AN INTRODUCTION TO THE DESIGN AND PHILOSOPHY OF MP, SEE - A FORTRAN
MULTIPLE-PRECISION ARITHMETIC PACKAGE (BY R. P. BRENT), 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 SECTION 6 BELOW.
0 A PREPROCESSOR (AUGMENT) WHICH FACILITATES THE USE OF THE MP PACKAGE IS
AVAILABLE. SEE - AN AUGMENT INTERFACE FOR BRENTS MULTIPLE PRECISION
ARITHMETIC PACKAGE (BY R. P. BRENT, J. A. HOOPER AND J. M. YOHE),
MATHEMATICS RESEARCH CENTER, UNIVERSITY OF WISCONSIN, MADISON, AUGUST 1978
(SUBMITTED TO ACM TRANS. MATH. SOFTWARE), AND SECTION 4 BELOW.
0 CORRESPONDENCE CONCERNING MP SHOULD BE ADDRESSED TO DR. R. P. BRENT,
COMPUTER CENTRE, AUSTRALIAN NATIONAL UNIVERSITY, BOX 4, CANBERRA, A.C.T.
2600, AUSTRALIA.
1 MP USERS GUIDE (AUGUST 1978) PAGE 2.1
02 SUMMARY OF MP ROUTINES
************************
0 BASIC ARITHMETIC
0 MPADD, MPADDI, MPADDQ, MPDIV, MPDIVI, MPMUL, MPMULI, MPMULQ, MPREC,
MPSUB
0 POWERS AND ROOTS
0 MPPWR, MPPWR2, MPQPWR, MPROOT, MPSQRT
0 ELEMENTARY FUNCTIONS
0 MPASIN, MPATAN, MPCOS, MPCOSH, MPEXP, MPLN, MPLNGS, MPLNI, MPSIN,
MPSINH, MPTAN, MPTANH
0 SPECIAL FUNCTIONS
0 MPBESJ, MPDAW, MPEI, MPERF, MPERFC, MPGAM, MPGAM, MPGAMQ, MPLI,
MPLNGM
0 CONSTANTS
0 MPBERN, MPEPS, MPEUL, MPMAXR, MPMINR, MPPI, MPPIGL, MPZETA
0 INPUT AND OUTPUT
0 MPDUMP, MPIN, MPINE, MPINF, MPOUT, MPOUTE, MPOUTF, MPOUT2
0 CONVERSION
0 MPCAM, MPCDM, MPCIM, MPCMD, CPCMDE, MPCMEF, MPCMI, MPCMIM, MPCMR,
MPCMRE, MPCQM, MPCRM
0 COMPARISON
0 MPCMPA, MPCMPI, MPCMPR, MPCOMP, MPEQ, MPGE, MPGT, MPLE, MPLT, MPNE
0 GENERAL UTILITY ROUTINES
0 MPABS, MPCLR, MPCMF, MPGCDA, MPGCDB, MPINIT, MPKSTR, MPMAX, MPMIN,
MPNEG, MPPACK, MPPOLY, MPSET, MPSTR, MPUNPK
0 ERROR DETECTION AND HANDLING
0 MPCHK, MPERR, MPOVFL, MPUNFL
0 TEST PROGRAMS
0 EXAMPLE, JACOBI, TEST, TESTV, TEST2
1 MP USERS GUIDE (AUGUST 1978) PAGE 2.2
0 AUGMENT INTERFACE ROUTINES
0 MPBASA, MPBASB, MPDGA, MPDGB, MPDIGA, MPDIGB, MPEXPA, MPEXPB, MPMEXA,
MPMEXB, MPSIGA, MPSIGB
0 MISCELLANEOUS ROUTINES USED BY THE ABOVE
0 MPADD2, MPADD3, MPART1, MPBES2, MPERF2, MPERF3, MPEXP1, MPEXT, MPGCD,
MPHANK, MPIO, MPLNS, MPL235, MPMLP, MPMUL2, MPNZR, MPSIN1, MPUPK,
MP40D, MP40E, MP40F, MP40G, TIMEMP
See: copyright.htm
Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page.
Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted.
To copy otherwise, to republish, to post on servers, or to redistribute to lists, requires prior specific permission and/or a fee.
Request permissions from Publications Dept, ACM Inc., fax +1 (212) 869-0481, or permissions@acm.org.
Collected Algorithms, Special Edition CD-ROM
Copyright © 2000 by the Association for Computing Machinery, Inc. 1-58113-338-5/00/11 ... $5.00