LIBRARY OF THE UNIVERSITY OF. ILLINOIS AT URBANA-CHAMPAIGN 510.84- UL6r no. 403-4-08 Cop. 2 The person charging this material is re- sponsible for its return to the library from which it was withdrawn on or before the Latest Date stamped below. Theft/ mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. To renew call Telephone Center, 333-8400 UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN DFC 3 RE B L161— O-1096 Digitized by the Internet Archive in 2013 http://archive.org/details/multipleprecisio404flec t r Report No, iOk sv/szz^c. coo- ii+6 9-0 167 June, 197O MULTIPLE PRECISION COMPUTATION OF LEGENDRE FUNCTIONS by Ruth Ann Potts Fleck THE LIBRARY OE THE NOV 9 1972 UNIVERSITY OF U.LINOIS AT URBANA-CHAMPAiGN DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMP mw NA, ILLINOIS MULTIPLE PRECISION COMPUTATION OF LEGENDRE FUNCTIONS BY RUTH ANN POTTS FLECK A.B., Ripon College, 1966 THESIS Submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science in the Graduate College of the University of Illinois, 1970 Urbana, Illinois Ill ACKNOWLEDGMENT The author wishes to thank Professor L. D. Fosdick for his supervision and guidance in the preparation of this thesis. Also sincere appreciation is extended to this author's husband for his patience, assistance and encouragement during the writing of this paper. The computational work described in this thesis was performed on the I.B.M. 360 75 at the University of Illinois. This research was supported by TR US AEC Contract 1U69. IV TABLE OF CONTENTS Page I . INTRODUCTION 1 II. CALCULATION OF ASSOCIATED LEGENDRE FUNCTIONS 5 III . MULTIPLE PRECISION PROGRAMS 16 IV. DESCRIPTION OF TESTS AND EXPERIMENTAL RESULTS kl V. SUMMARY AND SUGGESTIONS FOR FURTHER STUDY 63 BIBLIOGRAPHY 65 APPENDIX 67 I. INTRODUCTION In the solution of problems arising in applied mathematics, it is often necessary to know the value of certain functions. Many of these functions are difficult or extremely tedious and time consuming to evaluate and so it has become traditional for tables of often used functions to be published. The increased use and efficiency of the digital computer has both solved and created many problems in this area of function evaluation. It is now possible to create tables much faster than by human calculation and consequently to tabulate functions for many more values of the argument and to a greater number of decimal places. However the use of these tables is now more difficult. A physicist writing a computer program to solve a problem does not want to stop the program at the point where he needs a function value, print out the argument, go to a table to look up the function value, read it into the computer and resume running the program. It is impractical to have the tables themselves stored in the computer because of the storage space consumed and to have the tables on punched cards implies large volumes of physical storage space for the cards plus time to read the cards each time a program is run. A good solution would seem to be storage of tables on magnetic tape which requires less physical room and is much faster on input, in fact the whole table would not have to be read in every time if there were section markers on the tape. Many decisions would have to be made to make this system workable in a universal sense. Perhaps the most difficult is how to store the numbers 2 on tape so they can be read by many different computers. Another im- portant consideration is for what values of the argument should the function be evaluated so that the table is adequate, with the use of interpolation, for the greatest number of cases. Another approach to this problem is simply to generate the values of the table which are needed as the program progresses. Here the disadvantage is the complexity of the code which may be required to evaluate the function and the time required for the evaluation. The scientist may not be well versed in methods of function evaluation, particularly in terms of error estimation. Most programs which are written to evaluate functions are not transferrable between computers because of differences in machine languages. One step toward the solution of some of these problems has been attempted here through the development of a group of FORTRAN programs which calculate various associated Legendre functions for arguments greater than 1. These programs are machine independent in the sense that FORTRAN is a language available on most computers and the modifi- cations necessary to run on a different computer are relatively simple. The programs can also be considered as the first step toward a tabulation on tape since it would only be necessary to make the decisions mentioned above as to values of the argument and format of the tape and then ob- tain the function values from the programs described here. Another problem concerned with tables of functions is the number of decimals to which the function should be calculated. Although a scientist usually does not require large numbers of decimals in specific 3 numerical problems , there are many areas where they are needed in order to have any significance in the final answer. This subject has been discussed at length in the Query-Reply section of Mathematical Tables and Other Aids to Computation . It is particularly true in the case of associated Legendre functions for arguments greater than 1 where the function value is increasing or decreasing at a rapid rate compared with an increase in the value of the argument. Therefore it was con- sidered desirable to allow for arbitrary precision in developing the programs since the type of use to which the function values will be put is the determining factor in the precision needed. FOOTNOTES "Ql, QR1: Tables to Many Places of Decimals," MTAC 1 (19^3), pp. 30-31. "QR2: Tables to Many Places of Decimals," MTAC 1 (19^3), pp. 67-68. "QR3: Tables to Many Places of Decimals," MTAC 1 (19^3), p. 99- "QR4: Tables to Many Places of Decimals," MTAC 1 (19^3), pp. 99-100. "QR2T: Tables to Many Places of Decimals," MTAC 2 (19^7), pp. 226-228. CALCULATION OF ASSOCIATED LEGENDRE FUNCTIONS The associated Legendre functions of the first and second kind, P (z) and Q (z), are solutions of the differential equation 2 u 2 (l-z^w'^zw' + [a(a+l) - ^ — - ] w=0 . 1-z These solutions are usually represented in terms of the hypergeometric function. We shall consider here only the cases P (x) and Q (x) where m is a non-negative integer, x is real and greater than 1 and a # -1, -2, -3, ... The Legendre functions can then be represented 2 by the following definite integrals. ^m, N r ( q +m+1 ) / , / > > ., , N a .,. P (x) = ' / ^ s / (x + /x*-l cos t ) cos mtdt uu rP-i \ i -, \ m r(a+l) / cosh mt ,. Q a (x) = ( " 1} rfe^l) V ~= — ^KL dt (x + /x z -l cosh t) These integrals arise in many areas of applied mathematics, among them 3 aerodynamics and radiation field studies. Both r (x) and Q (x) satisfy identical three-term recurrence a a k relations of varying order m and degree a. Of particular interest will be 6 U* +1 (x) + -^^ lAx) + (m+a)(m-a-l)U m " 1 (x)=0 2.1 Ot r~ "h ~ Ot Ot vx -1 (a-m+1) U m ^(x)-(2a+l)xU m (x) + (a+m)U m , (x)=0 2.2 a+1 a a-1 where U stands for either P or Q. The Legendre functions for x > 1 can be generated in a straight forward manner from these recurrence relations, they have in fact been used for previous table computation, since the recursion is stable in the forward direction for P^x) and in the backward direction for Q (x). a a However, this method requires that starting values are available and these are difficult to obtain in some cases, particularly when a is not an integer or half-integer. The procedure to be described here, which replaces the condition of knowing two values of the function to begin with by a much more general condition in which at most one value must 6 be known, was developed by Gautschi based on the relationship between three-term recurrence relations and continued fractions. The algorithm is mathematically equivalent to the backward recurrence algorithm of J. C. P. Miller which has been used to obtain Legendre functions although not as extensively or efficiently as by Gautschi. A continued fraction can be related to every three-term recur- rence relation in the following manner. Let A and B be solutions of n n y ...n + a ±1 y + b xJ i = ° (n=0,l,2,. . . ) 2.3 n+1 n+l J n n+l J n-1 ' ' ' determined by the initial values A_ ± = 1, A Q = 0; B_ 1 = B Q = 1 . The A and B are then the numerators and denominators of the continued n n fraction 1 2 3 or equivalently , a_ — ^p - ^^j" Alternatively, let r n = ^ (n=-l,0,l,2,...) n The recurrence relation (2.3) can then be written as Vi r + a L + — = 0, n n+1 r n-1 -Vi or n-1 a , _ + r n+1 n Repeating this for increasing n gives us ln_ _ -Vl b n + 2 b n+3 n - X " y n-l a n + l" a n+2" W 2.1+ In order to apply this relationship we need information on the convergence of the continued fraction and for what particular solution the ratios are to be formed. For this purpose we intro- duce the concept of distinguished solutions. A distinguished solution , f , is a solution of a difference equation (2.3) with the property that lim f n->- °° — = y n for any other linearly independent solution y of the difference 8 equation. The following theorem based on a result of Pincherle is then available. Theorem 1 The continued fraction (2.U) converges if and only if the re- currence relation (2.3) possesses a distinguished solution f , with f f 0. In case of convergence, moreover, one has £n_ = - Vi V> V3 . (n=0 1 2 } n-1 n+1 n+2 n+3 provided f i for n = -1,0,1 ,2. . . . Now write (2.3) in the form y_,,+ay +by =0(n = 1,2,3,...) 2.5 n+1 n J n rr n-1 v ' ' ' ' and assume a nonvanishing distinguished solution, f , specified uniquely by one condition I m=0 X f m m s (A Q * 0, s i 0), 2.6 where s, and X , A , ..., are given and the series converges and is normalized so that \ n = 1. If A =0 for all m > this is the special m case f = s. The algorithm for computing f for n = 0, 1, 2, . . . , N, based on Theorem 1 is then given as follows: Step 1 : Select an integer v > N, and let = (n = 0,1, Step 2 : Calculate f^ (n = 0,1,..., N) according to: n (v) (v) _ n r = , r = t — r v n-1 . (v) a +r n n V .(v) . s (v) . (v) (A + s (v) , v n-1 n-1 n=v,v-l ,N) i + s v n n I .1, ^ 2.7 .(v) 1+s, (v) ,(v) 'n+1 = r (v) f (v) (n=0,l,...,N-l) n n Step 3 : If the N+1 values of f^ obtained in Step 2 do not agree with the current values of to within the desired accuracy, then redefine (v) (v) (v) by = f (n = 0,1,..., N), increase v by some integer, say nn n 5, and repeat Step 2; otherwise accept f as the final approximations to f (n = 0,1,. . . ,N). n The convergence of this algorithm is stated by Gautschi in the following theorem. 10 Theorem 2 Suppose the recurrence relation (2.5) has a nonvanishing distinguished solution, f , for which (2.6) holds. Let g he any- other solution of (2.5). Then the algorithm (2.7) converges in the sense lim f U) = f n n if and only if lim _v±i y- A = o v -> oo °v+l m=0 is satisfied. One remaining tool for applying this algorithm, coming from the asymptotic theory of difference equations, -will allow us to show that a given solution of a three-term recurrence relation is a distinguished solution. The asymptotic structure of the general solution of a difference equation (2.5) whose coefficients satisfy a ^an , b ^bn (ab ^ 0; a, B real; n -> °°) depends on the Newton-Puiseux diagram formed with the points P (0,0), P 1 (l,a), P 2 (2,B). This is PqP-^ if P 1 is above or on P Q P 2 , or P Q P 2 otherwise, a is the slope of P P (a=a) and x is the slope of P-,P ? ) . Gautschi now makes use of a theorem from work by Perron and Kreuser. 11 Theorem 3 (a) If the point P is above the line segment PJ?^ (i.e. o>t), the difference equation (2.5) has two linearly independent solutions, y _ and y _, for which n,± n,2 y n+l,l a y n+l,2 b t , v z — ^ - an , 2 — ^ - — n (n-> °°) . n,l n,2 (b) If the points P , P , P are collinear (i.e., a=x=a), let t ,t be the roots of t +at+b=0, and |t | >„ |tp|. Then (2.5) has two linearly independent solutions, y . and. y ~, such that n,l ^n,2' n+1,1 ^ t,n a °n+l,2 'n.1 ' r ' ' y n,2 ^ t n (n^ °°) , provided|t | >|t |. If |t. | = |t | (in particular, if t ,t„ are complex conjugates, then lim sup y. (n!)' 1/n for all solutions of (2.5) (c) If the point P lies below the line segment P n P p then 1/n - /FT lim sup n->- °° y. (n!) n 3/2 for all solutions of (2.5) 12 12 It can be shown that in case (a), and case (t>) where |t | > |t_|, the solution f = y is a distinguished solution of (2.5). Also case (b) with |t | > |t p | gives us ^+1 / rs\ lim = t (r = l, or r= 2), ''n where r = 2 for the distinguished solution, and r = 1 for any other solution. To illustrate how this algorithm may be applied to the calculation of Legendre functions, consider the case of recurrence with respect to the order m where the degree a is not an integer. From (2.1) we see that r (x) and Q (x), as functions of m, are solutions of a a y -, + — — • y + (m+a)(m-a-l)y =0 2.8 °m+l /~2 — ~ °m /J m-1 /X " 1 (m=l,2,3,...) . The Newton-Puiseux diagram of (2.8) is a straight line segment with slope 1. The roots of the equation 2 2x XT + — t + 1 = /x+1 . -1 and since x > 1, |t | > |t p |. Thus by case (b) of Theorem 3 we know that (2.8) possesses a distinguished solution f with 13 f i • m+1 m r 2 nr* °° m P?(x) Let F = / . .n \ — = P / ■-, \ / [x+/x 2 -l cos tJ a cos mt dt, F +1 P^U) Then ^ (m-* °°) m mr (x) a which has a finite limit, t or t_, as m.-* °°. If it were t n , then F 12 1 ' m 1 would tend to °°, since |t | > 1, but this is impossible by the definition of F . Therefore m P m+1 (x) lim = t m-> °° mr (x) a so Fix) is the distinguished solution of (2.8) and can be calculated for m = 0,1,2,... by the algorithm described, with the condition (2.6) 1^ given by the following identity, oo . . P a<*> + 2 V -&&?, ftM - (x ♦ /x-^T f. 2.9 ) r ( a+m+l ) J " a m=l In the actual program, the algorithm is applied to the recurrence Ik f + 2mx_ F + m±^_ m+1 (m + a + l)^H m m+a+1 m - 1 P^x) ct which has F = —, — r — as a distinguished solution. This simplifies m r(a+m+l) & * the identity (2.9) to - (x +' / x 2 -l ) a + 2 A F m ' r(a+l) oo I m=l We then have A = 1 and A = 2 (m = 1,2,...) in our algorithm and can apply Theorem 2 to show convergence. The application of the algorithm to other cases of the Legendre functions is similar and will be indicated in the next chapter as part of the description of the FORTRAN programs. It should "be noted that this theory is valid for x = z, a complex number outside the interval (0,1 ) with Re z > 0, although the programs here restrict x to a real number. Also of interest is the fact that the algorithm as applied above to the Legendre functions P (x) does not require the knowledge of any starting values which simplifies computation. 15 FOOTNOTES 1. M. Abramowitz and I. A. Stegun, ed. , Handbook of Mathematical Functions , (Dover, New York, 1965), P- 332. E. W. Hobson, The Theory of Spherical and Ellipsoidal Harmonics , (University Press, Cambridge, 1931), pp. 89-118, 178-292. 2. H. Bateman, Higher Transcendental Functions , Vol. 1, (McGraw-Hill, New York, 1953), p. 157- 3- W. Gautschi, "Algorithm 229 Legendre Functions for Arguments Larger than One," Comm . Assoc . Comp . Mach . 8 (1965), p. U91. k. Bateman, pp. l60-l6l. Abramowitz and Stegun, pp. 333-33^. 5- NBS Mathematical Tables Project, Tables of Associated Legendre Functions , (Columbia University Press, New York, 19^5). 6. W. Gautschi, Strength and Weakness of Three-Term Recurrence Relations , (Unpublished). 7. A Rotenberg, "The Calculation of Toroidal Harmonics," Math . Comput . Ik (i960), pp. 27^-276. Gautschi, Strength and Weakness , p. 11. Ibid . , p. 21 . Ibid . , p. 2k. Ibid . , p. 17. Ibid . , pp. 15-16, 18. Bateman, p. 166. 16 III. MULTIPLE PRECISION PROGRAMS In order to perform computations in multiple precision in a digital computer, it is necessary to define a system for storage of the multiple precision numbers and to have available subroutines to perform ordinary arithmetic operations such as add, multiply, etc. and also to evaluate common functions such as cosine or natural logarithm, since the regular operations and functions are, by implication of the name multiple precision, less accurate than needed. The basic arithmetic programs described here are a translation into FORTRAN by Henry C. Thacher, Jr. of Argonne National Laboratory of an ALGOL algorithm. Their great advantage is that, being written in FORTRAN, they are easily transferred from one machine to another while most packages for multiple precision arithmetic are written in the assembly language of a specific computer and thus require a considerable expenditure of effort to use elsewhere. This advantage is somewhat offset by a large increase in the computer time necessary to run a specific program. Each multiple precision number is stored in an integer array of arbitrary length K. This arbitrary array length is what allows the user to specify the precision desired. The number is stored in the N\b 2*N form a*(10 ) where N is chosen so that 10 does not cause integer overflow. In our specific applications on an IBM 360 the maximum integer is 21^7^836^7, and so N is set equal to k. The number is N ;tandardized so that 1 .< a < 10 . The multiple precision number ' ' 1.0 ) ) is then stored as follows: IT INTEGER**! X(K) X(l) = sign of the number: +1 if positive -1 if negative if zero (In this case, X(2) to X(K) are disregarded) X(2) = exponent b of 10 N X(3) = integral part of the mantissa, 1 ,< X(3) < 10 X(U) to X(K) = fractional part of the mantissa, stored N decimal digits per location Numbers stored in this manner have at least N*(K-3)+l significant decimal digits. The maximum magnitude of a number which we can store on the 360 with r = k is 9999.99 . . .*(lo' t ) 2:LU7U8361,7 - 9.99 ...no 858 "* 591 . This greatly extends the range of ordinary real numbers where the 75 maximum magnitude is approximately 10 The arithmetic operations are performed in floating point decimal arithmetic by a group of 13 subroutines. Listings and detailed descriptions of use are included in the Appendix. Since each sub- routine which performs operations on the multiple precision numbers needs internal arrays for the temporary storage of numbers , it is found convenient to set a maximum value on K, rather than to pass each temporary array as a parameter. This maximum value was chosen as 50, which limits the number of decimal digits to U*(50-3)+l=l89, which was considered adequate for our purposes. This maximum K could be increased in the basic arithmetic subroutines by merely changing the DIMENSION statements, but would involve increasing the precision of certain constants in the function subroutines such as cosine, etc. 18 LNGCNS - This is the first subroutine of the multiple precision package. It initializes certain constants which are used by many of the other subroutines and stores them in labelled COMMON. Constants which must be supplied are N (=h in our case) as described above, W(=8) the number of significant decimal digits in the representation of a real number, rounded up to the next higher integer if it is not an exact integer, and E(=75) the maximum allowable decimal exponent of a real number, rounded down to the next smaller integer if it is not an exact integer. This subroutine must be called before any of the basic arithmetic or function subroutines are used. COPY - The parameters are A, B, and K where A and B are multiple precision numbers of array length K. This subroutine performs the operation B=A. SET - This subroutine performs the operation A=I, where A is a multiple precision number and I is an ordinary integer. An error return occurs if the array A is not long enough to hold the integer I. LNGTHN - This subroutine performs the operation A=X, where A is a multiple precision number and X is real (single precision). The accuracy of A is limited to single precision, i.e. the number of significant digits in A is the same as X, the remainder of the array being filled with O's. I - A multiple precision number A is converted to a single pre- cision number X. An error return occurs if overflow would be produced. ADJUST - A subroutine called internally by ADD and SBTRCT. It is essentially equivalent to a right shift. 19 READJT - A subroutine called internally by ADD and MLTINT. It re- adjusts a multiple precision number into standard form. ADD - This subroutine performs the operation C = A + B, where A, B, and C are multiple precision numbers. SBTRCT - Given A and B, this subroutine performs C = A - B, where A, B, and C are multiple precision numbers. MLTPLY - This subroutine performs the operation C = A * B where A, B, and C are multiple precision numbers. MLTINT - Given A, a multiple precision number, and I, an ordinary N integer less than 10 , this subroutine performs the operation B = A * I. It should be used instead of MLTPLY when possible since it is much faster (the time to execute MLTINT is proportional to K, the array length of a multiple precision number, while the execution time of MLTPLY is proportional to K ) . An error return occurs if I >, 10 M . RECIP - This subroutine finds the multiple precision reciprocal, B, of a multiple precision number, A, by an iterative procedure based on Newton's method for finding the root of a function, in this case 1/B - A = 0. B is obtained by an ordinary division, 1.0/AS, where AS is the single precision equivalent of A obtained from the sub- routine SHORTN. Therefore the regular system failure for division by zero will occur if A = 0. The exponent is handled separately so that underflow or overflow will not occur when A is within the limits of a multiple precision number, but AS would not be within the limits of an ordinary real number. The iteration proceeds by the formula 20 B ., = B * (2 - B *A) and B , n is accepted as the final approximation n+1 n n n+1 if B , = B or B , = B n+1 n n+1 n-1 DIVIDE - Given the. multiple precision numbers A and B, this subroutine will perform the operation C = A/B by using the subroutines RECIP and MLTPLY. For programming considerations, note that DIVIDE, by- using RECIP, requires an iterative procedure while MLTPLY does not, so the programmer should avoid divisions whenever possible. For example, rather than dividing a multiple precision number by 2, multiply it by .5. Besides these basic arithmetic operations, it is also necessary to be able to evaluate certain common functions. For the applications here it was necessary to evaluate vA, In A, A , A , e , cos A, and r(A), where A and B are multiple precision numbers and I is an ordinary integer. The first two programs were written by Thacher and modified only slightly by this author. The remaining are original programs. SQRT - This subroutine finds Y = i/x where X and Y are multiple precision numbers, by Newton's method for finding roots applied to the function 2 1 - X/Y = 0. This gives an iteration of Y.,, = Y. +6., where l+l l l 2 6. = Y. (X - Y. )/2X, which minimizes divisions. Y is obtained from l l i ' the single precision square root routine and Y is accepted as the final approximation if 6. =0 or exp(Y ) - exp(6.) > K - 3 (exp the exponent, b, of 10 in the multiple precision representation of a number stored in an array of length K). If exp(Y. ) - exp(6. ) = and 6. = -6. . then Y.,, = Y. + 6./2 is accepted as the final l l-l l+l li ^ 21 approximation to Y. If X < 0, the normal system error will occur in the single precision square root routine. An error message is printed if none of the conditions described are met for i ^ 50. LNGLOG - Y = In X is evaluated by LRGLOG for X and Y multiple pre- cision numbers. An error return occurs if X is nonpositive. The series r 1 i i x-i \ 2 k+1 k=0 --iMar /*-(&T' is used for .3^ X s < 3. For other values of X, the range is reduced p by using two constants, a tabulated value of In 10, and k * In 10. The series is summed as follows: -K*9 V S. = S. _ +6. where l l-l l x - J— /x-l\ i ' 2i+l \ X+l ] 2i+l until 5. = or exp(S. ) - exp(6.) >^K - 2. If this criterion is not met for i < 1500 an error message is printed. OUTPUT - This is a subroutine used internally by LSQRT and LNGLOG to print error messages. 22 EXPO - This subroutine evaluates Y = X where X and Y are multiple precision numbers and I is an ordinary integer. Y is obtained by multiplying X by itself |l| times and then taking the reciprocal of the result if I < 0. An error message is printed if X = I = 0. LEXPO - This subroutine evaluates Y = x where X, Y, and A are multiple precision numbers. If A is an integer and less than or equal o to 20*10 , the subroutine EXPO is used as it requires less execution time. Subroutines EXP and LNGLOG are used to find Y from the equation A*ln X Y = e .An error message is printed if X = A = or if X < and A is not an integer. LEXP - Given the multiple precision number X, this subroutine will find Y = e X for X v < 20*10 8 . For < X < .5, e X is evaluated from the continued fraction 3 X e = b x + b 2 + b 3 + 1_ X X X X X X 1- 1+ 2- 3+ 2- 5+ 2- by the following algorithm u 1 = l v = w = — 1 1 b. u i+1 = . v i+i = V^+i-U^+i-V^i+i 1 + i+1 b."b. ., i i+I u, (1=1,2,3,...) 23 X A n A i If we write e = lim :r— then w. = r-~ and i is increased until B 1 B. n-> °° n i exp(w ) - exp(v ) >, K - 2. Then w. is taken as the final approximation to e . If this criterion is not met for i ^ 1000 an error message is printed. An error message is also printed for |x| > 20*10 . For other values of X, the range is reduced by using a tabulated value of e . I!0S - This subroutine finds Y = cos X for X and Y multiple precision numbers. X is first reduced to the interval v < X >< 2tt and then, using trigonometric reduction formulas, cos X is found by evaluating one of the following infinite series for >< X ,< r- . cos X = 1 - — + jj- - £j- +... , r , r A A A , sin X = X - — + jj- -—+... The reductions are made using a constant value of 2tt obtained from a tabulated value of ir. The series is summed until exp(Sum) - exp(term) ^ K - 2. If this criterion is not met by the 1000th term an error message is printed. The series is summed ten terms at a time with the addition made from the smallest value to the largest in each group of ten to cut down on loss of significant figures. LGAM - Y = r(X), where X and Y are multiple precision numbers of array length K, is evaluated to D significant digits. D should be less than or equal to (K-3)*^+l and X must be greater than 0. r(X) is computed from the equations 2k n c, In r(X+j) = (X+j -h * ln(X+j)-(X+j)+ lnv^+ V X , . 2 £i (x + j) 2l_1 3.1 (-1) 1 " 1 B. where c. = . . . . — — r— and B. is the ith Bernoulli number, l 2i(2i-l ) i In r(x+j) and r(X) = T-r: • 3.2 J i=l The values of j and n are chosen so as to give the desired D significant decimal digits. 7 The first equation above comes from the Sterling asymptotic expansion In r(t) = (t-^-)lnt- t + ln/27 + 0(t) 0(t) " 2- 2i(2i-D ^1 + R n (t) ' i=l l a with R (t) < J?" 1 " 1 , Wo — -,• -i — - for t > . n v ' 2(n+l)(2n+l) 2n+l x» It is seen from this that the error in r(X) will depend on n and t. (t) decreases as t increases it is only necessary to determine the minimum value t , corresponding to a given n, for which R (t) is small enough (i.e. R (t) < 10~ ) and then R (t) will also be small 25 enough for all t >, t . For a precision of up to 50 decimal digits, o Filho and Schwachheim determined, by a rough empirical relation, that computation time for In r(t) could be minimized by using n=20 and t = J if D < IT or t = D-10 otherwise. These values of n and t are used by LGAM for D ^ 50, but to extend D to 189, the maximum D since K is limited to 50 in these programs, it was desired to find the values of n so that the formula t = D-10 would still apply. Since we wanted to find n such that n+1 2(n+l)(2n+l) , t *x2n+l .< 10 -(t +10) 3.3 graphs were made of log y = (2n+l)log t - t - 10 - loj B n+1 2(n+l)(2n+l) for various n. For the values of t where log y >, the inequality (3.3) is satisfied. TABLE 1 is a summary of the values of n and t* chosen for use in LGAM. TABLE 1 1 to IT 20 18 to 50 20 51 to 80 31 8l to 110 kl 111 to 137 51 138 to 165 61 166 to 189 71 7 D-10 D-10 D-10 D-10 D-10 D-10 26 Given D, LGAM selects the proper n from TABLE 1. Then if X is greater than or equal to the corresponding t given in TABLE 1, Y = r(X) can be computed directly from (3.1). For X less than the corresponding t , equations (3.1) and (3.2) are used with j the least integer which satisfies X + j >, t . The constants needed in equation (3.1) are lnv2v, which was 10 calculated from tabulated values of In 2 and In tt , and the c. , ' i which were calculated using a table of the numerators and de- nominators of the Bernoulli numbers . The function subroutines were not extensively tested individually since the accuracy of the programs which calculate the Legendre functions using these function subroutines reflects their accuracy. (The tests of the Legendre programs are described in the next chapter. ) However it was considered informative to check the function 12 subroutines using a few known constants. Two tests were run for each subroutine. Table 2 shows the relative error of the computed answer. D is the number of significant decimal digits which can be stored in a multiple precision number of array length K. The seven subroutines which calculate the Legendre functions are translations of ordinary precision ALGOL programs written by 13 Gautschi into multiple precision FORTRAN programs using the multiple precision subroutines described above. The logic of the programs remains exactly the same with one exception. The starting value of Ik 'n algorithm (2.7) was determined empirically by Gautschi 27 o CO LA H II II ^ Q CO CO CO t— CO vo CO t- rH t— ^t -H; 3 CO CO CO CO CO CO co On co CO CO H H H H H H H H H H H H o O O O O O O O O O o O H H H H H H H H H H H rH * * * * * * * * * * * * o O H -d- CO H O J- c— O O CO rH CM -^ H J- rH ir\ VO on C\J CO J- CM ON t— CO QO t— CO NO C— C— O CO LTN LTN ON -=r -d" -=t -=1- _=r _* -3" -3- _=t- LTN ,3- _=!- -=1- o -J- H H H rH H H H rH H H rH H H -=t r-l 1 1 1 1 I 1 1 1 1 1 1 1 1 II II o o o o o O O O o o o O O W M O H H H rH H H H rH H rH H H H * * * * * * * * * * * * * CO o o ^1- CO -=f ON LA CO t— CM ON on H on -3- LTN H CO H LA J- CO CM J" LTN ffi o K On P-* o o W on rH ll II w « « > ON t— CO O O O H H H O O o H H rH * * * CO O o H rH -H t— CO t— t— t— ON CO LTn r— O o O o o O O o o H H H rH rH H H H H O O O o O O O O o H H H rH H rH H H H * * * * * * * * * CM CM on LTN VO LTN CO LTN on on vo t— LA -3" H -3" CM t— CO O ON VO CM NO 1 II II O W Q H * CM on t>- CO CO t— t~- VO t— t— ON ON LTn la NO VO VO VO VO VO VO VO VO VO vo vo o o o o o o O O O O o o rH H H H rH H rH H H H H H * * * * * * * * X * * * o o rH on H -3- LTN CM LTN VO H H H H -=f CM H CM LTN CM on cm on h CO ON CO CO t— t— co CM CM CM CM CM CM CM CM CM CM on CM CM CM o ON 1 1 1 1 1 1 1 1 1 1 1 1 1 1 rH CM O o o O O O O O O O o O O O II, II H H H H H H H H H rH H H H r-1 W (H * * * * * * * * * * * * * * CM CO o o _=1- CO CM CO LTN LTN VO LTN CO on on rH CM rH LTn H rH H CM on On _=f CM _=r CO W Eh a M Eh O « CO O o en w — lt= H O CM H On o o o >* 1 H (L> o II ^^ — ^ O "*^*- V >^ > — «* v„ • 1 1 -3- Pi G to en rH CO C a , and X >, 1, where m and n are integers and X is a real, multiple precision number. The resulting P's are multiple precision. For the special cases X = 1 or m = 0, P (x) = 1 and P n (X) = m m for n=l,2, . . . , nmax. P n (X) m For other cases define F = —, t-t . The F then satisfy the n ( n+m ) ! n recurrence relation F A + 2nX F + J2Z2ZL1 F .(, (n-1,2,3, n+1 (n + m + l)^Q- n (n+m+l) a- 1 m identity F n + 2 \ F = K . ; / n m! in n=l 29 Since P(x)=Oifn>m,F =0 for n > m and so the recurrence m n relation is related to the finite continued fraction (n+m) n-1 (n+m)(m+l-n) (n+m+l)(m-n) (n+m+2) (m-n-1 ) 2m» l r v nX x + (n+DX^ (n+2)X x + " 'mX^ u - n ^ where X = 2X/x 2 -l. The algorithm (2.7) can now be applied wi v = m and no iteration is needed since r = , m F m+1 th = 0. The program proceeds as follows: r = o. r = (m+1 7 n) r- m n-1 nX n +(n+m+l)R 1 r S =0,S n =R n (2 + S) m n-1 n-1 n m ? n=m,m-l, ... ,1 J P m (X) = m!F = m! (X+/X^l) m m! 1+S rt 3 n+l m = (m+n+1)! F = (m+n+l) P n R (n=0, 1, . n+1 m n m- -1) If nmax < m the last line above is calculated only for n = 0,1,..., nmax-1. If nmax > m, P is set equal to for n = m+1, m+2, ...,nmax. m Accuracy is affected if X is very close to 1 since the computation of X is subject to cancellation of significant digits. An error re- turn occurs if X < 1 or m < or nmax < 0. 30 ILEG2 - This subroutine generates to d significant digits the associated Legendre functions of the second kind Q (X) for n=0,l, . . . ,nmax given m >, and X > 1, where m and n are integers and X is a real multiple precision number. The resulting Q's are multiple precision. The program first generates Q_(X) from the recurrence relation Q i +1 (x) + - 2±X_ Q i (x) + (i+n) ( i _ n _ l)Q i-l (x ) = o (i=l,2,...,m-l) /FIT n n with n = and the starting values Q„(x) = — In ^rr- and £L(X) = — — — - 2 X-l /%?~± Now let F = Q (X) and apply algorithm (2.7) to the recurrence relation (n-m+1) F n+1 - (2n+l)XF n + (n+m) F^ = (n = 1,2,3,...) 00 with the identity F + ) • F. = Q m (x). O £__, 1 O i=l Step 1: v = 20 + integer part of (5/i+*nmax) inc = 10 = (n=0,l,. . .nmax) Step 2: R (v) = 0, R (v ] = (n+m) , , ^ v n-1 7 — r (n=v,v-l, . . . ,1;, (2n+l) X - (n-m+1) IT Vj n F< v) - <$(x). F { n ll = R< v) F< u) (n-0. 1 imax-1) Step 3: If — 7 — \ > .5*10 for any n = 0, 1, . .., nmax n then (v) = F (v) (n=0,l , . . . ,nmax) , n n v = v + inc , inc = 2*inc, go to Step 2. Otherwise accept F as the final approximation to Q (X) for n = 0,1, . . . ,nmax. Convergence of this algorithm is slow for X near 1. An error return occurs if X v < 1 or nmax < or m < 0. ILEG 3 - This subroutine generates to d significant digits the associated Legendre function of the second kind Q (X) for m=0,l,..., mmax given n >, and X > 1 where m and n are integers and X is a real multiple precision number. The resulting Q's are multiple precision. The program first calls ILEG2 twice with parameters X, m=0 and nmax = n and X, m = 1 and nmax = n to obtain Q (X) and Q (X). The remaining Q's are generated from the recurrence relation Q m+l (x) + 2mX_ Q m (x) + ( m+n ) ( m _ n _ l)c f-l (x) = ( m =i, 2 , . . . ,mmax-l) n fitt n using the initial values Q°(X) and Q (X) obtained from ILEG2. 32 An error return occurs if n < or mmax < or X v < 1. LEG1 - This subroutine generates to d significant digits the associated Legendre functions of the first kind „n, > r(A+n+l) f rv . rprrr +1 A . ,. A (X) = Trr(A+l) J •■ COS - 1 C ° S for n = 0,l,...,nmax given X >, 1 and A where n is an integer and X and A are real multiple precision numbers and A is not an integer. (ILEG1 should be used when A is an integer). The resulting P's are multiple precision. For the special case X = 1, P A °(X) = 1 and P?(X) = for n = 1,2,.... nmax. If A < - - , A can be replaced by -A-l since p"(X) = P n . (X). This substitution is made to avoid loss of accuracy when X is large. P A (X) Let F = —7- -T- and then algorithm (2.7) is applied to the n r(A+n+l) & * recurrence relation F n + 1 + T^ST) X l F n + (n^l) F n-1 " ° ("-4.2,3....) u 2X where X^ = "1 Z / - 1 A F n= ' MA*!)" n=l Step 1: v = 20 + integer part of f[ 37.26+. I 37.26 + 1283(A+38.26)X] 1283(A+l)X ) 33 *nmax inc =10 0^ V) = (n=0,l,. . . ,nmax) Step 2: R ( ^ = 0, R (v ] = -^=£1 n nX n +(n+A+l)R 1 n \ (v) \ n = v, v-1, . . . , 1, S (V) = 0, S (v ? = R (v > (2 + S (V) ) v n-1 n-1 n ,(v) (x+»ftZ-i) J r(A+i) 1 ; & , F ( ^» = R (v) F (v) (n=0,l nmax-1) n+1 n n Step 3: If i F (v) _ 0(v) n n > . 5*10 for any n = 0,1, .. . ,nmax then (v) = F^ (n = 0,1, . . . ,nmax) , n n v =v+inc, inc = 2*inc, go to Step 2. Otherwise P*J(X) .H. T(A+n+l)*F (n=0,l, n . ,nmax) Accuracy is affected if X is very close to 1 since the computation of X is subject to cancellation of significant digits. The rate of convergence of the algorithm decreases as X increases. An error return occurs if X < 1 or nmax < or A is an integer. 3U LEG 2 - This subroutine evaluates to d significant digits the associated Legendre functions of the first kind P., (X) for n = 0,l,...,nmax A+n given X >„1, A, and m >^ where m and n are integers and X and A are real multiple precision numbers and A is not an integer. The resulting P's are multiple precision. The program first calls LEG1 twice with parameters X, A and m and X, A+l and m to obtain P (x) and P (X). The remaining P's are generated from the recurrence relation m , » 2n+2A+l x p m m . n+A+m m , , _ , F A+n+l Uj n+A-m+1 X P A+n U} + n+A-m+1 P A+n-l (X} " ° tn-0,1,..., nmax-l) using the initial values P„(X) and PT (X) obtained from LEG1. An error return occurs if m < or X < 1 or nmax < or A is an integer. CONIC - This subroutine generates to d significant digits Mehler's conical functions P , (X), a special case of the associated Legendre function of the first kind, for n = 0,1,..., nmax given X >, 1 and T where n is an integer and X and T are real multiple precision numbers. The resulting P's are multiple precision numbers. For the special case X = 1, P° l/2+iT (X) = 1 and P^ ±/2 + iT (X) = ftr n= 1,2, . . . ,nmax. (x) 35 —1/2 + iT For other X, let F = ; and apply algorithm (2.7) to the recurrence relation nX x (n - |) + T' F n+1 + (n+l) F n + n ( n+1 ) F n-1 = ° (n = 1 ' 2 "" ) where X. = 2X 1 v6^I and the identity F + Y X F = cosr . T ln(X+v6i^l)] n=l where X = n! n r(f + iT) r(- + iT + n! r(|- it) r(| - iT + n) Step 1: v = 30 + interger part of { [l+( .1U0+. 02U6T) (X-l) ]nmax} inc = 60 ^ = (n=0,l,. . . ,nmax). n Step 2: Obtain x' v ' » JT V ' from v+1 V (v) _ 1 2~ ' A 2 (v) 3-^ (£ + T 2 )(|+ T 2 ) 36 (v) n+l -I 5-2-7^2 (2 V X „-1» ("=2-3. --.v) t 1 + sr» + v n-1 1 2 T 2 ^n n X 1 + (1 + i-)R (v) n' n S (v) = , S (v > = R (v > (A + S (V) ) v n-1 n-1 n n (v) l n-l = 2A - 1 2 T 2 n n+l \ > n=v , v-1, . . . ,1 J v .). co 5 [ T 'in ( X+. £5iL, F (v) =R .5*10 for any n=0,l, . . . ,nmax then (v) = F (v) (n=0,l,...,mnax), n n v = v + inc , inc = 2*inc, go to Step 2. Otherwise p _ l/2+iT (X)= n!F n (n=0,l, . . . ,nmax) 37 This algorithm converges slowly when X and T are both large. It may not converge at all due to accumulation of roundoff errors if d is too large. Therefore the iteration is terminated and an error message is printed if v becomes greater than 1500. An error return occurs if X < 1 or nmax < 0. TORPID - This subroutine generates to d significant digits the toroidal functions of the second kind Q ,_ (x) for n=0,l, . . . ,nmax given X > 1 a multiple precision number and m an integer. The result- ing Q's are multiple precision. Let F = Q m , (X). Algorithm (2.7) is then applied to the recurrence relation (n-m+ |) F n+1 - 2nXF n + (n+m- |) F = (n=l,2,3,... and the identity F n + 2 ) F £_, n=l I (-if y/f T(m+ |)(X-l) m 2 (X+I^ 2 Vx-i / Step 1: v = 20 + integer part of { [l. 15+( . 0lh6+. 00122m)/ (X-l ) ]nmax} inc = 10 = (n=0,l,. . . ,nmax). Step 2: R (v ^ = 0, ,(v) n-1 n+m - (v) = 0, S (v) n-1 2nX-(n-m+ ^-)R (v) 2 n R (V ] (2 + S (v) ) n-1 n n = v, v -1, ... ,1, 38 I— . - 1 fx+i\ ,(v) = (-i) m y| r( m+ | )(x-i) 2 Vx-ij F VW = y-x, y 2 M^ g MA-x . XA -x; , } ( v ) ( v ) ° , . q (v) ' F n + 1 - R n F n 1 + S (n=0,l, . . . ,nmax | F ( ^- j£ v) | _ d Step 3: If 1 — \ > .5*10 for any n=0 ,1, . . . ,nmax |f (v) I 1 n ' then (v) = F (v) (n=0,l» • • • ,nmax) , n n v = v + inc, inc = 2*inc, go to Step 2. Otherwise accept F as the final approximation to .m (X) for n=0,l, . . . ,nmax. - 1/2 + n Convergence is slow for X near 1. An error return occurs for X v < 1 or nmax < 0. 39 FOOTNOTES 1. I.D. Hill, "Algorithm 3^ Procedures for the Basic Arithmetical Operations in Multiple-Length Working," Computer , and x >>1, where m and n are integers and x is real. For this special case where m is a non- negative integer, P n (x) . 1x^1) m ^m . n/2 ,m+n 2 m! dx m+n (x2-l) m m+n Thus if n > m, P n (x) = since - — — (x 2 -l) m = 0. Because ILEG1 takes m n m+n dx advantage of this fact and sets P (x) equal to for n > m, we need only consider the upper right triangular matrix, which contains the non-zero values of P m (x), for our testing purposes. An 11x11 matrix was chosen for study (0 ^ n ,< m ,< 10). p . » . 1 . 10 10 ,10 10 U3 The absolute error of the diagonal elements was determined by comparison with values computed by means of the explicit expression: P*(x) = (x 2 -l) n/2 -jfj- (2n-2i+l) i=l Since the maximum length of the arrays used to hold the long precision numbers in ILEG1 is 50 locations, the P (x) were calculated by a pro- n gram which used arrays of length 52 to eliminate the accumulation of rounding error in the 50th location. Using arrays of length 52 assures that the error in P (x) is at most 1 in the least significant digit of the number, this coming from truncation of the final number. The error in the other elements of the matrix can be related to the error of the diagonal elements by means of (U.l). The following table shows this relationship, where e is the absolute error of P (x), determined from the value of P (x) computed by ILEG1 and by the ex- plicit expression as described above. /x z -l x e l x 2 -l 2x 2 c £^I (x 2 -l) 5 10!x i 5 10 x 2 -l 2x 2 '10 1 '10 '10 kk This table was arrived at in the following manner. Consider for example 7*7 that {P Q } = P Q + e and that all other values are exact. By (U.l) 9 9 we have : Similarly. {P 9 } * ■ 7— [2x{P^}* - 16 pj] {p o } * = — — t 2x ( p o + 0-16 pj; 9 ^2Ii 9 8 8 * 8 Pv { p 9> = P9 + 7=Z * * 1 8 * 8 {Pp = -^- [x{P°} - 17 P„] 9 1^=1 9 8 [r 9 3 * - 1 [x(P® + 2x e) - 17 Pq] 9 ^^ ^cr

