## A Fortran multiple-precision arithmetic package (1978)

Venue: | ACM Trans. Math. Softw |

Citations: | 60 - 3 self |

### BibTeX

@ARTICLE{Brent78afortran,

author = {Richard P. Brent},

title = {A Fortran multiple-precision arithmetic package},

journal = {ACM Trans. Math. Softw},

year = {1978},

volume = {4},

pages = {57--70}

}

### Years of Citing Articles

### OpenURL

### Abstract

A collection of ANSI Standard Fortran subroutines for performing multiple-precision floating-point arithmetic and evaluating elementary and special functions is described. The subrou-tines are machine independent and the precision is arbitrary, subject to storage limitations. The design of the package is discussed, some of the algomthms are described, and test results are given. Key Words and Phrases ' arithmetic, multiple precision, extended precision, floating point, elementary function evaluation, Euler's constant, gamma function, polyalgorithm, software package, Fortran, machine-independent software, special function evaluation, Bessel func-

### Citations

183 |
Multiplication of multidigit numbers on automata
- Karatsuba, Ofman
- 1963
(Show Context)
Citation Context ...uth [32]. We use R*-rounding [15, 33] after postnormalization, with four guard digits. 6.2 Multiplication The classical O(t ~) method is used because it seems difficult to implement faster algorithms =-=[29, 32, 46]-=- in a machine-independent manner in Fortran. Also, the faster algorithms would actually be slower for small t. The subroutine MPMUL could be modified without affecting the other routines in MP. We den... |

89 | Fast multiple-precision evaluation of elementary functions
- Brent
- 1976
(Show Context)
Citation Context ...fficiently small, and otherwise to use the relation exp(x) - 1 = [exp(x/2) - 1] [exp(x/2) + 1] to reduce Ix I. The time required is O(tl/2M(t)). -l'Iethods which require time O(log(t)M(t) ) are known =-=[11, 14]-=-, but turn out to be slower for small and moderate t. 5IPEXP uses the identity exp(x) = e ~ exp(f), where n is the integer part of x, and f = x - ~, to evaluate exp(x) (using 5~IPEXP1 to evaluate exp(... |

