Multiple Precision Complex Arithmetic and Functions

David M. Smith Loyola Marymount University

This paper describes a collection of Fortran routines for multiple precision complex arithmetic and elementary functions. The package provides good exception handling, exible input and output, trace features, and results that are almost always correctly rounded. For best eciency on dierent machines, the user can change the arithmetic type used to represent the multiple precision numbers. Categories and Subject Descriptors: G.1.0 [Numerical Analysis]: Generalcomputer arithmetic; G.1.2 [Numerical Analysis]: Approximationelementary function approximation; G.4 [Mathematics of Computing]: Mathematical SoftwareAlgorithm analysis, eciency, portability General Terms: Algorithms, Performance, Reliability Additional Key Words and Phrases: Complex arithmetic, multiple precision, accuracy, function evaluation, oating point, Fortran, mathematical library, portable software

1. INTRODUCTION The ZM package is a collection of Fortran subroutines that performs oating point multiple precision evaluation of complex arithmetic and elementary functions. These routines use the FM package [Smith 1991] for real multiple precision arithmetic, constants, and elementary functions. Brents MP package [Brent 1978] did not support complex arithmetic, and Baileys more recent MP package [Bailey 1993; Bailey 1995] provides complex arithmetic and some complex elementary functions, and contains a Fortran-90 module dening multiple precision derived types. ZM also provides a Fortran-90 module that denes three multiple precision data types and provides the interface routines for overriding arithmetic operators and intrinsic Fortran-90 functions. This allows a program to declare variables as multiple precision real, integer, or complex, and then to describe operations on these variables using the normal syntax for arithmetic expressions. ZM versions are available for the numerical Fortran-90 intrinsic functions, and they provide good speed, rounding, and exception handling. They support exible input and output conversion, and have the capability for automatic tracing of arithmetic and functions.
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 prot or direct commercial advantage and that copies show this notice on the rst page or initial screen of a display along with the full citation. 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, to redistribute to lists, or to use any component of this work in other works, requires prior specic permission and/or a fee. Permissions may be requested from Publications Dept, ACM Inc., 1515 Broadway, New York, NY 10036 USA, fax +1 (212) 869-0481, or [email protected].

D.M. Smith

There are three multiple precision data types: FM IM ZM (multiple precision real) (multiple precision integer) (multiple precision complex)

The following table lists the operations dened in the interface module for the three derived types. For each of the operations =, +, -, *, /, **, .EQ., .NE., .GT., .GE., .LT., and .LE., the interface module denes all mixed mode variations involving one of the three multiple precision derived types and another argument having one of the types integer, real, double, complex, complex double, FM, IM, ZM. This insures that if A is type FM, mixed mode expressions such as A = 12 A = A + 1 IF (ABS(A).LT.1.0D-23) THEN are handled correctly. Not all the named functions are dened for all three multiple precision derived types, so the list below shows which can be used. The labels real, integer, and complex refer to types FM, IM, and ZM respectively. For functions that accept two or more arguments, like ATAN2 or MAX, all the arguments must be of the same type. In addition to these types, the three conversion functions TO FM, TO IM, and TO ZM accept character strings for input arguments, (e.g., TO FM(3.45)), and those functions also accept any of the machine precision data types integer, real, double, complex, or complex double. For converting a multiple precision number to one of the machine precision types, there are ve conversion functions: TO INT, TO SP, TO DP, TO SPZ, and TO DPZ. Table 1. = + /
.EQ. .NE. .GT. .GE. .LT. .LE. ABS ACOS real real real real real real real real real real real real real real

Available operations for multiple precision derived types

integer integer integer integer integer integer integer integer integer integer integer integer integer complex complex complex complex complex complex complex complex complex complex complex complex complex complex HUGE INT LOG LOG10 MATMUL MAX MAXEXPONENT MIN MINEXPONENT MOD MODULO NEAREST NINT PRECISION real real real real real real real real real real real real real real integer complex complex integer integer integer integer integer integer integer complex complex complex complex complex

Multiple Precision Complex Arithmetic and Functions