>" = P* + &L . * x 2 -l 2x 2 y * m Now we know the value of e g = [{P^>* - P*] so our conclusion is that if (U.l) is satisfied by the values in the matrix, then the error 7 x 2 -l e 9 2 € 9 ' y 2x 2 y The table can also be expressed by the formula: , 2 n \(m-n) /2 n (x -1 J e = £ . . r, (m-n) m' (m-n) !x Since x > 1, (x 2 -l) (m " n)/2 < x m " n and the following inequality holds n 1 £ < m (m-n)! m Then if P . = mm (m-n)! P. ■in _^ ^ we can c alculate the maximum relative error for each column, r^, by £ m n_ = m P . mman Table 3 shows the results of tests run. K is the number of locations in the array used for each multiple precision number and D is the corresponding number of digits which can be stored in an array with K locations. Maximum relative residual refers to the check of recurrence relation (U.l)- It is defined as: CO CO CO t- C— t— CO t— t— C^ CM CM CM CM CM CM CM CM CM CM O o O O O O O O O O H H H H H H H H H H * * * * * * * * * * o o H 0\ vo CM -st H vo 00 oo o VO o\ oo H H CM t- H CM CM \s\ CM II CO t— t— t— t— tr- t— [— r— t- X CM CM CM CM CM CM CM CM CM CM O O o O O O O O o O H H H H H H H H H H 9 * * * * * * * * * * CO H on _st LTN vo CM CM LT\ CM * fl CO H H H H H H H H H CM II CO r— CO vo vo VO o o O H H H F| * * * o ^o H CO H H VO * * • t- r— CO VO vo VO o o o H H H 3 * * * o H -st a H H 0\ CO ON r— t— CO CO [— t- C-- CO CM CM CM CM CM CM CM CM CM CM O O o O o o O O O O H H H H H H H H H H * * * * * * * * * * O O -st o CO CM H CM ltn -st -st -st H H ITv H 00 H H CO On CM oo II H « W n i-l w m ►j < H H O H CO CO r— CO CO CO CO co CO CO CM CM CM CM CM CM CM CM CM CM O o o O o O O O O O H H H r-\ H H H H H H * * * * * * * * * * O VO -st o o\ O t— H CM VO r-\ t— H t- t- ON t— vo o\ t— CO CO CO CO CO t- r— CO CO CM CM CM CM CM CM CM CM CM o o o O O 1 O o 1 O O H H H H H H H H H * * ^ * * * * * * H CM CO t- H CM H O -st -si- CM VO H CM rH CM t— LT\ o o ON VO II Q o CM CO CO CO VO vo vo o o o H H H * * * O r— CO H On vo CO t— CO VO VO VO o o o H H H * * * H -st 00 ON rH CO CO CO CO vo . vo VO o o O H H H * * * H CM VO -st CM ON II CO C— t— CO t— t— t- t— CO r— CM CM CM CM CM CM CM CM CM CM o O O O O 1 O o 1 o 1 O i o H H H H H H H H H H * * * * * * * * * * H H O VO H -Si- CM -St O O vo H H ON H nn H H CO H II t- CO CO vo VO vo a o o o H H H a * * * -St CM t- H ON On CO CO c— C— C— 1— r- h- h- I— CM CM CM CM CM CM CM CM CM CM O O o O O O O O 1 o 1 o H H H H H H H H H H * * * * * * * * * * t— -St H CM OO O ON VO OO O CO o\ CM OO H CM CM OO OO OO CO CO OO VO VO VO o o O H H H * * * vo H ON LTN H vo II t— t— t~- t- vo t- t— r- h- 1— X CM CM CM CM CM CM CM CM CM CM § O o o O 1 1 1 1 o o o o 1 O O H H H H r-\ H H r-\ r-< H a * * * * * * * * * * ^ CM vo H CM H VO 0O -st VO UA -St LTN CM -st H H OO CM -St CM CO t- t- VO VO vo | o o o H H H a * * * i o c— H t- CM CM CMOO-st LP.VO C— CO Os CM 00 CO fr- r— CO fr- fr- r— NO no NO NO NO NO no o o o O O O o H r-\ H H H H H a * * * * * * * cr oo ir\ O CO no H no CO CM CM -3" CM 00 CM CO fr- fr- co t— fr- fr- no NO no NO no no NO m r-{ o o o o o o H r-\ H H H H a * * * * * * * ^ oo -3- CM On J" ir\ NO ON rH H 00 H H H fr- t— CO t — t — t 1 — fr— NO NO NO NO NO NO NO I I I I I I I o o o o o o o £ * * * * * * * cr oo ^t H LTN O UA O H H fr- H CO H 00 hi fr- t— CO r— r— fr- fr- NO NO NO NO NO NO NO O O O O o O O * r-{ H H H H H H CO * * * * * * * T> S CM CM H LP, -d- CM OO ON , and x > 1, where m and n are integers and x is real. It was mentioned in the beginning of this chapter that a program could be tested by independently calculating every value but that generally an explicit expression was known only for a very few values of m and n. In the case of Q (x), the explicit expressions have been generated for 0$n,m^J+ by Suschowk. These formulae were used, with k the exception of Q, which is incorrect and should read 5 + 630x U - - 105x T + 385x 5 - 511x 3 + 2T9x] Qax) = - — [(l05x 8 - l+20x 6 + 630x 1| -l*20x 2 + 105 )*ln[ (x+1 )/(x-l ) ]/2 4 (x 2 -l) 2 U9 by a program which used arrays of length K+2 for storage of the multiple precision numbers so that accumulation of rounding error would not affect the accuracy of the numbers calculated. The r m ->* results, 1Q i , truncated to multiple precision numbers of array length K can then be considered as exact values and the relative error of the Q calculated by ILEG2 can be determined by n {Q m }* - Q m m n n n = {Q m } n Table h gives the maximum relative error for each column generated by ILEG2, max m 1 = « i 1 m 0£n$4 n for various values of x and K. D is the number of significant digits asked for by the calling program. The largest relative -27 error observed was . U7*10 when x = 10.0, K = 10 and D = 28 indicating that when K is so close to the total number of digits which can be stored in the multiple precision number of array length K ((l0-3)*U+l =29), there may only be D-l significant digits in the result. I LEG 3 - ILEG3 generates to d significant digits the associated Legendre function of the second kind Q (x), for m=0,l , . . . ,mmax given n^O and x>l, where m and n are integers and x is real. 50 m 1 2 3 k x=l.l .13*10 .15*10 .88*10 .17*10 -27 -27 -28 -27 TABLE k ILEG2 K=10 D=28 x=5.0 x=10.0 n m .67*10 .78*10 .11*10 .77*10 92*10 -28 -28 -27 -28 -28 r ! m .1+7*10' .1+0*10 .86*10 .10*10 .15*10 ■27 -28 -28 -27 -27 x=25.0 m -27 .1+3*10 ' .15*10" 2T .10*10" 2T .19*10" 2T .1+1*10~ 2T x=l.l K=15 x=5.0 D=l+8 x=10.0 x=25.0 m 1 2 3 1+ n m .66*10 .9^*10 .96*10 .15*10 .26*10 -1+8 -1+8 -1+8 -1+7 -1+7 m .1+6*10 .1+7*10 .1+9*10 .86*10 .Il+*10 -1+8 -1+8 -1+8 -1+8 -1+7 m .12*10 .28*10 .1+9*10 .81*10 .21*10 -1+7 -1+8 -1+8 -1+8 -1+7 m .1+3*10 .77*10 .6i+*io .63*10 .13*10 -1+7 -1+8 -1+8 -1+8 -1+7 K=20 D=68 x=5.0 m 1 2 3 1+ m .33*10 .95*10 .88*10 .23*10 .31*10 -68 -68 -68 -67 -67 51 These are the same values as those generated by ILEG2, the difference being that ILEG2 generates the function values a column at a time and ILEG3 generates them a row at a time. Therefore the relative error of the Q calculated by ILEG3 can be determined exactly as described above for ILEG2. Table 5 gives the maximum relative error for each row generated by ILEG3, _ max m n n " 0^m$l+ n n for various values of x and K. The largest relative error observed —27 here is also .1*7*10 when x=10.0, K=10 and D=28 so the conclusion stated for ILEG2 is also true here: when K is so close to the total number of digits which can be stored in the multiple precision number of length K, there may only be D-l significant digits in the result . LEG1 and LEG2 - LEG1 generates to d significant digits the associated Legendre function of the first kind P (x), for m=0,l,..., mmax given x >, 1 and a where m is an integer and x and a are real, a not an integer. LEG2 generates to d significant digits the associated Legendre functions of the first kind P (x), for n=0 ,1, . . . ,nmax given m >, 0, x >. 1 and a where m and n are integers and x and a are real, a not an integer. 52 x=l.l n_ n -27 .12*10 ' • ll|*10" 2T .13*10" 2T .13*10" 2T .85*10" 28 x=l.l .19*10' .15*10 .18*10 .iU*io .1U*10 •hi -hi -hi -hi -hi TABLE 5 ILEG3 K=10 D=28 x=5.0 x=l0.0 n PR .92*10 .ii*io" 2T % -27 .35*10 ' .1+0*10~ 2T . 61+ *io~ 28 .1+5*10~ 2? .78*10" 28 .U7*10" 2T .67*10~ 28 .1+7*10" 2T K=15 D=l+8 x=5.0 x=10.0 n n .13*10 .11*10 .67*10 .1+2*10 .1+7*10 -hi -hi -1+8 -1+8 -1+8 .12*10 .lU*10 .12*10 .7^*10 .23*10 -hi -hi -hi -1+8 -1+8 x=25.0 n -27 .1+0*10 ' .36*10~ 2T .39*10" 2T .1+0*10~ 2T .1+3*10" 2T x=25.0 n n .1+0*10 .38*10 .37*10 .39*10 .1+3*10 -hi -hi -hi -hi -hi K=20 x=5-0 n n .31*10 .17*10 .ll+*10 .13*10 .12*10 D=68 -67 -67 -67 -67 -67 53 For these cases of P (x) where a is not an integer or restricted to 9. a half-integer an independent calculation is very difficult since there are no explicit expressions available and the infinite series representation is extremely complicated and slow to converge. It is for this very reason that Gautschi's algorithm is so useful since it does not require that any values of the function be known in advance. In fact this author found no previous table computations made for a other than integer or half-integer values. It does however make the testing of these pro- 2 grams difficult. Zhurina and Karmazina have recorded, as part of their work on the conical functions, the formulae for P. ,_ . (x) and P _ /n , . (x), -1/2+1T -1/2+1T valid for |l-x|<2. These were used with t = in a program which stored the multiple precision numbers in arrays of length K+2. The results, {P , p (x)} and {P , (x)} , truncated to multiple pre- cision numbers of length K were then considered exact. Test 1 of Table 6 shows the relative errors ni = {P m -1/2 } * - {P* 1 I 1 1 1 -1/2 1 m * {P -l/2 } m=0 or 1 and n 2 {P m -1/2 i 1 -1/2 1 (P m Y 1 -1/2* m=0 or 1 of {P^n/p} 1 and {P 111 , } 2 generated by LEG1 and LEG2 respectively. In the case of LEG2 which uses LEG1 to obtain values of P and P ^, , the a a+1 5h calling program specified a = -2.5 and n = 2 and then tested the relative error of P so that we would not merely be retesting values generated by LEG1 as would be the case if a = -.5 and n=0 was specified. -26 The maximum relative error observed was .3^*10 ' for LEG2 with x=2.5» K=10 and D=25- This indicates that there were 25 significant digits generated as specified by the calling program. To check other values of m and a, an "internal consistency" test was run. First a 5x5 matrix was generated from both LEG1 and LEG2, a o (P 1 } l a o {P^ I 1 ... a o 'V 1 (P° } ± 1 V u <- Ji where i=l or 2 for values of LEG1 or LEG2 respectively, and then the relative difference at each point was determined m n = a {p m } l _ {p m } 2 a a {P m } X a , 0$m$l+, a Q x 1 and x where m is an integer and x and x are real. * 1 * Exact values, {P (x)} and {P n ,_. . (x)} , were obtained -1/2+ix -1/2+ix from the formulae of Zhurina and Karmazina, as described above under LEG1 and LEG2, for testing values of P _ ,_ . (x) and -1/2+ix P _ ._, . (x) generated by CONIC. Test 1 of Table 7 shows the -1/2+ix ° J maximum relative error observed for each combination of x and x, n max max 0 1 and m where m and n are integers and x is real. The testing of T0R0ID follows the same pattern as that of ILEG1 : determining the error of key elements in a matrix of values generated by T0R0ID and then finding an approximation to the errors of the other elements by means of the recurrence relation (U-l). A 5x6 matrix was chosen for study (0 ,+1/2+n (x2-l) lA (xt^-l) n 00 E k=0 J[(^ + i)(|-^i) k!(n+k)! (-t) J where t = 59 x - /x 2 -l 2/x z -l *5 m derived from a formula for Q n (x) valid for |x| > 1, was used in a program which stored the multiple precision numbers in arrays of k * length K+2. The results, {Q -, i . } 9 n=0, . . . ,4, truncated to multiple precision numbers of length K were then considered exact and used to determine the absolute error "" k > k of the Q , generated by TOROID. The relationship of the errors in the other elements to the errors of the elements in the 5th row, arrived at as explained in the section about ILEG1 , is given by the following table. 2JL (x 2 -!) 2 k i_ (x 2 -!) 2 k 2 k (x 2 -!) 2 k h , 2 n ,2 , 5-3 h e i 5-3 k e l 3^7 h~ e 3 '" - (x T l} e k 2 p x 2 9-7-5-3~5 e 9 ? 2 3 x (x 2 -l)2 U £ ( x 2 -l)? l+ 2 3 (x 2 -l ) ? U ... £ (x 2 -l)2 k 5'3 2 3 £ 1 5-3 3 £ 1 3 3 £ 3 7^3 3 £ 9 2 X 2 X 2 X 2 £x(izD h i_ (x 2 -l) k 2 2 (x 2 -l) U ... 2 2 (x 2 -l) U 5* x 2 | 5-3 x 2 £ | 3 x 2 £ | F3 — £ | 6 ^^l f ^-J 2(Al)l .* 2 (x 2 -l)l k 7 | 5 X | 3 x a — J — e 2 2 - i 2 e Q 222 I TABLE 8 TOROID 60 m 1 2 3 4 x=1.5 m .i6*io~ 25 .15*10~ 25 .59*10" .6o*io~ 2T .47*io~ 28 r max .3U*10" 2T x=1.5 K=10 x=3.0 \ .8l*l0~ 25 .6i*io~ 25 .21*10" 25 .20*10" 26 .15*10" 2T r max .22*10~ 2T K=15 x=3.0 D=25 D=i+5 x=10.0 m ,97*10 ,70*10 ,23*10 ,32*10 ,32*10 -25 -25 -25 -26 -27 max 29*10 x=10.0 -27 m 1 2 3 1+ m m m 82*10 75*10 30*10 30*10 49*10 -k6 -k6 -h6 -hi -48 ,41*10 ,30*10 ,11*10 98*10 ,92*10 -46 -45 -45 -47 -48 28*10 18*10 60*10 11*10 80*10 -45 -45 -46 -46 -48 max 37*10 -47 max 18*10 -47 max 48*10 -47 m 1 2 3 4 K=20 x=3.0 n. D=65 m ,81*10 61*10 21*10 20*10 15*10 -65 -65 -65 -66 -61 rri.'j./ 25*10 -67 61 Table 8 gives the maximum relative error for each row generated by TOROID m max n m " 0 1,2,3 3 All) = IL = -IL 2 All) = A(l) + 1 IFU.EQ.O) RETURN 1 P = IL J = 3 5 P = P/T IF (P ) 6,6, 7 7 J = J + 1 PROCEDURE SFT IN THAT X IS REAL. IT IS IMPORTANT TO REALIZE THAT THE PROCEDURE IS THE LONG EQUIVALENT OF A = X . IT DIFFERS FROM KEAL X INTEGER All) SUBROUTINE LNGTHN ( X , A, K ) END RETURN 17 4121 = J - 3 15 A (P) = 16 no IS P = I|_,K IF ( IL - K ) 16,16,17 1 1 IL = J + 1 (.0 TO 10 P = P - 1 IL = Q A ( P) = IL - T*G 9 0= IL/T 10 IF! P-2 ) 11,11,9 P = J h I F( J.GT.K) RFTURN 1 ',0 TO 5 RESULTING ACCURACY CANNOT BE ANY GREATER THAN THAT OF THE OR- DINARY REAL REPRESENTATION. COMMnN/LONG/M,w,T ,S,R,TR,FX,RM 79 10 15 16 20 25 30 35 INTEGER N,W,T,S,EX REAL R» TRtRM REAL Y,Z,XL,XT INTEGER J»M,V M = -1 XL = -X IF(XL) 3t2,l M » XL = -XL M = M + 1 A(l) = M IF(M) 5»35,5 V = 1.0 J = IFIXL - TR)7,8,8 XL = R*XL J = J ♦ 1 GO TO 6 IF(XL - 1.0) 9, 10,10 XL = TR*XL J = J - 1 GO TO 7 A(2 ) = J V = W/N + 3 I F ( V.GT ,K ) V = K DO 15 M = 3tV J = XL AIM) = J XT = J XL = T*( XL - XT) M = V + 1 M = M - 1 I F (M.EQ.3.0R. A(M) , AIM) = A(M-1 ) = A(M-1 ) + GO TO 16 IF (A(3 ) .LT.T ) GO TO 25 A ( 3 ) = 1 A(2 ) = A(2 ) + 1 L = V + 1 00 30 M = L,K AIM) = RETURN END SUBROUTINE SHOR TN ( A , X , * , K ) INTEGER All) THE RESULT IS THE ORDINARY RFAL R EPRESENT AT I ONOF HELD IN THE ARRAY A. RETURN TO THE LABEL PARAMETER FLOW WOULD BE PRODUCED. COMMON/ LONG /N,W,T,StRtTR,FX,RM REAL R, TR,RM INTEGER N,W,T,S,EX REAL X,Y INTEGER J,V IF (N*( A(2 l-t-1 ) .GE.EX ) RETURN 1 CALL UNDFRZ ( 'OFF • ) V = W/N + 2 LT.T) GO TO 20 1 THE NUMBER OCCURS IF OVER 10 20 K - 1 AS IF IT COULD BE TIDIED UP SOME 80 X = A(3) IF(V.GT.K - 1) IF(V.LT,4) V=4 DO 10 J = 4 t V X = Y*A(J) + X Y * R*Y THIS LOOKS J = V + 1 Tl = IFU.LT.K .AND. A(J+1).GE.S) Tl * 1 X = Y*( A(J ) + Tl) + X Tl = TR**A(2) IF(T1*RM*X.GT.1.0)G0 TO 20 X = T1*X Tl = All) X=T1*X CALL UNDERZ( • ON • ) RETURN CALL UNDERZ( 'ON' ) RETURN 1 END SUBROUTINE AD JUST ( D ,H ,K ) INTEGER 0<1), H THIS PROCEDURE DEST ANDARDI Z ES A LONG NUMBER READY FOR ADDING OR SUBTRACTING. IT IS NOT INTENDED THAT IT SHOULD BE CALLED EXCEPT BY PROCEDURE ADD OR PROCEDURESBTRCT . H IS THE DIFFERENCE IN EXPONENTSt I.E. THE NUMBER OF PLACES D IS TO BE SHIFTED RIGHT. COMMCN/LONG/N,W,T, S,R,TR,EX,RM INTEGER REAL R, INTEGER REAL RT HL = H D(2) = D(2 ) HL = K - HL 1 F (HL - 2 ) 1 0(1) = RETURN I F(D(3) - 1(1) = RETURN D(K ) = 1 L = K - 1 I F (L-3 ) 10,6,6 DO 8 J = 3,L D( J ) = RETURN RT = 2*D(HL + 1 ) T 1 = R*kT D(K ) = D(HL ) + HL = HL - 1 L = K - 3 !0 Tl = 1,L rHIS COULD BE J = V - T 1 N,W,T,S,EX TR,RM J,HL,L,T1 T2 HL 2,3 S) 4,5,5 Tl T IDIFD UP IFIHL - 2) 15,15,16 8l 16 D( J) = D(HL) HL * HL - 1 GO TO 20 15 D(J) » 20 CONTINUE 10 RETURN END SUBROUTINE RE AD JT ( D, E ,K ) INTEGER E,0(1) C THIS PROCEDURE READJUSTS A LONG NUMBER IN CERTAIN CASES WHERE C ROUNDING CAUSES A FORM OF OVERFLOW. IT IS NOT INTENDED THAT IT C SHOULD BE CALLED EXCEPT BY PROCEDURE ADD OR PROCEDURE MLTINT. COMMON/LONG/N,W,T,S,R,TR,EX,RM INTEGER N,W,T,S,EX REAL R, TR,RM INTEGER J,EL,L,T1 REAL DT EL = E 0(2 ) = D( 2) + 1 DT = 2*D(K) D(K) = R*OT D(K ) = D(K ) + D(K-1 ) L = K - 4 00 5 Tl = l,L J = K - Tl 5 D( J ) = D( J-l ) D( 3) = EL J = K + 1 10 J = J - 1 IF INTEGER 4P< 1 ) ,BP( I) , C( 1 ) THIS PROCEDURE IS THE LONG EQUIVALENT OF C = A + B COMMON/ LONG/N,W,T, S,R,TR,EX,RM INTEGER N,W,T,S,EX REAL R, TR,RM INTEGER A( 50) ,B( 50) C THESE DIMENSIONS LIMIT K TO A MAXIMUM OF 50. INTEGER J,E,T1 CALL COPY(AP,A,K) CALL COPY(BP,B,K) IF(B(1 ) )2, 1,2 1 CALL COPY(A,C,K) RETURN 2 I F( A< 1 ) )4,3,4 3 CALL C0PY(8,C,K ) RETURN 4 I F (A ( 1 ) .EQ.BI 1 ) ) GO TO 50 10 50 15 25 30 35 .LT. T) GO TO 35 - T CALL READJT(C,E,K) IFU(l)) 10,10,5 B(l ) = 1 CALL SBTRCT(A,B,C,K) B(l) = -1 RETURN A(l) = 1 CALL SBTRCT(B,A,C,K> A(l) = -1 RETURN C<1> = A(l) IF(A<2) - B(2) ) 25*30,15 CALL ADJUST(B,A(2)-B(2), K) IF(B(1M 30,1,30 CALL ADJUST(A,B(2) - A(2),K) IF(A(D) 30,3,30 C<2) = A(2) E = 00 35 Tl = 3,K J = K + 3 - Tl C< J ) = A( J ) + B< J) + E E = IF(C(J) E = 1 C(J ) = C(J ) CONTINUE IFtE.EO.l ) RETURN END SUBROUTINE SBTRCT ( AP , BP ,C ,K ) INTEGER AP( 1),BP( 1) , C( 1) THE PROCEDURE IS THE LONG EQUIVALENT COMMON/LONG/N,W,T,S,R,TR,EX,RM INTEGER K,N,W,T,S,EX REAL R, TR,RM INTEGER J,E,JJ,M INTEGER A(50),B(50) THIS DIMENSION LIMITS K TO A MAXIMUM OF CALL COPY(AP,A,K) CALL COPY(BP,B,K) I F ( B ( 1 ) ) 1,2,1 CALL COPY(A,C,K) LABEL 2 IS ALGOL RETURN IF(A ( 1 ) ) 4,3,4 B( 1 ) = -B( 1 ) LABEL 3 IS ALGOL CALL COPY(BrCK) RETURN I F( A( 1 ) - B( 1 ) ) JJ = B ( 1 ) B ( I ) = All) CALL ADD(A T B,C,K) B< 1 ) = JJ RETURN C ( 1 J = A ( 1 ) 1 F (A i? ) - R( 2) ) 7,8,9 CALL ADJUST(B,A( 2) - B(2), K) 82 OF C = A - B. 50. LABEL EE. LABEL FF 5,6,5 11 12 20 14 18 16 17 23 22 21 24 26 25 30 31 I F ( S C 1 ) ) 8*2,8 CALL ADJUSTUtBl 2) - A(2>,K) 83 IF (A( 1 ) ) 8,3,8 C<2) = A(2) E = DO 20 JJ = 3tK J = K + 3 - JJ CIJI = AIJI - BIJ) + E IF(CU)) 11*12*12 E = -1 C(J) = C(J) ♦ T GO TO 20 E = CONTINUE IF(E) 14,21*21 C ( 1 ) = -C ( 1 ) E = K IF(C(E)) 17,16*17 E = E - 1 GO TO 18 C (E) = T - C(E) J = T - 1 E = F - 1 IF( E - 3) 21*22*22 C(E) = J - C(E) GO TO 23 E = J = 2 J = J + 1 I F(C( J ) ) 25*26*25 E = E + 1 IF(J-K) 24,28,24 C( 1 ) = RETURN IF(E.LE.O) RETURN M = K - E DO 30 J = 3,M I = J + E C ( J ) = C( 1 ) M = M + 1 DO 31 J = M,K C( J ) =0 C(2 ) = C<2 ) - E RETURN END SUBROUTINE MLTPL Y ( A , B ,C *K ) INTEGER A( 1 ) *B( 1) ,C( 1 ) THE PROCEDURE IS THE LONG EQUIVALENT OF C = A*B. COMMON / L ONG/N * W , T , S * R , TR * E X * R M INTEGER N,W,T,S,EX REAL R, TR,RM INTEGER G( 100) , I ,J,M,H,L THE MAXIMUM VALUE OF K IS DETERMINED BY THE DIMENSION OF G. WITH THE CURRENT VALUE, K MUST BE NO GREATER THAN 50. IF(A( 1 ) .NE.0.AND.B( 1) .NE.O) GO TO 1 C( 1 ) =0 RETURN 8U 10 11 12 IS 16 19 13 20 25 C(l) = A( 1)*B( 1) C(2) = A(2) + B(2) L = 2*K - 1 00 2 I « 3,L G( I ) » DO 10 I * 3,K 00 10 J = 3,K M - I ♦ J - 1 G(M) = G(M) ♦ A( I )*BU) M = M + 1 M = M - 1 IF(G(M) - T) 10,4,4 H = G(M)/T G(M) = G(M) - T*H G(M-1 ) = G(M-l) + H GO TO 3 CONTINUE J = K + 1 IF(G(4> .EQ.O) J = J + 1 IFU.EQ.3 .OR. K.E0.3) GO TO 11 IF(G( J+l) .GE.S) G( J) = G( J) + 1 J = J + 1 J = J - 1 IF(G( J ) - T) 16,15, 15 G( J ) = G( J) - T G( J-l ) = G( J-l ) + 1 GO TO 12 IF(G(4) ) 18, 19, 18 1 = 5 GO TO 20 C(2 ) = C( 2) + 1 I = 4 DO 25 J = 3, K C( J) = G( I ) 1 = 1 + 1 RETURN END SUBROUTINE MLT I NT ( A, I , B ♦ *,K ) INTEGER A(1),B(1),I THE PROCEDURE IS THF LONG EQUIVALENT OF B = A*I. IT FROM MLTPLY IN THAT I IS AN INTEGER, AND NOT AN INTEGER IS MUCH FASTER THAN MLTPLY AND SHOULD THEREFORE ENCc IF POSSIHLE. THE ABSOLUTE VALUE OF I MUST A JUMP TO THF LABEL ARGUMENT WILL OCCUR IF THIS LATED. COMMON/ LONG/N,W,T,S,R,TR,EX,RM INTEGER N,W,T,S,EX REAL R, TR,RM INTEGER J,M,E,IL,JJ IL = I ABS( I ) I F ( IL .GE.T ) RETURN 1 IF( I ) 1,2,3 B(l) = / J = 2,K • ) = ■ kN H( 1 ) r -All) DIFFERS ARRAY. IT BE USED IN PREFER- BE LESS THAN 10**N CONDITION IS VIO- 10 IS 1 1 GO TCI A 3 Rill = Alll 85 U I F ( B ( 1 ) .EQ.O) GO TO 5 F = B<2) = A(2) DO 10 JJ = 3, K J = K + 3 - JJ M = IL*A(J ) + E F = M/T 10 BUI = M - T*E IF(E.NE.O) CALL RF AOJT ( B , F ,K ) RFTURN FNO SUBROUTINE RFCIP( A, RP.K ) INTEGER A( 1 ) , RP ( 1 ) THF PROCFDURF IS THF LONG EQUIVALENT OF B = l./A. THF ORDINARY DIVISION OPERATION IS FIRST USFD UPON THE SHORTFNFO RF PRF SFNT AT I ON OF A TO GET A FIRST APPROXIMATION. SUCCESSIVE APPROXIMATIONS ARF THEN GENERATED RY THF ITERATIVE PROCESS B = R*(2.0 - R*A). IF A=0 THE ORDINARY FAILURE FOR A DIVISION, BY ZFRO OF THE SYSTEM IN USE WILL OCCUR. C.'iMMON/LONG/N,W,T,S,R,TR,EXtRM INTEGER N , W , T , S , E X KEAL R, TR,RM INTEGER C(SO) ,D( 50) ,F(50> ,1 t J»JJ?B( 50) THF DIMENSIONS HERE LIMIT THF PROCEDURE TO K NOT GREATER THAN 50. 1 =" A(2) A ( 7 ) = CALL SHORTNf A, X,£2,K) X=l ,0/X CALL LNGTHN< X,B,K ) I A( 7 ) A (2 ) R(2 ) 0(1) 0(2) D( 3 ) CALL DO 1 D ( J ) i,i i r 5 IS CALL CALL CALL CALL CALL DO 1 = I = B( 2 ) - I = 1 = = 2 COPY(R,F,K) J = 4,K = 5 ALGOL HH. 4 IS ALGOL GG COPY (B,F,K ) COPY(C»R,K) *LTPLY( A,R,C,K) SRTRCT /(X+1), AND REM = 2/(2NU+l)*RATIO**(2NU+l)/(l-RATIO**2). FOR OTHER X, THE RANGE IS REDUCED USING VALUES OF LOGIO AND LOGT COMPUTED ON THE FIRST ENTRY COMMON/ lONG/N,W,T,S,R,TR, EX, RM REAL R,TR,RM INTEGER N,W,T,S,EX INTEGER X(1),Y(1) INTEGER SUM (50) , LOG 10 ( 50 ) t LOGT < 50 ) , TERM ( 50 ) , RATSOR ( 50 ) ,RAT 10 ( 50 ) THESE DIMENSIONS LIMIT K TO 50. INTEGER X2, J,ED,SW,XEXP OAT A LOGIO/ 1,0, 2 ♦3025,8509,2994, 0456, 840 It 7991,4546*8436 ,4207, 1 6011,0148,8628,7729,760 3,332 7,9009,675 7,2609,6773, 2 5248,0235,99 72,0 508,9598,2983,4196,77 84,0422,8624, 3 863 3,4095,2546,5082,8067,5666,628 7,3690,9878,1689, 4 4829,0720,8 3 25,5 546,8084,3 799,8948,26 23/ DATA LOGT / 1, 0, 9,2103,4037,1976,1827,3607,1965,8187,3745, 1 6830,4044, 595,4515, 919, 413,3311,6038,7029, 438,7094, 992, 943, 2 9888,2035,8393,1933,6787,1136,1691,4499,4533,6381, 186, 331,2270, 3 2 666,5149,476 3,9512,675 7,9316,2883,3302,2187,2 337,5199,5793, 49 2/ HERE TO COMPUTE LOG(X). IF{ X( 1 ).LE.O) RETURN 1 XEXP = X(2 ) X(2 ) .= CALL MLTINT(L0GT,XEXP,SUM,£11 ,K ) 11 J = -1 ED = X(3) FIND NO OF POWERS OF 10 IN INTEGER PART OF X. 12 ED = ED/10 J = J + 1 IF(ED) 13,13,12 13 IF(X(3) .GT.3*10**J> J = J + 1 IF (J ) 17,17,14 17 CALL COPY( X,TERM,K) GO TO 21 14 CALL MLTINT( LOGIO, J, TERM, £15, K ) 15 CALL ADD(SUM, TERM, SUM, K ) CALL SET(TERM,10**J,£16,K ) 16 CALL 01 VIDE(X, TERM, TERM, K ) 21 X ( 2 ) = XEXP IS TERM = 1 IF(TERM(3) - 1) 20,18,20 18 00 19 J = 4,K IF(TERM( J) )20,19,20 19 CONTINUE GO TO 170 20 CALL SET (RATSOR, 1, £28, K ) 28 CALL StJTRCT( TERM, RATSOR, RATIO, K ) CALL ADDITERM, RATSOR, RATSOR ,K ) CALL DIVI0E(RATI0,RATSQR ,RATIO,K ) 50 CALL MLTPLY(RATIO, RATIO, RATSOR, K) CALL MLTINT(RATI0,2,RATI0,£51 ,K ) 51 CALL ADD(RATIO,SUM, SUM,K) THIS LOOP SUMS THE SERIES FOR THE SCALED RATIO. 9h DO 60 J = 1,1500 CALL MLTPLY(RATIO,RATSQR,RATIO,K) RATIO WILL NOW BE 2*( (X-l )/ < X + l ) )**< 2J + 1 ) CALL SET(TERM,2*J+1,652,K) 52 CALL DIVIDE(RATIO,TERM,TERM,K) IF(TERMQ).EQ.O) GO TO 170 IF(SUM(2 ) - TERM(2) -K ♦ 2) 53,54,54 53 CALL ADD(SUM, TERM, SUM, K) 60 CONTINUE WRITE<6,999) 999 F0RMAT(30H NO CONVERGENCE FOR LOGARITHM ) CALL 0UTPUT(X,6,K) WRITE(6,998) 998 F0RMAT(15H+ X = ) 54 CALL SET(RATIO, 1,S55,K) 55 CALL SBTRCT(RATIO,RATSQR, RATIO, K) IF(RATIOd)) 56,57,56 56 CALL DIVIDE(TERM, RATIO, TERM, K) CALL ADD(SUM,TERM,SUM,K) 57 CONTINUE 170 CALL COPY(SUM,Y,K) RETURN END 95 ENTRY NAME OUTPUT IDENTIFICATION *** This subroutine is used internally by LSQRT and LNGLOG to print error messages. OTHER SUBROUTINES USED *** IBCOM 96 SUBROUTINE OUTPUT < A, UNI T,K ) INTEGER A(l),UNITtIT,J C THIS SUBROUTINE OUTPUTS THE LONG NUMBER. A, INTHE FOLLOWING C FORMAT,. 15 SPACES FOR DESCRIPTION, (DECIMAL EXPONENT), INTEGER C PART WITH SIGN, DECIMAL POINT, FRACTIONAL PART IN 5-DIGIT BLOCKS, C SEPARATED BY SPACES, AND WITH INTERNAL SPACES SUPPRESSED. COMMON/LONG/N,W,T,S,R,TR,EX,RM INTEGER NtW,T,S,EX REAL R, TR,RM 999 FORMAT ( 15 X,1H( , 16 , 1H ) ,2X , 1 4 ,1H . , 10 ( I 5 , IX ) , ( /24X , 1 1 ( I 5 , IX ) ) ) IT = A(1)*A(3) WRITE(UNIT,999)A(2) ,IT,(A(J) , J ■ 4,K) RETURN END 97 ENTRY NAME EXPO IDENTIFICATION *** Multiple Precision Exponentiation with Integer Exponent PURPOSE *** Given a multiple precision number X of arbitrary array length K(see description of Multiple Precision Basic Arithmetic Package for storage details) and an ordinary integer I, EXPO will evaluate X**I by successive multiplications. OTHER SUBROUTINES USED *** IBCOM, SET, SBTRCT , COPY, MLTPLY and RECIP. USAGE *** INTEGER X(KX) ,Y(KX) ,1 ,K CALL LNGCNS(K,l+,8,75) CALL EXPO(X,I,Y,K) where the parameters in the calling sequence are X - The multiple precision argument. I - The integer exponent . input input Y - On output Y=X**I where Y is a multiple precision number. output K - The array length of the multiple precision numbers X and Y. K.LE.KX.LE.50 input 98 ERROR MESSAGE *** If X.EQ.I.EQ.O control returns immediately to the calling program after the following message has been printed: *** ERROR IN EXPO (0**0) *** ACCURACY *** The maximum relative error observed in tests run was .32*10**(-((K-3)*^-D). EXECUTION TIME *** Using the FORTRAN G compiler the execution time in seconds is given approximately by .69*10**(-U)*K**2*I. This is a good approximation for small I, but when I is very large it over estimates the time required. Note that EXPO requires much less execution time than LEXPO and should be used whenever possible. 99 SUBROUTINE EXPO( XX , A , Y ,K ) C C *** Y = XX ** A WHERE A IS AN INTEGER t K.LE.50 C COMMON /LONG/LN,LW,LT,LS,LR»LTR,LEX»LRM INTEGER LN,LW,LT,LS,LEX REAL LR»LTR»LRM INTEGER XX(l),ZfY(l)fKtSWfNtA ,X(50) CALL SET(X f l,65tK) 5 CALL SBTRCT(XX,X,X,K) IF (X( 1) .EO.O)GO TO 40 CALL COPY(XX,X,K) IF (XX( 1 ) .FQ.O)GO TO 80 SW = CALL SET(Y,l,£10tK) C c *** X**0 = 1 10 IF(A.EO.O)RETURN IF (A.GT.O)GO TO 20 SW = 1 A = -1 * A 20 DO 30 N=1,A CALL MLTPLY( Y,X,Y,K) 30 CONTINUE IF (SW.EO.O)RETURN CALL RECIP( Y,Y,K) A = -1 * A RETURN C C *=** 1**A = 1 40 CALL SET( Y,l f fi50tK) 50 RETURN C C *** 0**0 IS AN ERROR 70 WRITE (6,1000) 1000 FORMAT ('- *** ERROR IN EXPO (0**0) *** ') RETURN 80 IF (A.EQ.O)GO TO 70 C C *** o**A = CALL SET( Y,0,690,K) 90 RETURN END 100 ENTRY NAME LEXPO IDENTIFICATION *** Multiple Precision Exponentiation PURPOSE *** Given two multiple precision numbers X and Y of arbitrary array length K(see description of Multiple Precision Basic Arithmetic Package for storage details), LEXPO evaluates X**Y by the equation exp(Y * In X) using the library subroutines EXP and LNGLOG. USAGE *** INTEGER X(KX) ,Y(KX) ,Z(KX) s K CALL LNGCNS(K,U,8,T5) CALL LEXP0(X,Y,Z,K) where the parameters in the calling sequence are X - The multiple precision argument. Y - The multiple precision exponent. input input Z - On output Z=X**Y where Z is a multiple precision number. output K - The array length of the multiple precision numbers X, Y and Z. K.LE.KX.LE.50 input SUBROUTINES USED *** IBCOM, SET, SBTRCT, SHORTN, EXPO, LNGLOG, MLTPLY a (P. 101 ERROR MESSAGES *** Control is returned immediately to the calling program if either of the following messages is printed: *** ERROR IN LEXPO ( 0**0 ) *** *** ERROR IN LEXPO (X**A WHERE X IS NEGATIVE AND A IS NOT AN INTEGER) *** ACCURACY *** The maximum relative error observed in tests run was .2U*10**(-((K-3)*U-2)). EXECUTION TIME *** Using a FORTRAN G compiler the execution time in seconds is given approximately by .0032*K**3. SUBROUTINE LEXPO ( X , A , Y ,K ) X ** A , K.LE.50 102 C *** Y C COMMON /LONG/LN,LW,LT,LS,lR»LTR,LEX,LRM INTEGER LNtLWtLTtLStlEX REAL LR»LTRtLRM INTEGER X(1),A(1),Y(1),T(50),K,AA REAL AS IF (All ).EO.O)GO TO 20 IF (X(D.EO.O)GO TO 50 IF (A(2).LT.0.OR.A(2).GT.2.OR.A(2).E0.2.ANO.A(3).GT.20)GO TO * CALL SET(T,1,G4,K) 4 CALL SBTRCT(X,T,T,K) IF (T( 1>.EO.O)GO TO 70 CALL SH0RTN( A,AS,£5,K) 5 AA = AINT(AS) IF (AS.NE.AA)GO TO 6 C C *** SPECIAL CASE WHERE A IS AN INTEGER CALL EXPOI X,AA,Y,K> RETURN 6 IF ( X( 1 ) ,GT.O)GO TO 9 C C *** ERROR RETURN WHEN X IS NEGATIVE AND A IS NOT AN INTEGER WRITE (6,1001) 1001 FORMAT (•- *** ERROR IN LEXPO (X**A WHERE X IS NEGATIVE AND A IS N 10T AN INTEGER) ***• ) RETURN C C *** x**A = EXP(A * LN(X)) 9 CALL LNGL0G(X,T,&10,K) 10 CALL MLTPLY(T,A,Y,K) CALL LFXP(Y,Y,K) RETURN 20 IF (X( 1 ).EO.O)GO TO 40 C C *** X**0 = 1 CALL SET( Y,1,£30,K) 30 RETURN C C *#* 0**0 IS AN ERROR 40 WRITE (6,1000) 1000 FORMAT ('- *** ERROR IN LFXPO ( 0**0 ) ***') RETURN C C *** 0**A = 50 CALL SET( Y,0, &60,K) 60 RETURN C C *** 1**A = 1 7 CALL SET(Y,1,£80,K) 80 RETURN END 103 ENTRY NAME LEXP IDENTIFICATION *** Multiple Precision Exponential Function PURPOSE *** Given a multiple precision number X of arbitrary array length K(see description of Multiple Precision Basic Arithmetic Package for storage details), LEXP uses a continued fraction to evaluate e**X. OTHER SUBROUTINES USED *** IBCOM, COPY, SET, SHORTN, SBTRCT , MLTPLY , MLTINT, DIVIDE, ADD, RECIP and EXPO. USAGE *** INTEGER X(KX) ,Y(KX) ,K CALL LNGCNS(K,U,8,T5) CALL LEXP(X,Y,K) where the parameters in the calling sequence are: X - The multiple precision argument. ABS(X) .LE. 20*10**8 input Y - On output Y=e**X where Y is a multiple precision number. output K - The array length of the multiple precision numbers X and Y. K.LE.KX.LE.50 input ERROR MESSAGES *** If ABS(X) . GT. 20*10**8 control is returned immediately to the calling program after the following message has been printed: *** X IS OUT OF RANGE WHEN CALCULATING E**X IN LEXP *** If the continued fraction does not converge in 1000 steps control is returned to the calling program after the following message has been printed: *** NO CONVERGENCE FOR EXP *** ACCURACY *** The maximum relative error observed in tests run was .61**10»*(-((K-3)*U-1)). EXECUTION TIME *** Using the FORTRAN G compiler the execution time in seconds is given approximately by . 0020*K**3. 105 SUBROUTINE LEXP(X,W,K) C C *»* W ■ E ** X FOR ABS(X) .LE. 20*10**8 ♦ K.LE.50 C COMMON /LONG/LN,LW,LT,LS,LR,LTR,LEX,LRM INTEGER LN,LW,LT,LS,LEX REAL LR,LTR,LRM INTEGER X(1),W( 1 ) ,U( 50) , V< 50 ) , SW, A < 50) ,B1 < 50) , B2 < 50) INTEGER ONE(50)tAA(50) ,SW2 ,E1 ( 50 ) ,HALF < 50) ,SW3 REAL** XS DATA El/ It 0» 2 t 71 82*8182 ,8459,0452,3536,0287,4713, 5266 t 2497,7572* 1 4709,3699,9595,7496,6967,6277,2407,6630,35 35,4759,4571,3821,7852, 2 5166,42 74,2 746,6391,9320,0305,9921,8174,1359,6629,0435,7290,0334, 3 2952,6059,5630,7381,32 32,8627,9434,9076,32 33,8298,8075/ C C *** CHECK IF X IS IN CORRECT RANGE IF (X(2) .GT.2.0R.X(2) .EO .2 . AND.X ( 3 ) .GT .20 ) GO TO 160 C C *** E ** X IS CALCULATED FROM A CONTINUED FRACTION FOR O.LT .X .LE .0.5 C FOR OTHER VALUES, X IS REDUCED USING A PRECOMPILED VALUE OF E ** 1 CALL COPY(X,AA,K) CALL COPY(X,A,K) CALL SET(0NE,1,£5,K) HALF(l) = 1 HALF(2) = -1 HALF(3) = 5000 DO 4 1=4, K 4 HALF( I ) = SW3 = 5 CALL SHORTN(X,XS,UO,K) 10 IF (XS.NE.O.O)GO TO 11 CALL COPY(ONE,W,K) RETURN 11 IF (XS.GE.l.O.OR.XS.LT.O.O)GO TO 14 SW2 = 2 12 CALL SBTRCT( A, HALF, A,K) IF ( A( 1 ) .GT.O)GO TO 13 CALL COPY(AA,A,K) GO TO 20 13 CALL COPY( AA, A,K) CALL MLTPLY( A, HALF, A, K) SW3 = 1 GO TO 20 14 IF (XS.NE.l)GO TO 110 CALL C0PY(E1,W,K) RETURN C C *** U( 1 ) = V( 1 ) = W( 1 ) = 1 20 CALL C0PY(0NE,U,K) CALL C0PY(0NE,V,K) CALL C0PY(0NE,W,K) SW = -1 CALL C0PY(0NE,B1 ,K) DO 70 1=1, 1000 C C *** U( 1+1 )=l/( 1+(A( I + 1)/(B( I )*B( 1 + 1 ) ) )*U( I ) ) IF (SW.GT.O)GO TO 40 io6 CALL SET(B2,It£30tK) 30 GO TO 50 40 CALL SET(B2,2,£50,K) 50 CALL MLTINT(A t -l,A,£60,K) 60 CALL MLTPLY(A,U»U,K) CALL MLTPLY ) J GO TO 120 CALL C0PY(B2,B1,K) sw=-sw 70 CONTINUE WRITE (6,1000) 1000 FORMAT (•-*** NO CONVERGENCE FOR EXP **»•) RETURN 110 SW2 = 1 INT = AINT(XS) CALL SET( A,INT,£120,K) CALL SBTRCT(AA,A,A,K) IF (A( 1).NE.0)GD TO 150 CALL COPY(ONE,W,K) 120 IF (SW3.NE.DG0 TO 130 CALL MLTPLY(W,W,W,K) 130 IF (SW2.EQ.2)RETURN CALL EXP0(E1, INT,A,K) CALL MLTPLY( A,W,W,K) RETURN 140 CALL C0PY(E1,W,K) RETURN 150 CALL C0PY(A,AA,K) GO TO 12 160 WRITE (6,1001) 1001 FORMAT (•-*** X IS OUT OF RANGE WHEN CALCULATING E**X IN LEXP ***• 1 ) RETURN END 107 ENTRY NAME LCOS IDENTIFICATION *** Multiple Precision Cosine PURPOSE *** Given a multiple precision number X of arbitrary array- length K(see description of Multiple Precision Basic Arithmetic Package for storage details), LCOS uses an infinite series to evaluate the cosine of X. OTHER SUBROUTINES USED *** IBCOM, COPY, SET, SBTRCT , MLTPLY , MLTINT, DIVIDE and ADD. USAGE *** INTEGER X(KX) ,Y(KX) ,K CALL LNGCNS(K,U,8,T5) CALL LC0S(X,Y,K) where the parameters in the calling sequence are: X - The multiple precision argument. input Y - On output Y=cos(X) where Y is a multiple precision number. output K - The array length of the multiple precision numbers X and Y. K.LE.KX.LE.50 input ERROR MESSAGES *** If the series for cosine does not converge in 1000 steps, control is returned to the calling program after the following message has been printed: 108 *** COSINE DOES NOT CONVERGE *** ACCURACY *** The maximum relative error observed in tests run was •20*10**(-((K-3)*^-l)). EXECUTION TIME *** Using the FORTRAN G compiler the execution time is given approximately by . 00068*K**3 seconds. SUBROUTINE LCOS(XA,Y,K) *** Y = COS I XA) , K.LE.50 109 1 INTEGER Y( 1 ) ,K ,NUM ( 50 ) t DENOM ( 50 ) , TERM ( 50 ) ,X2(50) ,FACT,FACT2, SW,0NF(50) , XA( 1) , X ( 50 ) , TEMP { 50 ) , S I GN INTEGER PI2(50),TERM2(50tl0) ,PIC(50) DATA PIC/ It 0, 6,2831,8530,7179,5864,7692,5286,7665,5900, 1 5768,3943,3879,8750,2116,4194,9889,1846,1563,2812,5724,1799,7256, 2 696,5068,4234,1359,6429,6173, 265,6461,3294,1876,8921,9101,1644, 3 6 345, 718,8162,5696,2234,9005,6820,5403,8770,4221,1119,2892,4588/ CALL COPY(XA,X,K) CALL SET(0NE,1,£10,K) 10 IF - The number of significant decimal digits desired. For maximum efficiency set D = (K-3)*U+1. input 113 K - The array length of the multiple precision numbers X and Y. K.LE.KX.LE.50 input ERROR MESSAGE *** If X is less than or equal to control is returned immediately to the calling program after the following message is printed: ***ERROR IN LGAM*** Z IS LESS THAN OR EQUAL TO WHEN CALCULATING GAMMA (X) ACCURACY *** The maximum relative error observed in tests run with D set equal to (K-3)* 1 4+l was .80*10**(-(K-U)*U) . If D is set so as to be less than (K-M*U there should be D significant decimal digits in the answer. EXECUTION TIME *** Using the FORTRAN G compiler the execution time in seconds is given approximately by .0059*K**3. Ill* ,TEMP3(50) ,C(50,71 ) , B(2I) ARE THE SUBROUTINE LGAM( I IN ,Y ,D,K ) *** Y = GAMMA(ZIN) FOR ZIN.GT.O, K.LE.50, D ( NO .S IGN.OIG I TS ) .LE . ( K-3 ) *4 COMMON /LONG/LN,LW,LT,LS,LR,LTR,LEX,LRM INTEGER LN,LW,LT,LS,LEX REAL LRfLTRtLRM INTEGER ZIN( 1)»Y( 1) , Z ( 50) , TEMP ( 50 ) , T6MP2 ( 50 ) 1 LN2PI (50) ,0NE(50) ,D,K , SW ,T ( 50 ) *** SET CONSTANTS Cdl = B ( 2 I ) / ( 2 I * ( 2 I- 1 ) ) WHERE BERNOULLI NUMBERS INTEGER C1(250),C2(250),C3(25 3),C4(2 50),C5(2 1 C8(250),C9(250) ,C10<250> ,C1 1(250) ,C12(250) , EQUIVALENCE ( C ( 1 , 6 ) , CI ( 1 ) ) , (C ( 1 1 1 1 ) t C2 ( 1 ) ) , ( 1(C( 1,21 ) ,C4( 1) ), (C( 1,26) ,C5(1) ), (C( 1,31) ,C6( 28009, 9502, ,9950, ,2487, ,5621, ,8905, ,4726 ,3682, ,4848i 4539, 3751, .6520, .1458, 1055, , 734, ,7819, .9014, . 226, .5149, .9687, , 985, 8536, 3741, .5180, .5872, 5652, ,9516, •1659, .4604, ,9595, .7740, .2933/ ,3668, 538, 7161, 7263, 2093, 5756, ,4499, ,7330, .1502, . 848, .7282, .1186, f 1362, 7669, 6163, 7763, 9570, 8725, , 983, .1819, .8179, , 794, 7, .5166, f 3357, 1418, 375, 4804, 9276, 4181, ,2511, 4155, .2511, .4155, 2511, ,4155, ,4155, 2511, 4155, 2511, 4155, 2511, ,2511, 4155, 2511, 4155, 2511, ,4155, ,6220, 6261, 8053, 7780, 1511, 8128, ,1111, 1111, 1111, 1111, 1111, 1111, ,1111, 1111, 1111, 1111, 1111, 1111, ,1111, 1111, 1111, 1111, 1111, 1111, ♦2089, 5654, 6242, 4808, 2412, 6356, ,1134, 5696, 6983, 1875, 9801, 5807, 7807, 9439, h291, 6235, 4714, OAT A 6567, 5215, 6156, R780, *398, 151, 761 , ,M56, «627, HI 71 , 3605, ^795, 15823, 6197, 8105, DATA 5167, 8184, 5747, •3019, 394, 2329, 1770, 8347, 3894, '8046, 2096, <-4fe7, "♦629, 3 44, ■ 647, DATA 1013, 2782, ■5320, '1208 , '9266, '8974, I '8896, '1440, 2640,6404 9171 ,2526 -1 8069,8304 6067, 724 8032,6391 C8/ 1 6296,2804 4296, 1971 9735,7330 -1 7044,7003 5768,4843 4785,6013 1 193, 840 8627,4509 4509,8039 -1 2347,5688 7439,4450 5423, 375 1 3629, 276 6751,8605 3257,9136 C9/ -1 999, 311 854, 1812 7983, 98 1 3480, 427 2652,3297 7491, 394 -1 3911,5939 4051,8631 5110,4397 1 2750,2271 2199,3127 3539,5189 -1 1381,9119 8829,7314 9132,7617 CIO/ 1 925,7734 5611,9648 7026,6666 -1 6697,9299 438,4197 7817, 605 1 1540,5955 4391, 399 ,5257 ,9747 t 12 ,2075 ,8281 ,9730 , 13 ,4555 ,3421 ,5567 » 13 ,8282 ,3095 ,4377 » 14 ,2505 ,8039 ,2156 , 15 ,8605 ,7147 ,4489 , 15 ,9152 ,7704 ,7644 , 16 ,9478 ,7971 ,9790 , 16 ,5712 ,4910 ,2652 , 17 ,6665 ,8722 ,9621 , 18 ,7730 ,1477 , 34 , 18 ,2049 ,5799 ,6102 t 19 ,3048 , 258 , 563 , 19 ,2240 ,3499 ,1169 , 20 ,7511 , 505 ,5353 ,3312 ,3167 ,3030 ,4451 ,5125 , 50 ,6213 ,5875 ,9508 ,8529 ,4140 ,5086 ,1244 , 150 ,4563 ,2156 ,8627 , 2 ,4688 ,1754 ,7963 , 540 ,4237 ,5140 ,8233 t 10 ,2603 ,3194 ,3514 ,2327 ,3111 ,3942 ,3297 * 51 ,6336 ,4331 ,4268 , 1 ,9854 ,6632 ,3642 , 286 ,3144 ,4284 ,4587 , 7 ,7969 ,2893 ,4211 ,1876 ,1628 ,4199 ,3228 , 50 ,6088 ,8596 ,1178 ,3015 , 226 ,3915 ,1053 ,4238 ,7000 t 731 ,8490 ,1560 ,9728 , 895 ,9312 ,1310 ,6417 ,9808 ,8627 ,4509 ,7893 ,1031 ,7593 , 35 ,9350 ,8260 ,4834 ,2098 ,9753 ,8689 ,6068 ,9247 ,4876 ,3439 ,6523 ,4910 ,5392 ,7888 ,1393 ,7813 , 1906 ,7424 ,3024 ,6116 ,6893 ,2655 ,2769 ,3072 ,1893 ,9797 ,5182 ,1326 , 693 ,4037 ,1817 ,3079 ,9049 ,2424 ,6473 ,4859 ,2206 ,6348 ,4461 ,1793 ,2582 ,6461 ,9401 ,7531 ,8111 ,2030 ,4544 ,1966 ,7502 ,2809 ,4127 ,4509 ,8039 ,4947 ,9303 ,2445 ,1696 ,4352 ,3718 ,5718 ,3638 ,3782 ,6852 ,1162 ,9152 ,2026 ,9270 ,2974 ,3942 ,9162 ,1319 , 353 ,7484 ,2102 ,5633 , 549 ,8384 ,8960 , 825 ,1253 ,1557 , 780 ,2830 , 64 ,2505 ,4305 ,8099 ,7932 ,9291 ,1469 ,7587 ,8169 »7052, ,6440, ,8666, ,6387, t 590, ,1552, »2lil, » 610, ,3600, ,4581, t 551, ,6374, ,699 7, ,2978, ,3405, ,4509, ,8039, ,2156, , 383, ,8537, »7621, ,8107, ,8604, ,8555, ,8594, ,4136, ,1508, ,4499, ,8088, , 536, ,1847, ,3512, t9103, ,6523, , 653, »6975, , 864, ,2464, ,3089, ,6568, ,8281, ,8797, ,2966, ,3986, ,9738, , 41, ,2337, ,5521, ,1145, ,1686, , 467, ,2271, ,4662, ,6895, , 751, ,1088, ,4110, 9370,1677 6921,1881 1827,4134 695,7705 7998,4964 624, 345 3734,3179 3426, 498 793,2966 7839,2582 8816,2084 8834,9557 2371,2423 687,7334 9857,6695 8039,2156 2156,8627 8627,4509 1636,8712 2100,9345 3765,7762 7064,4640 1500,5763 3658, 622 4427,4831 7599,4035 5198,5501 9327,5855 2238,3290 2013,8726 9173,4786 7249,1039 9426,5232 2974,9103 2139, 194 52,9993 6045,6158 4980,7715 226,4576 1311,7102 7369,4158 2508,5910 7369,6226 3500,1121 8223,6708 8526,7011 2986,4392 9845,9798 9453,8510 8607,8039 1847,8954 6173,4594 6366,2453 6463,3327 704,8132 5216,6093 6660,7328 651 1467 9556 1960 8804 813 2648 7852 4336 6763 52 3840 3826 7493 1173 8627 4509 8039 8838 5145 8788 2267 5618 9009 2380 104 6788 3503 4683 5217 4103 4265 9749 9426 6552 7743 1996 83 8481 5008 756 6529 4205 9220 5193 5496 8443 4974 7165 1492 6272 6573 4408 2272 9151 4162 2895 ,4903 ,8168 ,7742 ,8668 ,4823 , 906 ,1531 ,6451 ,3577 ,8455 ,2162 ,8707 ,8191 ,9119 ,6037 ,4509 ,8039 ,2156 , 1686 , 866 ,6019 , 312 ,7188 ,9069 ,7116 ,5972 ,7261 ,2343 ,9685 , 524 ,2052 ,2329 ,1039 ,5232 ,1217 ,9293 , 717 ,1505 ,6176 ,5910 , 137 ,2096 ,4190 ,3375 ,3678 ,3981 ,7139 ,8451 ,3628 ,6860 ,6582 ,8315 , 112 , 889 ,8692 ,4240 ,7128 117 6215 9277 5613. 6379 5877 8015 7487 4849 7673. 548 2788 1170 7106 6051 9879. 8039 2156. 8627 3127 3271 5370 7837 4152 2524 1002. 5746 7079 5429 1375 5110 1849 7491 4265 9749 1717 9556 9850. 2999 8713 6529 4570 2199 772 1860 2163. 2466, 3137 7706 6083. 2311 6766 225 3526 517 8962 6469 225 648 58,3296,2468,4764,8 716,807 5,95 31,2268,2094,6447 C -1, 21, 1,4351,4288,2843,2171, 936 D 752,9850,2618,6805,5957,7080,3259,5029, 670,4580 66 306,1110,6478,9732,2205,4882,3059,988 2,2 317,4848 F 1378, 942 7, 4487,5797,5 194, 1032, 4036, 138 7, 65 88, 895 G 1, 21, 420, 905,7506,6658,6940,7247 H2330, 6286, 4274, 3146, 729, 680,2763,3238,3912,4076 16486,2654,5397,4051,718 5,4894,4498,1050,3043,2694 J 2 64 5, 11 15,938 7,3243, 7 36,6371,4052,443 7,9099,4032 DATA CI 1/ -1, 22, 12,7585,8127,2247,5759, 843 12 045,3193,1983, 98,9091,7026,45 30,6582,1855,8047 28557,3617,3869,7649,5244,4591,2821,2202,2217,8 774 33312, 963,1430,4216, 795,1271, 218,7002,7416,7074 4 1, 22,4017,7568,4036,2181,8804,5903 55528,7105,4463, 534,9215,232 2,3640,7637, 969,5373 65939,4848,6444,7253,8 277,8 480,1236,1286,6975,698 8 73 14 5,1046,495 2,942 8,2904,902 3,7 392,89 22,6014,8897 8 -1, 23, 131,1013,6312, 68,3458,7473 978 79,6212,605 8,3407,38 28,5 32 7,8464,8866,2 82 7,8032 A9 5 68, 5028, 4498, 5 384, 4264 ,308 5, 2 370, 4249, 5 70 1,30 18 B4618, 299,6806,6814, 506, 181,7735,2001,9651,1913 C 1, 24, 4,4299,9565,3697,8259,6377 D5691, 7838, 4230, 3837, 6938, 826, 411,4234,9237,6322 61628,3203,5794,7727,9202,2792, 227,9202,2 792, 227 F9202,2792, 227,9202,2792, 227,9202,2792, 227,9202 G -1, 24,1549,2140,6951,7355,3881,1542 H2 547, 1418, 1057,3373,9 601 ,2 370,8171 ,7806,5382,5930 12430,9461,8694,4563,3172,7198,5339, 496,4913,4777 J9200, 4079, 920 1,4 147, 582,2863,9475,2408,6179,3857 DATA C12/ 1, 25, 56, 376,4855,6273,5197,4950 18800,4366,9592, 900,6681,2345,8103,9077,2556, 238 29690,62 5 6,8348,87 83,8643,2506,88 70,52 34,1597,7961 33 415,9779,6143,2 506,88 70,52 34,1597,7961,43 25, 6 88 4 -1, 26, 2, 953,9241,4858, 979,2706 52855, 5944, 3198, 7352, 4201, 5953, 21 59, 9949, 5931, 8851 65634,1485,8564,4000,5889, 834, 108,4010,8401, 840 78401, 840,1084, 108,4010,8401, 840,1084, 108,4010 1, 26, 809,5290,3005,6367,8084, 63 9 5997,2820, 1881,8398,3044,5464, 3 64,1512,5125,8306 A 60 75,63 70, 3 024, 4291, 7977, 2 38 6, 6286, 19 7,7379,7685 B2 79 0,4718,4731, 978,6496,7147,2187,3847,4849,8207 -1, 27, 32,2963,3556, 269,5213,4559 D 9 58,98 12,3402,9130,9 321,2400,352 7,6342,43 9 8,2642 61956,9890,5 726,7601,3392,4421,7179,5596,4373,9385 F9 542, 998 3, 168,2 877,8 755,5967,2687,9 728,2692,6046 G 1, 28, 1,3298,6391,7366,8232, 869 H 658,5413,4002, 922,8311,5936,1392,4138,9612,6198 1 583, 284, e428, 4891, 9580, 9530, 6273,6393, 42,4805 J 341 3,6859,45 19,66 8 5,4771,7544,6996,1464,45 81,8024 DATA C13/ -1, 28, 564,9116,4660,5107,5898,6024 12188,1037,6286,7139, 178,7531,9053, 990,9775,4518 2 6028, 701 1 ,302 7,87 5 2,3561,4806,358 3,2 703,1641,8 990 39 3 59,89 17,5 354,9 089, 1743,1518,4985,3654,4714,5494 4 1, 29, 24,7437,7846,3950,7093,8854 53718,9292,7520,2293,8585,6543,58 42,85 72, 5 30,4828 6 559,2257,2516,5851,8167,5704,1250, 319,7374,6867 7 7 468,6716,7919,7994,9874,6867,1679,1979,9498,7468 118 1146,7141,3737, 5153,4591,8310, 9720,1892,8967, 7458,9287,8560, 3621,9098,7433, 2242,8288,5867, 2551,2366,9039, 7906,7426,6324, 1370,1360,7079/ 2594,5728,4945, 8594,3538,3133, 6568,3026,7234, 8635,5119,5921, 1259,3158,7661, 1292,9279,3966, 3410,5913,7519, 3170,3890,9959, 7825,1695,5141, 4125,9381,8085, 9142,7167,7720, 5347,5804,4706, 1152, 544,8820, 4420,1146,7108, 9202,2792, 227, 2792, 227,9202, 8813,4953, 339, 5723,1506,3202, 4580,7540,9157, 8358,8558,9707/ 6206,7876,8770, 1118,5400,4688, 4325, 688,7052, 7052,3415,9780, 4754,4170,7191, 1765,8764,7819, 1084, 108,4010, 8401, 840,1084, 6463,6778,8443, 8966,6171, 379, 8432,7087,3416, 8470, 99,8163, 7210, 298,7523, 2899,6957,3269, 5179,8672,2247, 89,5476,3008, 8228,8729,1823, 2809,3169,9242, 4842,7645,8890, 179,1750,2739/ 4699,2770, 809, 4469,3938,9072, 6961,2807,5939, 7996,8980,5411, 9333,8990, 862, 4601,9346,2012, 1679, 1979,9498, 6716,7919,7995, 119 1000 10 20 40 so 60 8 -1, 26, 4,1736, 993, 9 959, 679,6456,7119,5997,1118,9800, A8 446, 94 78, 2 023, 65 5 0,649 3, 9477, 4 104, B6 19 5, 286 1,952 8, 6 195, 286 1,9 528, 61 95, C 1, 30, 519,4881,6925, 02830,2753,5169,1125, 769,9054,2487, El 990, 4337, 5511, 1986,2829,4235,7 561, F490 1,1 5 14, 7249,4491 ,2616,4868, 301 , G -1, 31, 24,8780* 657, H66 50, 11 84, 7376, 2 7 56, 7422, 3284, 4736, 1 645,1038,8936,8521,6462,1414,5918, J1297,7222, 900,4839, 445,6636,9069, DATA C14/ 1, 32, 1,2263, 550, 17272,9834, 300,8689,6615,9526,1452, 2 6414,7680,9821,55 21,4343,2515,4458, 3 5011,82 03,3096,9 26 7,1394,7990,5437, DATA LN2PI/ 1, -1,9189,3853,3204 1 8613, 9747, 3 637, 7834, 1281, 7151 ,5 404 2 9863, 5954, 19 76,2200,5646,6246,3433 3 58 75,9152, 2268, 139 3,6035,6074,2547 IF(ZIN( 1 ).GT.O) GO TO 4 WRITE (6,1000) FORMAT ( •- ***ERROR I 1LCULATING GAMMA( X) • ) RETURN CALL COPY( 2 IN,Z,K) SW = CALL SET(0NE,1 ,f.5,K) I T = D - 10 IEND=20 IF (D.GT. 50) IEND = 3 I F ( D.GT.80) IENO = 4 IF (D.GT. 110) IEND = IF(D.GT.137) IEND=61 IF (D.GT. 165) IEND = IF (D.GT.16)Gn TO 10 IT = 7 CALL SFT( T, IT,C20,K ) CALL SBTRCT ( Z,T,TEMP2,K) IF (TEMP? ( I ) .GF.O)Gn TO 60 SCALE Z SO CALCULATION IS LN GAMMA(Z) AS 1=0, K-l OF (Z+I ) ) CALL COPY(Z,TEMP,K> CALL ADD(Z ,ONF, Z,K) CALL ADD( TFMP2,0NE ,TEMP2,K ) IF , SUM ( 50 ) ,R ( 50 ) , S ( 50 ) , 1 RR< KX,1), TEMP(50) , TEMP2 ( 50 ) , Z ERO ( 50 ) , ONE < 50 ) , TWO ( 50 ) CALL LNGCNS(K,4,8,75) INITIALIZE CONSTANTS FOR LONG ARITHMETIC CALL SET! ZERO, 0, CIO, K) CALL SET(ONE, 1,620,K) CALL SET(TW0,2,&30,K) CHFCK IF INPUT PARAMETERS ARE IN CORRECT RANGE CALL SBTRCT( X, ONE, TEMP, K) IF(TEMP( 1) .LT.O.OR.A.LT.O.OR.NX.LE.O) RETURN 1 40 45 * ** 50 60 65 90 IF X = l OR A=0 THEN P(0)=1 AND P(N)=0, IF (TEMPI 1 ) .NF.O.AND.A.NE.O) GO TO 50 CALL COPY(DNE,P ( 1, 1 ) ,K) IF (NX.LT.2)G0 TO 45 00 40 N=2,NX CALL COPY! ZFRO,P ( 1 ,N) ,K) CONTINUE RETURN P(N)=0 FOR N=A+1,NMAX IF (A+2.GT.NXJG0 TO 65 AP2 = A + 2 00 60 N=AP2,NX CALL COPY(ZERO,P ( 1 ,N) ,K) CONTINUE XI = S0RT(X**2-1) CALL MLTPLYI X,X,TEMP ,K) CALL SBTRCTI TEMP, ONF, TEMP, K ) CALL LSQRT(TEMP,X1 ,K) C = A FACTORIAL CALL COPY(ONE,C,K) I F ( A.LT.2 )G0 TO 95 00 90 N=2,A CALL MLT INT (C,N,C, £300, K ) CONT INUE N=1,NMAX *** SUM = (X + XI ) **A / C 95 CALL A00( X, XI ,TEMP ,K ) CALL FXPOf TEMP, A,TEMP,K ) CALL DIVIDE(TEMP,C,SUM,K) l r)o X1=2*X/X1 CALL NLTINT( X,2,TFMP,r,100,K ) CALL 01 VlDE< TEMP ,X1 ,X1 ,K ) 12U c c *** c c *** 110 120 140 c c *** c c *** 150 C C *** 160 170 C C * * * 200 C c #** son l ooo R = S = CALL COPY(ZERO,R,K) CALL C0PY(ZERO,S,K> 00 150 NN=1,A N = A + 1 - NN R = (A+l-N) / (N*X1 + (N+A+l) * R) CALL MLTINT(Xl,N,TEMP,£300 t K) NPAP1 = N + A + 1 CALL MLTINT(R,NPAP1,TEMP2,£300,K) CALL ADD(TEMP,TEMP2,TEMP,K) AP1MN = A + 1 - N CALL SET(TEMP2,AP1MN,£300,K) CALL DIVI0E(TEMP2,TEMP,R,K) S = R * (2 + S) CALL ADD( S, TWO, TEMP, K> CALL MLTPLY(RtTEMPtStK) IF (N.GE.NX)GO TO 150 RR(N-1 ) = R CALL COPY(R,RR< 1,N) ,K) CONTINUE P(0) = C * SUM / ( 1 + S) CALL MLTPLY(Ct SUM, TEMP, K) CALL A0D(ONE,S,TEMP2,K) CALL DIVIDE(TEMP,TEMP2,P( 1,1 ) ,K) IP (NX.GT.A+1 )G0 TO 160 END = NX-1 GO TO 170 END = A 00 200 N=1,EMD MP1 = N+l P(N+1) = (N+A+l) * RR(N) * P(N) CALL MLTPLY(P( 1,N) ,RR( 1,N) ,TEMP,K) NPAP1 = MPl + A - 1 CALL MLTI MT( TEMP, NPAPl,P(l,NPl),f» 300, K) CONTINUE RETURN FRROR RETURN WRITE ( 6, 1000) FORMAT ( ' - *** RETURN1 END ERROR IN ILEG1 *** NUMBERS HAVE BECOME TOO LARGE* > 125 ENTRY NAME ILEG2 IDENTIFICATION *** Multiple Precision Associated Legendre Function of the Second Kind for Integer Order and Degree. PURPOSE *** This subroutine generates, in multiple precision, to D significant digits the associated Legendre function of the second kind Q(X), for integer order M.GE.O and degree N=0,1, . . . ,NX-1 and X.GT.l. OTHER SUBROUTINES USED *** FIXPI , IBCOM, LNGCNS , SET, SBTRCT , MLTPLY, LSQRT, ADD, DIVIDE, LNGLOG, COPY, RECIP and MLTINT. USAGE *** INTEGER X(KX) ,M,NX,D,Q(KX,NXX) ,K,KX,QQ(KX,NXX) ,RR(KX,NXX) "CALL ILEG2(X,M,NX,D,Q,K,&L,KX,QQ,RR) vhere the parameters in the calling sequence are: X - The multiple precision argument stored in the form X.EQ.A*(l0**U)**B where ■ 1 . LE . A . LT . 10**U as ' follows : X(l) = sign of the number: +1 if positive -1 if negative if zero (in this case, X(2) to X(K) are disregarded) X(2) = exponent B of 10**U X(3) = integral part of the mantissa, l.LE.X(3).LT.10**U X(h) to X(K) = fractional part of the mantissa, stored h decimal digits per location. The value of X must be greater than one. input 126 M - The integer order of the Legendre functions which will be generated. M.GE.O input NX - The Legendre functions will he generated for degree 0, 1, . .., NX-1. NX.GT.O, NX.LE.NXX input D - The number of significant digits desired in the function values. A suggested value to use is (K-U)*U+3. input Q - On output, the Jth column, Q(l,J), 1=1, ...» K will contain the multiple precision Legendre function of order M and degree J-l stored in the same manner as X. output K - The number of locations used in the arrays to store the multiple precision numbers. K.LE.KX.LE. 50 input L - Control returns immediately to the statement labelled L (an actual integer) if X.LE.l or M.LT.O or NX.LE.O. Also see ERROR MESSAGES. input KX - The maximum value of K.KX.LE-50 input QQ,RR - Arrays used for temporary storage in ILEG2. ERROR MESSAGES *** In some cases the program will be unable to continue execution and control will return to the label parameter L 127 specified in the calling sequence after the folio-wing message has been printed: *** ERROR IN ILEG2 *** NUMBERS HAVE BECOME TOO LARGE This can occur if M or NX are too large or if the number of significant digits D asked for is too great and the algorithm is failing to converge. ACCURACY *** There should be D significant digits in the generated function values. However if D is set equal to (K-3)*^ there may only be D-l significant digits. Convergence is slow for X near 1, 128 SUBROUTINE I L EG2 ( X ,M ,NX ,D ,0 ,K , *, KX , QAPPX,RR) C *** THIS SUBROUTINE GENERATES* IN MULTIPLE PRECISION, TO SIGNIFICANT C *** DIGITS THE ASSOCIATED LEGENDRE FUNCTIONS OF THE SECOND KIND QUI, C *** FOR INTEGER ORDER M.GE.O AND DEGREE N«0, 1 , . . . ,NX-1 AND X.GT.l INTEGER X(l) ,M,NX,D,K,Q(KX,1),N»NU,NN,QAPRX(KX,1) ,RR(KX,1 ), 1 X1(50),Q0(50),QH50),Q2(50),R(50),EPS(50),TEMP(50), 2 TEMP2(50) ,ZERO(50>,ONE(50),HALF(50) CALL LNGCNS(K,4,8,75) C C *** INITIALIZE CONSTANTS FOR LONG ARITHMETIC CALL SET(ZER0,0,£20,K) 20 CALL SET(ONE,1,G30,K) 30 HALF(l) = 1 HALF(2) = -1 HALFO) = 5000 DO 35 1=4, K 35 HALF( I ) = C C *** CHECK IF INPUT PARAMETERS IN CORRECT RANGE CALL SBTRCTJX, ONE, TEMP, K) IF (TEMP( 1) ,LE.0.OR.NX.LE.O.OR.M.LT.0)RETURN 1 C C *** XI = SORT( X**2-l ) 40 CALL MLTPLY(X,X,TEMP,K> CALL SBTRCT(TEMP, ONE, TEMP, K) CALL ISQRT(TEMP,X1,K) C C *** 01 = .5 * LN< (X+1)/(X-1 ) ) CALL AOD( X, ONE, TEMP, K) CALL SBTRCT(X,0NE,TEMP2,K) CALL 0IVI0E(TEMP,TEMP2,TEMP,K) CALL LMGL0G( TEMP, TFMP , £50, K ) 50 CALL MLTPLY(TEMP,HALF ,01 ,K ) IF (M.NE.O)GO TO 60 C C *** 0(0 ) = 01 CALL C0PY(01,0( 1 ,1) ,K) GO TO 160 C C *** 02 = -1/X1 60 CALL RECIP(X1,02,K) 02(1 ) = -1 * 02( 1) C C *** XI = 2*X/X1 CALL MLTINT( X ,2 ♦ TEMP , £80 ,K ) 80 CALL MLTPLY(TEMP,02,X1,K) XI ( 1 ) = -1 * Xl( 1 ) MMl = M-l I F (MMl .EO.O)GO TO 155 DO 150 N=1,MM1 C C *#* 00 = 01 CALL C0PY(01 ,00,K) c c *** 01 = 02 CALL C0PY(02,01,K) 129 ««* 02 = -N*X1*Q1-N*(N-1)*Q0 NNM1 = N * (N-l) CALL MLTINT(00,NNMl f TEMP2»£500 t K) CALL MLTINT(X1,-N,TEMP,C500,K) CALL MLTPLY(TEMP,Q1,TEMP,K) CALL SBTRCT(TEMP,TEMP2tQ2tK) 150 CONTINUE *** 0(0) = 02 155 CALL C0PY(02»0(1,1) t K) 160 DO 170 N=ltNX *** OAPRX(N) = CALL COPY(ZERO,OAPRX( 1 t N) ,K) 170 CONTINUE *«* EPS = .5 * lO**(-0> FPS( 1 ) = +1 EPS(2) = -D/4 - 1 FPS(3) = 5 * 10**(-4*EPS(2) - D 00 180 I=4 t K IPO FPS( I ) = NU = 20 + I F I X ( 1.25*(NX-1 ) ) NUINC = 10 - 1) *** R = 210 CALL COPY(ZFRO,R»K) DO 250 NN=l f NU N = NU-NN+1 *** R = (N+M)/ ( ( 2*N+1 )*X-(N-M+1 ) *R ) N2P1 = 2*N+1 CALL MLTIMT( X, N2P 1 , TEMP , £500 ,K ) NM1 = N-M + l CALL MLTINT(R,NMl,TEMP2t&500,K) CALL SBTRCT (TEMP ♦ TEMP 2 tT EMP2 ,K ) NM = N+M CALL RECIP(TEMP2tTEMP2,K) CALL MLTINT(TEMP2tNM,R ? &500tK) IF (N.GE.NX)GO TO 250 *** RR(N-1 ) = R CALL COPY(R,RR( 1,N) »K) 250 CONTINUE NXM1 = NX-1 *** Q(N+1 ) = RR(N) * 0(N) DO 260 N=1»NXM1 CALL MLTPLY(RR( l t N) ,0( 1»N) t 0( 1»N+1 ) t K) 260 CONTINUE DO 280 N=l t NX *** IF (ABS(O(N)-OAPRX(N) ) .GT . EPS*ABS ( ( N ) ) GO TO 290 CALL SBTRCT ( 0( 1,N) ,OAPRX( 1,N) ,TEMP,K) TEMPI! ) = +i CALL MLTPLY(EPS,0( ltN) »TEMP2,K) TEMP?( 1 ) = +1 130 CALL SBTRCT(TEMP,TEMP2,TEMP,K) IF(TEMP(1 ).GT.O) GO TO 290 280 CONTINUE RETURN 290 DO 300 N*1,NX C C *** OAPRX(N) = 0(N) CALL COPY(Q(l t N) t OAPRX(l f N),K) 300 CONTINUE NU = NU ♦ NUINC NUINC = 2 * NUINC GO TO 210 C C *** ERROR RETURN 500 WRITE (6,100> 1 00 FORMAT (•- *** ERROR IN ILEG2 *** NUMBERS HAVE BECOME TOO LARGE 1 ) RETURN1 END 131 ENTRY NAME ILEG3 IDENTIFICATION *** Multiple Precision Associated Legendre Function of the Second Kind for Integer Order and Degree. PURPOSE *** This subroutine generates, in multiple precision, to D significant digits the associated Legendre functions of the second kind Q(X), for order M=0,1, . . . ,MX-1 and integer degree N.GE.O and X.GT.l. OTHER SUBROUTINES USED *** IBCOM, ILEG2, COPY, MLTPLY, SET, SBTRCT , LSQRT, MLTINT and DIVIDE. USAGE *** INTEGER X(KX) ,N,MX,D,Q(KX,MXX) ,KX,QQ(KX,MXX) ,RR(KX,MXX) ,K CALL ILEG3(X,N,MX,D,Q,K,&L,KX,QQ,RR) where the parameters in the calling sequence are : X - The multiple precision argument stored in the form X.EQ.A*(10**>0**B where l.LE.A.LT.10**U as follows: X(l) = sign of the number: +1 if positive -1 if negative if zero (in this case, X(2) to X(K) are disregarded) X(2) = exponent B of 10**U X(3) = integral part of the mantissa, l.LE.X(3) .LT.10**U X(U) to X(K) = fractional part of the mantissa, stored k decimal digits per location. The value of X must be greater than one. input 132 N - The integer degree of the Legendre functions which will he generated. N.GE.O input MX - The Legendre functions will he generated for order 0,1, . . . ,MX-1. MX.GT.O, MX.LE.MXX input D - The number of significant digits desired in the function values. A suggested value to use is (K-U)*H+3. input Q - On output, the Jth column, Q(l,J), 1=1, ...K will contain the multiple precision Legendre function of degree N and order J-l stored in the same manner as X. output K - The number of locations used in the arrays to store the multiple precision numbers. K.LE.KX.LE.50 input L - Control returns immediately to the statement labelled L (an actual integer) if X.LE.l or N.LT.O or MX.LE.O. Also see ERROR MESSAGES. input KX - The maximum value of K. KX.LE.50 input QQ, RR - Arrays used for temporary storage in ILEG3. ERROR MESSAGES *** If the parameters N and MX are too large the program will be unable to continue execution and control will return to the label parameter L specified in the calling sequence after the follow- ':ssage has been printed: 133 *** ERROR IN ILEG3 *** NUMBERS HAVE BECOME TOO LARGE Also see ERROR MESSAGES in the description of ILEG2. ACCURACY *** There should be D significant digits in the generated function values. However if D is set equal to (K-3)*U there may- only be D-l significant digits. Convergence is slow for X near 1. 13fc SUBROUTINE I LBG3( X,N,MX t D,Q f K ,* t KX t QAPRX,RR ) C *** THIS SUBROUTINE GENERATESt IN MULTIPLE PRECISION, TO D SIGNIFICANT C *** DIGITS THE ASSOCIATED LEGENDRE FUNCTIONS OF THE SECOND KIND O(X), C *** FOR ORDER M=0 , 1 ♦ . . . ,MX-1 AND INTEGER DEGREE N.GE.O AND X.GT.l COMMON /LONG/LN,LW,LT,LS,LRtLTR,LEX,LRM INTEGER LN,LW,LT,LS,LEX REAL LR.LTRtLRM INTEGER X(1),N,MX,D,K,KX,0(KX,1),QAPRX IF (TEMP(1 ).GT.O)GO TO 270 240 CONTINUE C C *** P1I0I = C*P1<0) CALL MLTPLY(CtPl(ltl)fPl(ltl)*K) IF ,N ,NU , M , 1 EPS(50).T(50)tXl(50)tX2(50)«SUM(50),LMBDAl(50)*LMIDA2(50)« 2 LMtDA(50) ,R(5O)tS(50)»ZERO(50)»ONE(50>tTWO(50)tTEMR(50)t 3 TEMP2(50),TEMP3<50),HALF<50) REAL XStTAUS CALL LNGCNS(K,4,8,75) C C *** INITIALIZE CONSTANTS FOR LONG ARITHMETIC CALL SET(ZER0,0,U0,K) 10 CALL SET(0NEtl,620,K) 20 CALL SET(TWO,2,S30,K) 30 HALF(l) = 1 HALF(2) = -1 HALF(3) = 5000 DO 35 1=4, K 35 HALF( I ) = C C *** CHECK IF INPUT PARAMETERS ARE IN CORRECT RANGE CALL SBTRCTU, ONE, TEMP, K) IF (TEMPd ).LT.O.OR.NX.LE.O)RETURN 1 C C *** IF X = l, P(0)=1, P(N)=0 N=1,NMAX IF (TEMPI 1 ) ,GT.O)GO TO 50 CALL COPY(ONEtP(lf 1 ItK) IF (NX.EQ.l)RETURN 00 40 N=2tNX CALL COPY(ZERO,P( 1 ,N) ,K) 40 CONTINUE RETURN C C *** T = TAU**2 50 CALL MLTPLY(TAU,TAU,T,K) C C *** PAPPROX(N)=0 FOR N=0,NMAX DO 60 N=1,NX CALL COPY(ZERO,PAPRX( 1,N) ,K) 60 CONTINUE C C *** FPSILON = ,5*10**<-D) EPSd ) = +1 EPSI2) = -D/4 - 1 EPS(3) = 5 * 10**(-4*EPS(2) - D - 1) on 70 1=4, K 70 EPS( I ) = C C *** XI = S0RT(X*«2-1) CALL MLTPLYi X,X,TEMP,K) CALL SBTRCT( TEMP, ONE, TEMP, K ) CALL LSORT( TFMP,X1 ,K) ll+9 C *** X2 = X + XI CALL ADD(X,X1,X2,K) C C *** SUM « C0S(TAU*LN(X2) )/SQRT(X2) CALL LNGL0G(X2,TEMP,S80,K) 80 CALL MLTPLY(TEMP,TAU,TEMP,K) CALL LCOS(TEMP,TEMP,K) CALL LSQRT(X2,TEMP2 f K) CALL 0IVlDE(TEMP,TEMP2,SUMtK) C C *** XI ■ 2*X/X1 * XI ■ 2*X/X1 CALL MLTlNT(X,2 f TEMP,C90,K) CALL DlVIDE(TEMP,XltXl,K) C C »** NU=30+ENTIER( ( l+( . 140+ .0246*TAU > *( X-i ) )*NMAX) CALL SHORTN(X,XS,£350,K) CALL SH0RTN(TAU,TAUS,£350,K) NU = 30+AINT( ( 1+(.140+.0246*TAUS)* CALL MLTPLYUMBDA,TEMP2tLMBDA,K) CALL MLTPLY(TEMP,HALF,TEMP2,K> CALL ADD(ONE,TEMP2,TEMP2,K) CALL MLTPLY(TEMP2,TEMP2,TEMP2,K> CALL MLTPLY(TEMP,TAU,TEMP,K) CALL MLTPLY(TEMP, TEMP, TEMP, K) CALL ADDITEMP2, TEMP, TEMP, K) CALL DIVIDEUMBDA, TEMP, LMIDA,K ) C C *** IF N.LT.NU THEN LAMBDA1«LAMBDA2 ; LAM10A2-LAMIDA ; N»N+1; GO TO LI IF (N.GE.NU)GO TO 200 CALL C0PY(LMB0A2,LMBDA1,K) CALL C0PY(LMBDA,LMBDA2,K) N = N+l GO TO 160 C c *** R = s = 200 CALL COPY(ZERO,R,K) CALL COPY(ZERO,S,K) C C *** L2:R=-( ( l-.5/N)**2+(TAU/N)**2 )/(Xl+(l+l/N)*R) 210 CALL SET(TEMP3,N,&500,K) CALL RECIP(TEMP3,TEMP3,K) CALL ADD(TEMP3,0NE,TEMP2,K) CALL MLTPLY(TEMP2,R,R,K) CALL ADD(R,X1,R,K) CALL MLTPLY(HALF,TEMP3,TEMP2,K ) CALL SBTRCT(0NE,TEMP2,TEMP2,K) CALL MLTPLY(TEMP2»TEMP2tTEMP2rK ) CALL MLTPLY(TEMP3,TAU,TEMP,K) CALL MLTPLY(TEMP, TEMP, TEMP, K) CALL A0D(TEMP2,TEMP,TFMP,K) CALL DIVI0E(TEMP,R,R,K) R(l ) = -R( 1 ) C C *** S = R * (LAMBDA2+S) CALL ADD(LMBDA2,S,S,K) CALL MLTPLY(S,R,S,K) C C *** IF N.LE.NMAX THEN RR(N-1)=R IF (N.GT.NX-1 )G0 TO 240 CALL COPY(R,RR( 1,N) ,K) C C *** LAMB0A1 = LAMB0A2 240 CALL C0PY(LMBDA2,LMBDA1 ,K) C C *** LAMBDA2 = 2*LAMBDA2-( ( 1 + ,5/N ) **2+ ( TAtJ/N ) **2 ) *L AMBDA/ ( 1+1/N) CALL MLTPLY(HALF,TEMP3,TEMP2fK) CALL AODlONE,TEMP2,TEMP2,K ) CALL MLTPLY(TFMP2,TFMP2,TEMP2,K ) CALL MLTPLY(TAU,TEMP3,TEMP,K) CALL MLTPLY( TEMP, TEMP, TEMP, K ) CALL ADD( TEMP2, TEMP, TEMP, K ) CALL MLTPLY{ TEMP, LMBOA, TEMP, K ) CALL A00(nME,TEMP3,TEMP3,K ) CALL 01 VIOF( TEMP, TEMP3, TEMP, K) 151 CALL MLTINT(LMBDA2t2fLMB0A2tt2 70,K) 270 CALL SBTRCTUMBDA2,TEMP,LMBDA2»K) C C *** LAMBDA = LAMBDA1 CALL C0PY(LMB0A1,LMB0A,K) C C *** N-N-l; IF N.GB.l THEN GO TO L2 N « N - 1 IF (N.GE.l)GO TO 210 C C *** P(0) « SUM/d + S) CALL ADD(ONE»S,TEMP t K) CALL DIVIDE(SUMtTEMPtP(l,l),K) C C *** FOR N=O t NMAX-l DO P ( N+ 1 ) = RR ( N > *P ( N ) IF (NX.EQ.l)GO TO 290 NXM1 = NX-1 DO 280 N=1,NXM1 CALL MLTPLY(RR(1,N) t P( l t N),P(l,N-H > t K) ?80 CONTINUE 290 DO 300 N=1,NX C C *** IF ( ABS( P(N)-PAPPROXIN) ) .GT .FPS I LON*ABS ( P ( N ) ) GO TO 320 CALL SBTRCT(P(1,N) »PAPRX( 1»N),TEMP,K) TEMP(l) = 1 CALL MLTPLY(EPStP( UN) ,TEMP2,K) TEMP2(1) = 1 CALL SBTRCT(TEMP,TEMP2,TEMP»K) IF (TEMP( 1) .GT.O)GO TO 320 300 CONTINUE C C *** T = 1 CALL COPY(ONEtT,K) C C *** FOR N=1,NMAX DO T=N*T; P(N)=T*P(N) IF (NX.EQ. 1 )RETURN DO 310 NN=2tNX N = NN-1 CALL MLT INT(T,N,T,£500 t K) 310 CALL HLTPLY, WITH INTEGER ORDER M AND C *** DEGREE (-1/2 ♦ N)t N«0 t 1 » • . . tNX-1 AND X.GT.l COMMON /LONG/LN,LW,LT,LS,LR,LTR,LEX,LRM INTEGER LN,LN,LT,LS,LEX REAL LR»LTR,LRM INTEGER X(1),M,NX,D,K,KX,Q(KX,1),QAPRX(KX,1),RR 20 HALF< 1 ) = 1 HALF(2) = -1 HALF(3) = 5000 DO 30 1=4, K 30 HALF! 11=0 C C *** CHECK IF INPUT PARAMETERS ARE IN CORRECT RANGE CALL SBTRCT(X,f)NE,TEMP,K) IF (TEMPI 1) .LE.O.OR.NX.LE.ORETURN 1 C C *** FOR N=0,NMAX DO QAPPROX(N)=0 CALL SET (TEMP, 0,840, K) 40 00 50 N=1,NX CALL COPY(TEMP,OAPRX( 1,N) ,K) 50 CONTINUE C C *** EPSILON = .5 * 10**(-D) EPS( 1 ) = 1 EPS(2) = -0/4 - 1 EPS(3) = 5 * 10**(-4*EPS( 2) -D -1) 00 60 1=4, K 60 EPS ( I ) = C C *** IF M.GE.O THEN FOR N=0,M-1 DO C=-(N+.5)*C IF (M.LT.O)GO TO 95 IF (M.E'J.O)GO TO 90 CALL C0PY(HALF,TEMP2,K ) DO 80 NN=1,M CALL MLTPLY(TEMP2,C,C,K) C(l ) = -C( 1) CALL ADD( TEMP2,0NE,TEMP2,K ) 80 CONTINUE 90 GO TO 120 156 no *** 120 *** IF M.LT.O FOR N=O t STEP -1»M+1 DO C > -C/1N-.5) 95 CALL C0PY(HALF,TEMP2,K) TEMP2(1) « -1 MM - -M DO 110 NN=1,MM CALL DIVIDE(C,TEMP2,CtK) C(l) » -C(l> CALL SBTRCT(TEMP2,ONE t TEMP2tK) CONTINUE SUM = C*( (X+l)/(X-l))** CALL DIVIDE(C,TEMP, SUM,K) CALL SET(TEMP,M,C300,K) CALL MLTPLY(TEMP, HALF, TEMP, K) CALL LEXP0(TEMP2, TEMP, TEMP, K) CALL MLTPLY(SUM, TEMP, SUM, K) *** XI = 2*X CALL ADD(X,X,X1,K) *** NU = 20 + ENTIER( ( 1.15*1 .0146*.00L22*M) / ( X-l ) )*NMAX) CALL SH0RTN(X,XS,£350,K) NU = 20 +AINT< (1.15+( .0146* .00122*M ) / { XS-1 ) ) * ( NX-1 ) ) NUINC = 10 *** LO: R = S = 140 CALL SET