62 | Multiple-precision zero-finding methods and the complexity of elementary function evaluation
- Brent
- 1975
(Show Context)
Citation Context ...ivisions except by n. The iteration is n y~+l = y~ q- yj(1 -- xy~ )/~. .lIPROOT implements this iteration, starting with a small value of t and approximately doubling it at each step, as described in =-=[14, 17]-=-. Thus the time required by ~[PROOT is O(M(t)). MPREC implements the special case n = 1. Division and square roots require one additional multiplication after the computation of --1 X--I/2, x or respe... |

23 |
Schnelle Multiplikation grosser Zahlen. Computing 7
- Schönhage, Strassen
- 1971
(Show Context)
Citation Context ...uth [32]. We use R*-rounding [15, 33] after postnormalization, with four guard digits. 6.2 Multiplication The classical O(t ~) method is used because it seems difficult to implement faster algorithms =-=[29, 32, 46]-=- in a machine-independent manner in Fortran. Also, the faster algorithms would actually be slower for small t. The subroutine MPMUL could be modified without affecting the other routines in MP. We den... |

8 | On the precision attainable with various floating-point number systems
- Brent
- 1973
(Show Context)
Citation Context ...1 10 (1°) PDP 10 (7) F10 (9) 36 11.0 10 (1°) Univac 1108 Fortran V (11) 36 4.42 17 Univac 1100/42 Fortran V (12) 36 3.39 12 (13) Univac 1100/42 Ascii Fortran ('4) 36 3.28 13 (13) Cyber 76 Fortran 4.2 =-=(15)-=- 48 (16) 0.43 6 (1) No floating-point hardware, 28K words memory (3) Fortran V004A,/ER,/ON,/SU under DOS. (3) Fortran IV G (level 19) under OS/]VIFT. (4) Fortran H (level 20.1), opt = 2, under OS/MFT.... |

7 | A System for Polynomial Manipulation - Collins, PM - 1966 |

7 |
Euler's constant to 1271 places
- Knuth
(Show Context)
Citation Context ...72... is usually defined by "r= ,+~lim(~,-, 1/i- ln(n)). We may use this definition to evaluate % estimating the remainder after a finite n by the asymptotic Euler-Maclaurin series; see, for example, =-=[31]-=-. However, this method requires the Bernoulli numbers which appear as coefficients in the EulerMaclaurin series, and the number of Bernoulli numbers required increases with t if we demand reasonable e... |

6 |
On the computation of Euler’s constant
- Sweeney
(Show Context)
Citation Context ...hod here because similar ideas may be used for other multiple-precision calculations, for example, in the computation of F(x) (see Section 6.10). The method was suggested (though not used) by Sweeney =-=[48]-=- and depends on the identity .y = S(n) - R(n) - ln(n), where nk(_l) k-1 S(n) = 25 k-1 kt k ' R(n) = f/exp(--u) du exp(-n) 2 kt u ~ n +-o (-n)+' ACI~I Transactions on Mathematical Software, Vol 4, No 1... |

6 |
A Portable Extended Precision Arithmetic Package and Library with FORTRAN Precompiler
- Wyatt, Lozier, et al.
- 1976
(Show Context)
Citation Context ...or checking purposes, 1000D results are available [12]. Table II shows how the execution time depends on the precision required. 2 2 Some of the results were checked using the package of Wyatt et al. =-=[50]-=-. MP was slgmficantly faster: the ratio of execution times ranged from 1.49 (for 20-place accuracy) to 10.1 (for 200place accuracy). ACM Transactions on Mathematical Software, Vol 4, No 1, March 1978s... |

4 |
A precise numerical analysis program
- Aberth
- 1974
(Show Context)
Citation Context ... transmitted to the MP routines in COMMON. 3. CAPABILITIES OF MP The present version of SiP contains 101 subroutines and four main programs. The capabilities of the subroutines include the following: =-=(1)-=- conversion of integer, real, and double-precision numbers to multipleprecision format, and vice versa; A precompiler using an extension of MORTRAN2 has been written and is currently being tested. ACM... |

3 |
Error analysis of a computation of Euler's constant
- Beyer, Waterman
- 1974
(Show Context)
Citation Context ... using MPPIGL as well as published values of ~r, and MPEUL was checked using published values of ~ [31, 48]. The testing of MPEUL actually revealed an error in the most accurate published value of ~, =-=[5, 6]-=-. 5/IPGA54Q was checked using published values [21] as well as various identities. MPGASI and MPLNGM were tested using MPGAMQ and various identities. MPEI, 5,IPERF, 5IPERFC, and MPBESJ were tested for... |

2 |
Language facilities for multiple-precision floating-point computation
- HULL, HOrBAUER
- 1974
(Show Context)
Citation Context ...ations on doublelength integers. Arithmetic operations on mp numbers are performed by subroutine calls. Thus, instead of Z = X ~ Y we need to write CALL MPADD (X, Y, Z). A precompiler in the style of =-=[27]-=- or [52] could generate the appropriate subroutine calls. 1 Sufficient working space for the MP routines must be declared in COMMON in the main program. The parameters b, t, etc., are also transmitted... |

2 |
The Elements of Programming Style
- PLAUGER
- 1974
(Show Context)
Citation Context ... programming error if the allowable exponent range is large. To make the SiP routines as intelligible and easy to modify as possible, we have followed most of the suggestions of Kernighan and Plauger =-=[30]-=-. 6. SOME ALGORITHMS USED IN MP 6.1 Addition This is straightforward, and is performed much as described in Knuth [32]. We use R*-rounding [15, 33] after postnormalization, with four guard digits. 6.2... |

1 |
Handbook of Mathematwal Functions
- STEGUN
- 1964
(Show Context)
Citation Context ...sing an extension of MORTRAN2 has been written and is currently being tested. ACM Transactions on Mathematical Software, Vol. 4, No. 1, March 1978sA Fortran Multiple-Precision Arithmetic Package • 59 =-=(2)-=- multiplication and division of mp numbers by small integers; (3) addition, subtraction, multiplication, and division of mp numbers. (4) powers and roots of mp numbers, (5) elementary functions of mp ... |

1 |
American National Standard Fortran (ANSI X3.9-1966
- Nat
- 1966
(Show Context)
Citation Context ...rary-precision subroutines for all the commonly used special functions. We invite others to contribute such subroutines. 4. DEPARTURES FROM ANSI FORTRAN We have attempted to use ANSI Standard Fortran =-=[3]-=- as far as possible. The only known violation of the standard is that arrays used as arguments of MP subroutines (or in COMMON) are declared with dimension 1 in the subroutines, e.g. SUBROUTINE A (X, ... |

1 |
Multiple-precision arithmetic and the exact calculation of the 3-j, 6-j, and 9-3 symbols
- RDLICH
- 1964
(Show Context)
Citation Context ...March 1978sA Fortran Multiple-Precision Arithmetic Package • 59 (2) multiplication and division of mp numbers by small integers; (3) addition, subtraction, multiplication, and division of mp numbers. =-=(4)-=- powers and roots of mp numbers, (5) elementary functions of mp numbers (log, exp, sin, tan, arcsin, arctan, sinh, eosh, tanh), (6) some special functions and constants (Bernoulli numbers, Bessel func... |

1 |
Decimals and partial quotients of Euler's constant and ln(2
- BEYER, WATERMAN
(Show Context)
Citation Context ... using MPPIGL as well as published values of ~r, and MPEUL was checked using published values of ~ [31, 48]. The testing of MPEUL actually revealed an error in the most accurate published value of ~, =-=[5, 6]-=-. 5/IPGA54Q was checked using published values [21] as well as various identities. MPGASI and MPLNGM were tested using MPGAMQ and various identities. MPEI, 5,IPERF, 5IPERFC, and MPBESJ were tested for... |

1 |
An extended arithmetic package
- BLVM
- 1965
(Show Context)
Citation Context ...onstants (Bernoulli numbers, Bessel functions of the first kind, error and complementary error functions, exponential and logarithmic integrals, Dawson's integral, gamma function, ~r, % ~'(~), etc.); =-=(7)-=- fixed and floating-point decimal output and free-field decimal input of mp numbers; (8) integer and fractional parts of mp numbers; (9) routines for error handling, testing, and debugging; (10) misce... |

1 |
Algorithm 524. MP, A Fortran multzple-precision arithmetic package
- BRENT
- 1978
(Show Context)
Citation Context ...[31]. MP has also been used successfully to verify the orders of certain root-finding methods and to determine their asymptotic error constants [16]. 8. AN EXAMPLE The main program, EXAMPLE, given in =-=[9, 13]-=-, is intended to give a simple example of the use of MP. It computes 7r and exp(r(163)I/2/3) and prints them to 100 decimM places. The correct output is and 7r = 3.14159265358979... 1170680 exp0r(163)... |

1 |
Computation of the continued fraction for Euler's constant
- BRENT
- 1977
(Show Context)
Citation Context ... TEST2 if required. TEST2 has been run on several different machines. Knuth's computation of the continued fraction for ~/has been verified and extended using MP. Details will be published separately =-=[10]-=-, but it is worth noting that only 38 seconds of Ul108 time were needed to compute and verify the 372 quotients given by Knuth [31]. MP has also been used successfully to verify the orders of certain ... |

1 |
Knuth's constants to 1000 decimal and 1100 octal places
- BRENT
- 1975
(Show Context)
Citation Context ...es with several different wordlengths. Details are given in Table I. The program TESTV is a slight modification of TEST to allow variable precision, For checking purposes, 1000D results are available =-=[12]-=-. Table II shows how the execution time depends on the precision required. 2 2 Some of the results were checked using the package of Wyatt et al. [50]. MP was slgmficantly faster: the ratio of executi... |

1 |
MP users guide Tech Rep
- BRENT
- 1976
(Show Context)
Citation Context ...[31]. MP has also been used successfully to verify the orders of certain root-finding methods and to determine their asymptotic error constants [16]. 8. AN EXAMPLE The main program, EXAMPLE, given in =-=[9, 13]-=-, is intended to give a simple example of the use of MP. It computes 7r and exp(r(163)I/2/3) and prints them to 100 decimM places. The correct output is and 7r = 3.14159265358979... 1170680 exp0r(163)... |

1 |
Some high-order zero-finding methods using almost orthogonal polynomials
- NT
- 1975
(Show Context)
Citation Context ... to compute and verify the 372 quotients given by Knuth [31]. MP has also been used successfully to verify the orders of certain root-finding methods and to determine their asymptotic error constants =-=[16]-=-. 8. AN EXAMPLE The main program, EXAMPLE, given in [9, 13], is intended to give a simple example of the use of MP. It computes 7r and exp(r(163)I/2/3) and prints them to 100 decimM places. The correc... |

1 |
The complexity of multiple-precision arithmetic
- BIENT
- 1976
(Show Context)
Citation Context ...ivisions except by n. The iteration is n y~+l = y~ q- yj(1 -- xy~ )/~. .lIPROOT implements this iteration, starting with a small value of t and approximately doubling it at each step, as described in =-=[14, 17]-=-. Thus the time required by ~[PROOT is O(M(t)). MPREC implements the special case n = 1. Division and square roots require one additional multiplication after the computation of --1 X--I/2, x or respe... |

1 | A floating-point technique for extending the available precision - DECKER - 1971 |

1 |
High accuracy gamma function values for some rational arguments
- GAANT, BYRD
- 1968
(Show Context)
Citation Context ...MPEUL was checked using published values of ~ [31, 48]. The testing of MPEUL actually revealed an error in the most accurate published value of ~, [5, 6]. 5/IPGA54Q was checked using published values =-=[21]-=- as well as various identities. MPGASI and MPLNGM were tested using MPGAMQ and various identities. MPEI, 5,IPERF, 5IPERFC, and MPBESJ were tested for sufficiently large arguments by using both the asy... |

1 |
Algorithm 236. Bessel functions of the first kind
- GAVTSCm
- 1964
(Show Context)
Citation Context ... p and multiple precision x. If x is sufficiently large, Hankel's asymptotic expansion is used [2, sect. 9.2.5]; otherwise, either the power series [2, sect. 9.1.10] or the backward recurrence method =-=[22]-=- is used, depending upon how much cancellation occurs with the power series. The time required is O(tM(t)). 7. TEST RESULTS 5Iost testing has been done using Fortran V on a U]~ivac 1108. Gross errors ... |

1 |
Acceleration of series. Memo 304
- GOSPER
- 1974
(Show Context)
Citation Context ...nd other system reqmrements. a 36-bit word, b = 65536, t = 10, and the precision is sufficient to give correctly rounded 40D results. TEST calculates ~(3) from the rapidly converging series of Gosper =-=[23]-=- : 5 ~ (- 1)k(k!) 2 = (k 4 The other constants each require only a few calls to MP routines for their computation. TEST has been run successfully on various machines with several different wordlengths... |

1 | Computer Approx~matwns - HART, AL - 1968 |

1 |
Oeuvres de Charles Hermzte
- I-IERMITE
- 1908
(Show Context)
Citation Context ...14159265358979... 1170680 exp0r(163)l/2/3) = 640320.00000000060486 ... 6477590. EXAMPLE also computes exp0r(163) 1/2) = 262537412640768743.99999999999925007... and prints it to 90 decimal places. See =-=[25, 37, 39]-=- for the reason these numbers are interesting. EXAMPLE has been run successfully on the same machines and with the same compilers as TEST (see Section 7). ACM Transactions on Mathematmal Software, Vol... |

1 | Algorithm 34, Procedures for the basic arithmetical operations in multiplelength working - HILL - 1968 |

1 | Algorithm 72. Multiple integer arithmetic procedures in Algol - JONES - 1972 |

1 |
The Art of Computer Programming, Vol 2" Seminumerical Algorithms
- KNVTH
- 1973
(Show Context)
Citation Context ...fy as possible, we have followed most of the suggestions of Kernighan and Plauger [30]. 6. SOME ALGORITHMS USED IN MP 6.1 Addition This is straightforward, and is performed much as described in Knuth =-=[32]-=-. We use R*-rounding [15, 33] after postnormalization, with four guard digits. 6.2 Multiplication The classical O(t ~) method is used because it seems difficult to implement faster algorithms [29, 32,... |

1 | L Basic Q-precision arithmetic subroutines including input and output - LAWSON - 1967 |

1 | Q-precision subroutines for the elementary functmns and aids for testing single-precision and double-precision function subroutines - LAWSON - 1968 |

1 | L Summary of Q-precision subroutines as revised in October 1968. Tech Memo 211, Jet Propulsaon Lab - LAWSON - 1969 |

1 | Fortran programs for arbitrary precision arithmetic - MAXIMON - 1971 |

1 |
Collected Papers of Srinwasa Ramanu3an Cambridge U
- RAMANUJAN
- 1927
(Show Context)
Citation Context ...14159265358979... 1170680 exp0r(163)l/2/3) = 640320.00000000060486 ... 6477590. EXAMPLE also computes exp0r(163) 1/2) = 262537412640768743.99999999999925007... and prints it to 90 decimal places. See =-=[25, 37, 39]-=- for the reason these numbers are interesting. EXAMPLE has been run successfully on the same machines and with the same compilers as TEST (see Section 7). ACM Transactions on Mathematmal Software, Vol... |

1 | A multiple precision arithmetic package for the IBM 360/370 systems - KNOBLE - 1974 |

1 | The PFORT verffier Software--Practice and Experience - RYDER - 1974 |

1 | The testing of mathematical function software in a multi-machine enwronment - SCHONFELDER - 1975 |

1 | Fortran multiple precision, Pt. l, 2 - SPIRA - 1973 |

1 | A set of procedures making real arithmetic of unlimited accuracy possible within Algol 60 - TIENAR - 1966 |

1 | Multiple precision arithmetic design with an implementation on the Univac - CRARY |

1 | 52 CRARY, F.D. The Augment precompiler, Pt. I--User information - Center, Wisconsin, et al. - 1971 |