AIMAG AINT ANINT ASIN ATAN ATAN2 BTEST CEILING CMPLX CONJ COS COSH DBLE DIGITS DIM DINT EPSILON EXP EXPONENT FLOOR FRACTION real real real real real real integer real real real real real integer complex complex complex integer integer integer complex complex real real integer complex complex complex complex complex real real real real real integer complex complex complex complex complex complex RADIX RANGE REAL RRSPACING SCALE SETEXPONENT SIGN SIN SINH SPACING SQRT TAN TANH TINY TO FM TO IM TO ZM TO INT TO SP TO DP TO SPZ TO DPZ real real real real real real real real real real real real real real real real real real real real real real integer integer integer integer integer integer integer integer integer integer integer integer integer

complex complex complex complex

complex complex complex complex complex complex complex complex complex complex complex complex complex complex

DOTPRODUCT real

2. EFFICIENCY Prior to Baileys package, most multiple precision packages stored the numbers in integer arrays. This was faster for machines of the time, but Bailey reported that for some supercomputers and scientic workstations it was better to use double precision. His routines stored the numbers in single precision arrays and used double precision arithmetic to accumulate the multiple precision results. The advantage of using double precision arithmetic is that a larger base can be used for the multiple precision arithmetic. This reduces the number of array elements required to store a number with a given precision. The base is usually restricted to be less than the square root of the largest representable integer. Using integer arithmetic on a 32-bit computer, each word holds about four decimal digits, compared to seven for double precision. Multiplication, division, and functions have running times no faster than O(t2 ) for t words of precision unless t is large. These times may be faster by a factor of about three if double precision arithmetic is as fast as integer arithmetic. At commonly used lower precisions (20 to 60 decimal digits), the overhead involved in a general multiple precision package reduces this potential speedup factor. However, some machines are faster using integer arithmetic. To increase the exibility of both the FM and ZM packages, they have been written so that the type of arithmetic used for the basic operations can be easily changed. The user can select the type of arithmetic for the packages to use internally. The default is to use double precision. To produce an equivalent integer version, the user duplicates

D.M. Smith

the source code le and uses a text editor to make two (global) string replacements in the code. The details are given in the documentation at the beginning of the new FM package. The new version of the FM package incorporates several changes that support the ZM routines. There are some new routines, and most of the old ones are faster. The division routine uses a new algorithm [Smith 1996] that is about twice as fast as the original. In addition, the package now includes a collection of routines for integer multiple precision arithmetic. Here are some timing comparisons at 50 signicant digits and 1000 signicant digits that illustrate the eect of using dierent arithmetic on dierent machines. For a precision level of 50 digits, the integer runs used 14 digits in base 104 , while double precision used 8 digits in base 107 . The reason 14 digits are needed in base 104 is that normalization can cause the rst word to have only one signicant digit, so only about 53 signicant digits can be assumed to be present.

Table 2.

Complex arithmetic timing (seconds per call on a 604 Macintosh) 50 S.D. ZM Integer ZM D.P. 1000 S.D. ZM Integer ZM D.P.

add subtract multiply divide sqrt(z ) exp(z ) ln(z ) sin(z ) arctan(z )

.000 .000 .000 .000 .000 .000 .002 .000 .002

008 009 055 142 228 858 130 908 650

.000 .000 .000 .000 .000 .000 .002 .000 .002

008 008 047 074 186 842 190 861 620

.000 .000 .003 .011 .006 .062 .102 .062 .121

053 052 130 500 330

.000 .000 .001 .001 .002 .027 .049 .027 .053

037 038 200 690 380

Table 3.

Complex arithmetic timing (seconds per call on a 68030 Macintosh) 50 S.D. ZM Integer ZM D.P. 1000 S.D. ZM Integer ZM D.P.

add subtract multiply

.000 280 .000 303 .002 750

.000 480 .000 523 .003 600

.001 620 .001 720 .209

.003 380 .003 500 .160

Multiple Precision Complex Arithmetic and Functions

divide sqrt(z ) exp(z ) ln(z ) sin(z ) arctan(z )

.005 .009 .032 .073 .033 .088

790 500 700 100 900 900

.006 .011 .045 .107 .045 .127

120 600 200 000 800 000

.450 .360 3.680 5.780 3.680 6.620

.307 .279 2.900 4.810 2.860 5.440

These results are average times based on many calls to each routine. The same compiler and optimization settings were used in all cases, as well as the same set of input arguments. These tables are meant to compare the double precision and the integer versions on dierent machines, and to give a rough idea of the timing for various routines. However, the times can be sensitive to changes in compilers or hardware. With the more recent (604) machine, double precision is slightly better than integer at 50 digits. For 1000 digits, multiplication and division are much faster in double precision, and this causes the elementary functions to be much faster as well. The older (68030) machine has slower oating point performance compared to its integer arithmetic. At 50 digits the integer version is much better. The two versions are about equal around 800 digits, and the double precision version is faster for precision much higher than that. These timing tests were also done using a Sun workstation, 68040 and 601 Macs, and a Pentium PC. The results were similar to Table 2. Although most current machines run the double precision version faster than the integer version, for a machine with double-length integer arithmetic that is at least as fast as its double precision oating-point arithmetic the integer version would be better. 3. ACCURACY The complex arithmetic and function routines return results that are nearly always correctly rounded. Precision is raised upon entry to each routine so that the predicted error in the result before the nal rounding operation is less than 0.001 units in the last place (measured in terms of the original precision). The philosophy used for measuring errors is to assume all input arguments are exact. When zero digits are appended to the input arguments as precision is raised at the start of a routine, the values can still be assumed to be exact. This assumption is also used by Hull, Fairgrieve, and Tang [Hull et al. 1994]. For cases where the the input values really are exact, making this assumption assures that accuracy is not lost needlessly. Except for cancellation error in addition or subtraction, most algorithms lose accuracy slowly and predictably. A simple guess for how many guard digits are needed will usually prove correct. In rare cases, unexpected cancellation error will mean that the goal of 0.001 ulps of error cannot be achieved using the number of guard digits originally chosen. With complex arithmetic, even multiplication or division can suer from cancellation error. The following example using 5 digits in base 10 shows how ZM

D.M. Smith

routines handle this problem, although the actual routines would initially use more than two guard digits in this case. (0.63287 + 0.52498 i) (0.69301 + 0.83542 i). If precision is raised to 7 digits and the simple formula for multiplication is used, the result is 6.4000E-6 + 8.9253E-1 i. The real part is correct to only two signicant digits. The multiplication routine detects this loss of signicance, raises the precision to a higher level, and tries again. Using at least 10 digits gives the result 6.4471E-6 + 8.9253E-1 i. This is the correctly rounded answer. Similar cancellation can occur in functions such as ln( .77266 + .63483 i ) = 6.3022E-6 + .68778 i. 4. SOME ALGORITHMS USED IN THE PACKAGE Many of the basic formulas are the ones used in [Wynn 1962]. Hull, Fairgrieve, and Tang [Hull et al. 1994] remarked that the formulas for complex arithmetic and elementary functions appear deceptively simple. Many special cases must be handled, including cases where one or more parts of the input arguments are zero, where exceptions (underow, overow, or unknown) are generated as output, and where exceptions are encountered as intermediate results but the nal results are normal numbers. As noted by Hull, Fairgrieve, and Tang, a good exception-handling facility makes designing the routines much easier. The treatment of exceptions in the FM package simplies exception-handling for the complex functions. Because the basic FM real arithmetic is faster than the equivalent complex arithmetic, the elementary complex functions use real functions whenever possible. Some of the complex functions must check for too much cancellation error and on rare occasions must raise the precision and do the calculation a second time. 4.1 Multiplication For relatively low precision the multiply routine uses the simple formula (a + b i)(c + d i) = (a c b d) + (a d + b c) i. At higher precision the routine rst computes P = (a + b)(c + d), then applies the formula (a + b i)(c + d i) = (a c b d) + (P a c b d) i. This replaces one O(t2 ) multiplication by three O(t) additions. At low precision, a real multiplication takes about twice as long as a real addition. Because complex multiplication and division each require several real operations, precision must be raised to provide guard digits. The result is that complex multiplication is about seven times slower than real multiplication at low precision, and between three and four times slower at high precision.

Multiple Precision Complex Arithmetic and Functions

4.2 Division This method is due to [Smith 1962]. If |c| < |d|, let Q = c/d and then use a+bi aQ + b bQ a = + i. c+di cQ + d cQ + d Otherwise, let Q = d/c and use a+bi bQ + a b aQ = + i. c+di dQ + c dQ + c This requires three multiplications, three divisions, and three additions. Making |Q| < 1 prevents overow in the intermediate results. Complex division is about seven times slower than real division at low precision, and about ve times slower at high precision. 4.3 Square Root The algorithm used is due to [Friedland 1967]. |a| + To nd c + d i =

a + b i, let

t= Then if a 0, c = t, Otherwise, c= b . 2t

a2 + b2 . 2

d=

|b| , d = Sign(b) t. 2t This is about three times slower than a real square root.

4.4 Exponential and Logarithmic Functions The exponential function uses the formula ea+b i = ea cos(b) + ea sin(b) i. One of the trigonometric functions is computed and the other can be obtained quickly using an identity. The real sine function is about as fast as the exponential, so the complex exponential takes slightly more than twice as long as the real exponential. For the logarithm, with |a + b i| = a2 + b2 and arg(a + b i) = atan2(b, a), then ln(a + b i) = ln( |a + b i| ) + arg(a + b i) i. The complex logarithm takes slightly more than twice as long as the real logarithm. The general exponential, y x , uses three complex operations: logarithm, multiplication, exponential. The integer power function, y n , uses the binary method with complex multiplication. The rational power function, y k/n , uses Newton iteration to get y 1/n , followed by a call to the integer power function. 4.5 Trigonometric Functions cos(a + b i) = cos(a) cosh(b) sin(a) sinh(b) i

D.M. Smith

sin(a + b i) = sin(a) cosh(b) + cos(a) sinh(b) i sin(2a) sinh(2b) tan(a + b i) = + i cos(2a) + cosh(2b) cos(2a) + cosh(2b) These functions are slower than their real counterparts by factors ranging from 2.5 or 3 at low precision to slightly more than 2 at high precision. The inverse functions call the complex logarithm and complex square root. The two equivalent forms are used to avoid cases where one or the other suers from cancellation. cos1 (z ) = i ln z + i sin1 (z ) = i ln tan1 (z ) = i ln 2 1 z 2 = i ln z i 1 z2

1 z 2 + i z = i ln i+z iz

1 z2 i z

These are two or three times slower than the real versions. 4.6 Hyperbolic Functions cosh(a + b i) = cosh(a) cos(b) + sinh(a) sin(b) i sinh(a + b i) = sinh(a) cos(b) + cosh(a) sin(b) i sinh(2a) sin(2b) tanh(a + b i) = + i cosh(2a) + cos(2b) cosh(2a) + cos(2b) Timing is similar to the trigonometric functions. 4.7 Conversion and Utility Functions Routines are provided for complex absolute value, argument, conjugate, real or imaginary parts, integer or nearest integer values, and conversion between real and complex FM formats. The input conversion routines do free-format conversion from character format to ZM format. Input strings can be in (1.23 , 4.56) form or 1.23 4.56 i form. The two numeric parts can be in any integer, xed, or exponential format. If one part is omitted, input of the form 3.45 or 5 i is correctly converted to ZM format with the missing part set to zero. The output conversion routines can be set for the desired complex format, the number of digits to display, and the real format to use for the individual parts of the complex number. 5. TESTING To test the accuracy of the functions in this package, many random arguments were generated and results compared at a sequence of increasing precisions. As independent tests for the functions, ZM results were compared to values tabulated in [Abramowitz and Stegun 1965] and to values generated by the Mathematica computer algebra system [Wolfram 1991]. Tests were also made comparing

Multiple Precision Complex Arithmetic and Functions

ZM results at low precision with values generated by the complex Fortran intrinsic functions.
REFERENCES Abramowitz, M. and Stegun, I. 1965. Handbook of Mathematical Functions. Dover, New York. Bailey, D. 1993. Multiprecision translation and execution of fortran programs. ACM Trans. Math. Softw. 19, 288319. Bailey, D. 1995. A fortran 90-based multiprecision system. ACM Trans. Math. Softw. 21, 379387. Brent, R. 1978. A fortran multiple-precision arithmetic package. ACM Trans. Math. Softw. 4, 5770. Friedland, P. 1967. Absolute value and square root of a complex number. Comm. ACM 10, 665. Hull, T., Fairgrieve, T., and Tang, P. 1994. Implementing complex elementary functions using exception handling. ACM Trans. Math. Softw. 20, 215244. Smith, D. 1991. A fortran package for oating-point multiple-precision arithmetic. ACM Trans. Math. Softw. 17, 273283. Smith, D. 1996. A multiple-precision division algorithm. Math. Comp. 66, 157163. Smith, R. 1962. Complex division. Comm. ACM 5, 435. Wolfram, S. 1991. Mathematica: A System for doing Mathematics by Computer. AddisonWesley, Redwood City, Calif. Wynn, P. 1962. An arsenal of algol procedures for complex arithmetic. Bit 2, 232255.