Algorithm 467: Matrix Transposition in Place

Algorithms L.D. Fosdick and A.K. Cline, Editors Submittal of an algorithm for consideration for publication in Communications of the A C M implies unrestricted use of the algorithm within a computer is permissible. Copyright © 1973, Association for Computing Machinery, Inc. General permission to republish, but not for profit, all or part o f this material is granted provided that A C M ' s copyright notice is given and that reference is made to the publication, to its date o f issue, and to the fact that reprinting privileges were granted by permission of the Association for Computing Machinery. Algorithm procedure tqlrat (n,macheps) trans: (d, e2); value tl, macheps; integer ii; real macheps; array d, e2; comment n macheps d[l :tl] e2[l:n] d[1 :n] e2[l:n] Algorithm 464 Eigenvalues of a Real, Symmetric, Tridiagonal Matrix [F2] C h r i s t i a n H . R e i n s c h [ R e c d . 11 M a r . 1971] Mathematisches Institut der Technischen Universit~t, 8 0 0 0 M i i n c h e n 2, A r c i s s t r a 21, G e r m a n y Key Words and Phrases: eigenvalues, QR Algorithm CR Categories: 5.14 Language: Algol Description This algorithm uses a rational variant of the QR transformation with explicit shift for the computation of all of the eigenvalues of a real, symmetric, and tridiagonal matrix. Details are described in Ill. Procedures tredl or tred3 published in [2] may be used to reduce any real, symmetric matrix to tridiagonal form. Turn the matrix end-for-end if necessary to bring very large entries to the bottom right-hand corner. begin integer i, k, m; real b, b2,f, g, h, p2, r2, s2; for i : = 2 step 1 until n do e2[i--l] : = e2[i]; e2[n] : = b : = b2 : = f : = 0.0; for k : = 1 step 1 until n do begin h : = macheps X mac/wps X (d[k]~2 + e2[k]); if b2 < h t h e n begin b : = sqrt(h); b2 : = h end; comment Test for splitting; for m : = k step I until n do if e2[m] _< b2 then go to contl ; contl : if m = k then go to root; comment Form the shift from leading 2 X 2 block; nextit: g : = d[k]; p2 : = sqrt(e2[k]); h : = ( d [ k + l ] - - g ) / ( 2 . 0 X p 2 ) ; r2 : = sqrt(hXh+l.O); d[k] : = h : = p2/(ifh<O.O then h - - r 2 else h + r 2 ) ; h := g--h;f:=f+h; for i : = k d- 1 step 1 until n do d[i] : = d[i] -- h; comment Rational QL transformation, rows k through m; g : = d[m]; if g = 0.0 then g : = b; h : = g; s2 : = 0.0; f o r i : = m -- 1 step --1 u n t i l k d o begin p2 : = g X h; r2 : = p2 d- e2[i]; e2[i-F1] : = s2 X r2; s2 : = e2[i]/r2; d [ i + l ] : = tt + s2 X (h+d[i]); g : = d[i] -- e2[i]/g; if g = 0.0 then g := b; h : = g X p2/r2 end i; e2[k] : = s2 X g X h; d[k] : = /i; if e2[k] > b2 then go to next#; References 1. Reinsch, C.H. A stable, rational QR algorithm for the computation of the eigenvalues of an Hermitian, tridiagonal matrix. Math. Comp. 25 (1971), 591-597. 2. Martin, R.S., Reinseh, C.H., Wilkinson, J. H. Householder's tridiagonalization of a symmetric matrix. Numer. Math. 11 (1968), 181-195. Input: order of the matrix, the machine precision, i.e. minimum of all x such that 1 + x > 1 on the computer, represents the diagonal of the matrix, represents the squares o f the sub-diagonal entries, (e211] is arbitrary). Output: the computed eigenvalues are stored in this array in ascending sequence, is used as working storage and the original information stored in this array is lost; root: h:=d[k]+f; comment One eigenvalue found, sort eigenvalues; f o r i := kstep --1 until2do if h < d[i-- 1] then d[il : = d [ i - 1] else go to coot2; i:=1; cont2: d[i] : = h end k end tqlrat; 689 Communications of the A C M November 1973 Volume 16 Number 11 Algorithm 465 Student's t Frequency [S 14] G.W. Hill [Recd. 24 Aug. 1971, 23 Feb. 1972, 10 July 1972] C.S.I.R.O., Division of Mathematical Statistics, Glen Osmond, South Australia Key Words and Phrases: Student's t statistic, density function, series approximation CR Categories: 5.12, 5.5 Language: Algol Description T h e frequency function for Student's t distribution, r(½n + ½) f ( t in) - (~n)½r(½n~ (1 W t2/n) -(~'~+~), is evaluated for real t a n d real n > 0 to a precision near that o f the processor, even for large values of n. T h e factor involving t is evaluated as e x p ( - l ~ b ) where b is c o m p u t e d as (n W 1)In(1 W t2/n) if F / n = c is large ( > c m a x , say) or, to avoid loss of precision for smaller c, by s u m m i n g t h e series for b = (t 2 W c)(1 - c / 2 W c~/3 - c~/4 W " " ) until negligible t e r m s occur, i.e. cr/(r W 1) < ~, where ~ is the relative m a g n i t u d e of processor round-off. T h e relative error up to ~/ c m a x in evaluating hl(l W c) a n d the a c c u m u l a t e d round-off error of order E ,,/R in s u m m i n g a m a x i m u m of R t e r m s of the series c a n be limited to a b o u t the same low level by choosing c m a x = R -~ where R-½n/R "~ ~. T h u s for R = 12, 16, 23, or 32, values of c m a x ~.~ 0.2887, 0.25, 0.2085, or 0.1762, respectively, c o r r e s p o n d to processor precision where e = 2 -24, 2 -36, 2 -56, or 2 - u , respectively. Evaluation of the ratio of g a m m a functions by exponentiating the difference of almost equal values o f their logarithms would involve considerable loss of precision for large n. This is avoided by use of the asymptotic series obtained by differencing the Stirling a p p r o x i m a t i o n s , changing t h e variable to a = n - ½, a n d exponentiating the result (see also [1]): r(½n + ½) r(½n) (½a)~ ~ C~(4a) -~, ~=o where Co = C~ = 1, C2 = - - 1 9 / 2 , C3 = 631/2, C4 = - 1 7 4 3 1 7 / 8 , C~ = 204 91783/8, C6 = - 7 3 3 4 8 01895/16, C7 = 185 85901 54455/16, Cs = - 5 06774 10817 68765/128, C9 = 2236 25929 81667 88235/128, C10 = - 2 4 80926 53157 85763 70237/256. T h e relative error of t h e s u m of the first s terms is negligible for n > nmin where I c~ I × [4 ( n m i n - ½)]-2, ~ ~, e.g. for s = 5 a n d ~ = 2 -24 or 2 -3~, nmin ~ 6.271 or 13.76, respectively, a n d for s = 10 a n d ~ = 2 - ~ or 2 -~, nmht ~ 15.5 or 40.89, respectively. For smaller n the ratio of g a m m a functions is obtained f r o m the ratio for s o m e N >_ nmin by the relation: r(½n W ½) r(½n) n (n q- 2) (n W 1) (n + 3) (N - 2) r(½N W ½) (N-- 1) r(1N) F o r large n, processor u n d e r t o w at line 21 is avoided by use of the n o r m a l a p p r o x i m a t i o n , which is adequate for values of n > 1/~, whose representation is unaffected by subtraction o f 0.5. Protection against negative or zero n is provided by returning the distinctive value, - 1 . 0 , which m a y be s u p p l e m e n t e d by a n error diagnostic process, if required. For double precision calculations speed is improved by evaluating higher order t e r m s o f the g a m m a ratio series using single precision operations. C o m p a r i s o n o f double precision (~ = 2 -84) 690 results with single precision results (e = 2 -36, nmin = 13.76, c m a x = 0.25) for a Control D a t a 3200 indicated achievement generally o f a b o u t ten significant decimal digits, dropping to a b o u t eight significant decimals for a r g u m e n t s b e y o n d the 10-~0 probability level. Valuable c o m m e n t s f r o m the referee are gratefully acknowledged. Reference 1. Fields, J.L. A note on the asymptotic expansion of a ratio o f G a m m a functions. Proc. Edinburgh Math. Soc. Ser. 2 15 (1966), 43-45. Algorithm real procedure t frequency (t, n); value t, n; real t, n; if n _< 0.0 then t frequency : = - - 1 . 0 else begin real a, b, c, d, e, nmin, cmax; comment for 36-bit precision processor; nmin := 13.76; c m a x : = 0.25; b : = t >( t ; c : = b / n ; a := d : = b W c ; if c > c m a x then b : = (nW1.0) X h1(1.0-t-c) else for e : = 2.0, e W 1.0 while b ~ d d o begin a : = - - a )< c; b : = d; d : = a / e W d end; a : = n; c : = 0.3989422804; comment 1/sqrt(2r) = 0.3989422804014326779399461 . . . ; for e : = a while e < nrn'in do begin e : = c X a / ( a W l . O ) ; a : = a W 2.0 end; a : = a -- 0.5; if a ~ n then begin c : = sqrt(a/n) X c; a : = 0.25/a; a : = a )< a; c := ((((-21789.625XaW315.5) Xa-9.5) XaW1.0) XaW1.0) ×c end; t frequency := e x p ( - O . 5 X b ) X c end Student's t-frequency Algorithm 466 Four Combinatorial Algorithms [G6] Gideon Ehrlich [Recd. 25 Aug. 1971, 4 Jan. 1972, and 12 Dec. 1972] Department of Applied Mathematics, Weizmann Institute of Science, Rehovot, Israel Key Words and Phrases: permutations and combinations CR Categories: 5.39 Language: PL/I Description Each o f the following algorithms produce, by successive calls, a sequence of all combinatorial configurations, belonging to the appropriate type. PERMU P e r m u t a t i o n s of N>_3 objects: X(1), X ( 2 ) , . . . , X(N). C O M B I C o m b i n a t i o n s of M natural n u m b e r s o u t of the first N. COMPOMIN C o m p o s i t i o n s of a n integer P to M W 1 ordered terms, I N D E X ( k ) , each of which is not less t h a n a given minimum MIN(k). Communications of the A C M N o v e m b e r 1973 Volume 16 N u m b e r 11 C O M P O M A X The same as C O M P O M I N but each term has its own maximum M A X (k). The four algorithms have in common the important property that they use neither loops nor recursion; thus the time needed for producing a new configuration is unaffected by the "size" (N, N and M, P and M respectively) of that configuration. Each algorithm uses a single simple operation for producing a new configuration from the old one, that is: P E R M U A single transposition of two adjacent elements. COMBI Replacing a single element x by a y having the property that there is no element between x and y belonging to the combination. C O M P O M I N ( M A X ) Changing the values of two adjacent terms (usually only by 1). The algorithms are written in PL1 (F). Special instructions for the user and notes. P E R M U (1) The mean work-time is actually a decreasing function of N since, on (N -- 1)/N of the calls, it returns by the first RETURN. (2) The procedure operates directly on any object vector x[l:N]. (3) For the first permutation one must call FIRSTPER; for other permutations P E R M U must be used. (4) Together with the last permutation, which is the original one, we will get D O N E = 'I'B. If we continue to call P E R M U , the entire sequence will repeat indefinitely. If at any stage we set D O N E = 'O'B, then at the end of the appropriate sequence it will become 'I'B. (5) The entire resulting sequence is the same as that of Johnson [1] and Trotter [2]. C O M B I Every combination is represented in two forms: (1) As a bit array of M ' l ' s and N -- M 'O's which is identical to A(1), A ( 2 ) , . . . , A(N). (2) As an array C of M different integers not greater than N. The M elements are ordered according to their magnitude. If the second representation is not needed one can omit Z, H and C together with the last line of the procedure. For the first combination we can use the following initialization (for other initializations see [3]): D E C L A R E A(O:N) B I T (1), (X, Y, T(N), F(0:N), I, L, Z, H(N), C(M)) FIXED; DO K = 0 TO N -- M; A(K) = '0~B; END; DO K = N - - M + 1 TO N; A(K) = ' I ' B ; E N D ; DO K = 1 TO M; C(K) = N - - M - F K; H ( N - - M - F K) = K; END; T ( N - M) = - 1 ; T(1) = 0; F(N) = N - - M-}- 1 ; I = N - M; L = N; (The initialization was not done in the body of the procedure COMBI only in order to simplify the procedures C O M P O M I N MAX:.) Instead of using such a large number of parameters it is possible to retain only A, I, L as parameters of the procedure and declare and initialize the other present parameters in the body of the procedure (as is done in PERMU). In such a case N, T, F, L, H must be declared as S T A T I C or CONTROLLED ('own' in ALGOL). C O M P O M I N Each of the M + 1 MIN(k), as well as P, can be any integer (positive, negative, or zero), but the sum S of all those minima cannot be greater than P. For the first composition set INDEX(I) = P - S + MIN(I) INDEX(k) = MIN(k), for k > 1. Set N = P - S + M, and declare and initialize all variables that also appear in COMBI in the same way as was done for COMBI. Together with the last composition, we will get I = 0 as a signal to halt. C O M P O M A X The instructions for C O M P O M I N are valid for C O M P O M A X provided: (1) M I N is replaced by M A X (S > P); and (2) N is initialized to N = S - P -t- M. The vector C (but not H!) has no use in COMPOMIN(M,4X), so one can omit all statements in which it appears. A justification for the four algorithms and for some others can be found in [3]. References 1. Johnson, S.N. Generation of permutations by adjacent transformations. Math. Comp. 17 (1963), 282-285. 2. Trotter, H.F. Algorithm 115, Perm. Comm A C M 5 (Aug. 1962), pp. 434-435. 3. Ehrlicb, G., Loopless algorithms for generation permutations combinations and other combinatorial configurations. J. A C M 20 (July 1973), 500-513. Algorithm FIRSTPER: PROCEDURE (X,DONE); DECLARE (X(*), (XN,XX) STATIC) DECIMAL, DONEBIT(1) (N,S,V,M,LjI,DI,IPI) BINARYSTATIC, (p(O:N),IP(N-I~D(N-I),T(N)) BINARYCONTROLLED; N=OIM(X,I); IF ALLOCATION(P) THENFREEP,IP,D,T; ALLOCATEP,IP,D,T; DOM=l TO N-l; P(M),IP(M)=M; D(M)=-]; END; XN=X(N); V=-I; S,P(O),P(N)=N; M,L=I; T(N)=N-I; T(N-I)=-2; T(2)=2; DONE='O'B; PERMU: ENTRY (X,DONE); IF S.=M THENDO; X(S)=X(S+V); S=S+V; X(S)=XN; RETURN;END; I=T(N); DI=D(1); IP(1),IPI=IP(1)+DI; M=P(IPI); IP(M)=IPI-OI; P(IPI-DI)=M; P(IPI)=I; M=IPI+L; XX=X(M); X(M)=X(M-OI); X(M-DI)=XX; L=I-L; V=-V; M=N+I-S; IF P(IPI+DI) < I THEN DO; IF I=N-l THENRETURN; T(N)=N-I; T(N-I) = -I; RETURN; END; D(1)=-DI; IF T(1) < 0 THEN DO; IF T(1)-~=l-I THENT(I-I)=T(1); T(1)=I-l; END; IF I ~ = N - l THENDO; T(N)=N-I; T(N-l)=-I-l; END; T(I+I)=T(1); IF I=2 & P(2)=2 THENDONE=')'B; END; COMBI PROCEDURE (A,N,X,Y,T,F,I,L,Z,H,C); DECLARE A(*)BIT(1), (N,X,Y,T(*),F(*),I,L,Z,H(*),C(*)) FIXED; IF T(1) < 0 THEN DO; IF-T(I)"~=I-] THEN T(I-I)=T(I); T(1)=I-]; END; IF--, A(1) THEN DO; X=I; Y=F(L); IF A(I-7) THEN F(I)=F(I-I); ELSEF(I)=I; IF F(L)=L THEN DO; L=I; I=T(1); GOTOCHANGE;END; IF L=N THEN DO; T(F(N))=-I-]; T(I+I)=T(1); I=F(N); F(N)=F(N)+I; GOTOCHANGE; END; T(L)=-I-I; T(I+])=T(1); F(L)=F(L)+~; I=L; GOTOCHANGE END; Y=I; IF I-I=L THEN DO; F(L),X:F(L)-I; F(I-I)=P(1); IF L=N THEN DO; IF I=F(N) -l THENDO; I=T(1); GOTOCHANGE;END; T(F(N)-])=-I-I; T(I+I)=T(1); I=F(N)-l; GOTOCHANGE; END; T(L)=-I-I; T(I+I)=T(I); I=L; GOTOCHANGE; END; X=N; F(L-I)=F(L); F(N)=N; L=N; IF I=N-I THENDO; I=T(N-]); GOTOCHANGE;END; T(N-I)=-I-I; T(I+I)=T(1); I=N-l; CHANGE; A(X)='I'B; A(y)='O'B; H(X),Z=H(Y); C(Z)=X; END COMBI; COMPOMIN: PROCEDURE(INDEX,A,N,X,Y,T,F,I,L,Z,H,C); DECLAREA(*) BIT(1), (INDEX(*),N,X,Y,T (*),F(*),I,L,Z,H(*),C(*)) FIXED; CALL COMBI (A,N,X,Y,T,F,I,L,Z,H,C); INDEX(Z)=INDEX(Z)+X-Y; INDEX(Z+I)=INDEX(Z+I)~Y-X: END COMPOMIN; COMPOMAX:PROCEDURE(INDEX,A,N,X.Y,T,F.I,L,Z.H,C); DECLAREA(*) BIT(1), (INDEX(*),N,X,Y,T(*),F(*),I,L,Z,H(*),C(*)) FIXED; CALL COMBI (A,N,X,Y,T,F,I,L,Z.H,C); INOEX(Z)=INDEX(Z)-X+y; INDEX(Z+I)=INDEX(Z+]).y+x; END COMPOMAX; Acknowledgment. I would like to thank Professor Shimon Even for guidance and encouragement. 691 Communications of the ACM November 1973 Volume 16 Number 11 Algorithm 467 Pr~oor. Both m and s are divisible by d, and therefore so is any subcycle element sn r (mod m). But n and m have no common factors (since m = nn2 - 1), so no divisor of m larger than d can divide sn r. [] THEOREM 2. For every subcycle beginning with s, there is Matrix Transposition in Place [F 1] another (possibly the same) subcycle beginning with m - N o r m a n B r e n n e r [ R e c d . 14 F e b . 1972, 2 A u g . 1972] M.I.T., Department of Earth and Planetary Sciences, Cambridge, MA 02139 s. PROOV. The elements of the second subcycle are just - s n r (rood m). It is the same subcycle if for some r, n r --- - 1 (mod m'), for m' = m / ( s , m). [] The next theorem gives the group representation of the integers modulo m. THEOREM 3. Factor m into powers o f primes, m = p~1 . . . pt,~l. Let ri be a primitive root o f pi ; that is, the powers ri k (mod Pi) f o r k = 0, 1 , . . . , p - 2, comprise every positive integer less than Pi • Define the generator gl = 1 -b R m / p i etl , where R =-- (ri -- 1) (m/p~. i) -1 (mod p~. i). Define the Euler totient function 4~(1) = 1; Key Words and Phrases: transposition, matrix operations, permutations, primitive roots, number theory CR Categories: 3.15, 5.14, 5.39 Language: Fortran otherwise 4~(k) = the number o f integers less than k having no common factor with it. Then, f o r any integer x less than m, there exist unique indices ji f o r which 0 <__ ji < 4~(p~i/(x, p ~ i ) ) and Description Introduction. Since the problem of transposing a rectangular matrix in place was first proposed by Windley in 1959 [I], several algorithms have been used for its solution [2, 3, 7]. A significantly faster algorithm, based on a number theoretical analysis, is described and compared experimentally with existing algorithms. Theory. A matrix a, of n~ rows and n~. columns, may be stored in a vector v in one of two ways. Element a~j (0-origin subscripts) may be placed rowwise at vk, k = i m + j, or columnwise at vk,, k' = i + jn~ . Clearly, letting n = ill and m = n~n~. - 1, k' =-- nk (mod m). (1) Transposition of the matrix is its conversion from one mode of storage to the other, by performing the permutation (1). This permutation may be done with a minimum of working storage in a minimum number of exchanges by breaking it into its subcycles. For example, for a 4 X 9 matrix, one subcycle representation is (0) (1 4 16 29 11 9) (34 31 19 6 24 26) (22 18 2 8 32 23) (13 17 33 27 3 12) (5 20 10) (30 15 25) (7 28) (14 21) (35). The notation for the sixth subcycle, for example, means that V5 4 - - - V 2 0 ~ Vl0 ~ - " ~)5 • For a subcycle starting with element s, the elements of the subcycle are sn ~ (mod m), for r = 0, 1 , . . . . The following theorems are easily established. THEOREM 1. All the elements o f the subcycle begbming with s are divisible by d = (s, m), the largest common factor o f both s and m. They are divisible by no larger divisor o f m. x --- (x, m)gil I . . . g~Z(modm). PROOf. In [4]; if any p~ = 2, replace g ~ by 4-5 jl, where 0 _< j~ < 4~(2'~1-2/(x, 2'~1-2)). [] For example, for m = 35, as in our example above, x ~22h31 j2 (mod 35) for (x,35) = 1 and for 0 < j~ < 4 and 0 < j2 < 6. Index notation is analogous to logarithmic notation in that multiplication modulo m becomes merely addition of indices. The following theorem solves the problem of the subcycle starting points. It is similar to the algorithm in [6]. THEOREM 4. Let n and m be defined as f o r (1). Then, f o r any hlteger x less than m, upper bounds Ji m a y be f o u n d so that unique flMices ji exist in the range 0 <_ .]i < Ji and x =-- + ( x , m ) , / O g ~ l . . . g~t (mod m). PROOF. Express n and - 1 in index notation. Then, compute from the indices of n the smallest e such that n ~ --- l (mod m). Initially, set each Ji = ¢ ( p ~ i / ( x , pTi)). Next, doing only index arithmetic, examine each power -4-nJ for nontrivial relations of the form g ~ ~ =tznig~ 1 . . . g~l (mod m / ( x , m)) where 0 _< jk < Jk for each k. Then set Ji = Ji • Stop when the product of the Ji and e equals 4~(m/(x, m)), which is the number of integers in subcycles divisible only by (x, m). [] Notice that the choice of J; by this method is not unique. For example, continuing from above, for (x, m) = 7, n = 4, x --- 7.4/°22 h (mod35), for 0 _<j0 < 2 a n d 0 < j ~ < 2. The relations found were ( - 1)1 --- 41 (mod 5), 222 ~- 41 (mod 5) and 311 ~ 4o (mod 5). Theorem 4 is more important in theory than in practice. The Timing Tests Hi 112 m (all times in msec) Alg.302 Alg.380 45 45 46 46 47 47 45 45 46 46 47 47 692 50 60 50 60 50 60 180 200 180 200 180 200 13.173 2699 112.19 31.89 34-29 2819 7.13.89 8999 17.487 9199 11.769 3.13.241 Alg.380 T~ T~ 1WRK= (nlWn2)/2 T8 350 558 367 425 383 483 1200 1767 1816 1700 1450 983 317 123 339 350 378 127 1050 408 1233 508 1133 1150 167 117 217 250 267 133 816 416 583 417 667 1067 IWRK=O XPOS NWORK=O Ti XPOS NWORK= (n1+n2)/2 T5 T1/ T4 T2/ T4 T3/ T5 133 90 106 133 72 90 517 283 267 383 383 550 67 100 83 83 67 100 300 300 267 317 267 467 2,62 6,20 3,46 3,19 5,18 5,36 2,25 6,25 6,41 4,44 3,78 1~69 2,38 1,37 3,21 2,63 5,23 1,41 2,03 1,44 4,63 1,33 2,96 2,09 2,50 1,17 2,60 3,00 4,00 1,33 2,72 1,39 2,19 1,32 2,50 2,29 Communications of the ACM November 1973 Volume 16 Number 11 tremendous labor in finding primitive roots for large primes (since a table of roots is very bulky) and in finding the index representation of n is not compensated for by time savings afterward; see the timing tests below. The same practical objection holds against the algorithm in [6]. Algorithm. An efficient program breaks naturally into two parts. First determine starting points for the subcycles and then move the data. In each part, the program below is significantly faster than Algorithm 380 in [3]. For each divisor d of m, the subcycles beginning with d and with m - d are done. If the number of data moved is still less than @(m/d), further subcycle starting points of the form sd are tried, for s = 2, 3, . . . . The most general test is that sd is acceptable if no element in its subcycle is less than sd or greater than m - sd. Since this test requires much time-consuming computation, it is much faster to look for sd in a table where marks are made to indicate that an element has been moved. In some applications, a bit within each datum may be used. For example, if the data are all biased positive, the sign bit may be used; or, for normalized, nonzero, binary floating point data, the high bit of the fraction is always one and so may be used. In general, a special table of length N W O R K is used. As in [3], N W O R K = (n~ q- m ) / 2 was found to be sufficient for most cases. However, when m has many divisors, Algorithm 380 must perform the time-consuming general test for many possible starting points when the new algorithm need not. The inner loop of the algorithm computes (1), moves data, marks in the table, and checks for loop closure. Since the major part of the time of the inner loop is calculating (1), time is saved over Algorithm 380 by moving elements vk and v,~-k simultaneously. In special cases, further savings may be made. For example, m is divisible by 2 only when both n~ and nz are odd. Then the subcycles beginning at m / 2 - s and m / 2 + s may be done simultaneously with the subcycles from s and m - s, thus reducing the number of times (1) is computed. Timb~g tests. A set of test matrices were transposed on the 360/65 with all programs written in Fortran H, OPT = 2. The new algorithm was always faster than both Algorithm 380 [3] and Algorithm 302 [2] when N W O R K = (m + n~)/2. When N W O R K = O, it was slower than Algorithm 380 (for 1 W R K = 0) and Algorithm 302 only for a few cases when n~n~ < 100. It was especially faster than Algorithm 380 when m = thn2 - 1 had many factors and there were hence many subcycles. An experiment was made for cases when m was prime. A known primitive root of m was then taken from a table [5] and was used to generate subcycle starting points. Since no time was wasted in finding the primitive root or in finding subcycle starting points, this test showed the maximum time savable by implementing Theorem 4. For N W O R K = (m + n2)/2 and m > 200, no improvement was found over the normal algorithm. For N W O R K = O, the gain in speed was never more than 25 percent. References 1. Windley, P.F. Transposing matrices in a digital computer. Comp. J. 2 (Apr. 1959), 47-48. 2. Boothroyd, J. Algorithm 302, Transpose vector stored array. Comm. A C M 10 (May 1967), 292-293. 3. Laflin, S., and Brebner, M.A. Algorithm 380: In-situ transposition of a rectangular matrix. Comm. A C M 13 (May 1970), 324-326. 4. Bolker, E. An bttroduction to Number Theory: An Algebraic Approach. Benjamin, New York, 1970. 5. Abramowitz, M., and Stegun, I. Handbook (~fMathematical Functions, Table 24.8. Nat. Bur. of Standards, Washington, D.C., 1964. 6. Pall, G., and Seiden, E. A problem in Abelian Groups, with application to the transposition of a matrix on an electronic computer. Math. Comp. 14 (1960), 189-192. 7. Knuth, D., The Art o f Computer Programming, Vol. I. Addison-Wesley, Reading, Mass., 1967, p. 180, prob. 12, and p. 517, solution to prob. 12. 693 Algorithm SUBROUTINE XPOSE(A, N I , N2, H I 2 , MOVED, NWGRK) C TRANSPOSITION OF A RECTANGULAR MATriX IN S I T U . C BY NORMAN BRENNER, MIT, 1 / 7 2 . CF. ALG. 3 8 0 , CACM, 5 / 7 0 . C TRANSPOSITION OF THE NI BY N2 MAT~IX A AMOUNTS TO C REPLACING THE ELEMENT AT VECTOR POSITION I ( O - O R I G I N ) C WITH THE ELEMENT AT POSITION N I * I (MOD H I * N O - 1 ) . C EACH SUBCYCLE OF THIS PERMUTATION IS COMPLETED IN ORDEr. C M O V E D IS A L O G I C A L WORK ARRAY OF LENGTH NWO~K. LOGICAL MOVED DIMENSION A(NI2), MOVED(NWO~K) C REALLY A(NI,NB), BUT NI2 = NI*N2 DIMENSION I F A C T ( 8 ) , IPQWEH(8), N E X P ( 8 ) , I E X P ( 8 ) IF (NI.LT.2 .OR. N 2 . L T . 2 ) RETURN N = NI M = NI*NB - 1 IF (NI.NE.NB) G O TO 3 0 C SQUARE MATRICES ARE DONE SEPARATELY FOR SPEED IIMIN = 2 DO O0 IIMAX=N,M,N 12 = I I M I N + N - I DO 1 0 I I = I I M I N , IIMAX ATEMP = A { I I ) A(II) = A(I2) A ( I 2 ) = ATEMP I2 = I2 + N IO CONTINUE IIMIN = IIMIN + N + I 20 CONTINUE RETURN C MODULUS M IS FACTORED INTO PRIME POWERS. EIGHT FACTORS C S U F F I C E UP T0 M = 2 . 3 . 5 . 7 " 1 1 " 1 3 " 1 7 " 1 9 = 9,767,520. 30 CALL FACTOR(M, IFACT, IPOWER, NEXP, NPOWER) D O 40 I P = I , N P O W E R IEXP(IP) = 0 40 C O N T I N U E C GENERATE EVERY DIVISOR OF M LESS THAN M/2 IDIV = ! 50 I F ( I D I V . G E . M / 2 ) GO TO 190 C THE NUMBER OF ELEMENTS WHOSE INDEX 15 D I V I S I B L E BY I D I V C AND BY NO OTHE~ DIVISOR OF M IS THE EULER TOTIENT C FUNCTION, P H I ( M / I D I V ) . NCOUNT : M / I D I V DO 60 I P = I , N P O W E R IF (IEXP(IP).EQ.NEXP(IP)) GO TO 60 NCOUNT = (NCOUNT/IFACT(IP))*(IFACT(IP)-;) 60 C O N T I N U E DO 70 I=I,NWO~K MOVED(I) = .FALSE. 7 0 CONTINUE C THE STARTING POINT OF A SUBCYCLE IS D I V I S I B L E ONLY BY I D I V C AND MUST NOT APPEAR IN ANY OTHER SUBCYCLE. ISTART = I D I V 80 MMIST = M - ISTART I F ( I S T A R T . E Q . I D I V ) GO TO IBO I F (ISTART.GT.NWOaK) GO TO 90 I F (MOVED(15TAKT)) GO TO 160 90 IS OID = I S T A R T / I D I V DO 1 0 0 I P = I , N P O W E R I F ( I E X P ( I P ) . E Q . N E X P ( I P ) ) GO TO l O G I F (MOD (IS OID , I F A C T ( I P ) ) . E Q . O ) GO TO 160 100 CONTINUE I F (ISTART.LE.NWORK) GO TO 120 I T E S T = ISTART 110 ITEST = MOD(N*ITEST,M) IF (ITEST.LT.ISTART .OR. ITEST.GT.MMIST) G O TO 1 6 0 I F ( I T E S T . G T . I S T A R T ,AND. I T E S T . L T . M M I S T ) GO T O 110 120 ATEMP = A ( I S T A R T + I ) BTEMP = A ( M M I S T + I ) I A I = ISTART 13U I A 2 = M O D ( N * I A I , M ) MMIAI = M - I A I MMIA2 = M - I A 2 I F ( I A I . L E . N W O R K ) M O V E D ( I A I ) = .TRUE. I F (MMIA I.L E .N W OR R ) MOVED(MMIAI) = .TRUE. NCQUNT = NCOUNT - 2 C MOVE TWO E L E M E N T S , THE SECOND FROM THE N E G A T I V E C SUBCYCLE. CHECK FIRST FOR S U B C Y C L E C L O S U R E . I F ( I A 2 . E Q . I S T A R T ) GO TO 140 IF (MMIA2.EQ.ISTART) GO TO 1 5 0 A(IAI+I) = A(IA2+I) A I M M I A I ÷ I ) = A(MMIAB+I) I A I = IA2 GO TO 1 3 0 = ATEMP 140 A(IAI+I) A(MMIAI*I) = BTEMP GO T 0 1 6 0 : BTEMP 150 A(IAI+I) A(MMIAI÷I) = ATEMP 160 ISTART = ISTART + IDIV I F (NCOUNT.GT.O) GO T O 80 DO 1 8 0 IP=I,NPOWER IF (IEXP(IP).EQ.NEXP(IP)) GO T O 1 7 0 IEXP(IP) = IEXP(IPI + I IDIV = IDIV*IFACT(IP) GO T O 50 170 IEXP(IP) = 0 IDIV = IDIV/IPOWER(IP) • Communications of tile ACM November 1973 Volume 16 Number 11 180 190 C C C CONTINUE RETURN END SUBROUTINE FACTOR(N, ~FACT,, I P O W E R s NEXP,,, NPOWER) FACTOR N INTO ITS PRIME POWENSs NPOWER IN NUMBER. E.G., FOR N = 1 9 6 0 = 2 . . 3 .5 .7..2, NPOWER=3* I F A C T = 3 , S , TJ I P O W E R : B , S J A g s AND N E X P = 3 , 1 , 2 . DIMENSION I F A C T ( B ) , IPOWER(S)J N E X P ( 8 ) IP = 0 IFCUR = 0 NPART = N IDIV = 2 10 IOUOT = N P A R T / I D I V IF" ( N P A R T - I D I V * I O U O T ) 60~ 20s 60 20 IF (IDIV-IFCUR) 40, 40~ 30 3 0 IP : IP + 1 IFACT(IP) : IDIV IPOWER(IP) = I D I V IFCUR = I D I V NEXP(IP) = I GO TO 5 0 40 IPOWER(IP) = IDIV*IPOWER(IP) NEXP(IP) = NEXP(IP) + 1 50 NPART = IQUOT GO TO 60 70 IF IF 10 (IQUOT-IDIV) (IDIV-2) SO, BO I D I V = 3 GO TO I 0 90 I D I V = I D I V I00 II0 120 140 70 + 2 GO TO 1 0 IF (NPART-1) lzlO, 140, 110 I F ( N P A R T - I F C U R ) 130, 130, 120 IP = IP + I IFACT(IP) = NPART IPOWER(IP) = NPART NEXP(IP) = I GO TO 130 t00, lOOJ 80, 90 140 IPOWER(IP) = NPART*IPOWER(IP) NEXP(IP) = NEXP(IP) + I NPOWER = I P RETURN END Algorithm 468 Algorithm for Automatic Numerical Integration Over a Finite Interval [D 1] T . N . L . P a t t e r s o n [ R e c d . 20 J a n . 1971, 27 N o v . 1972, 12 D e c . 1972, 26 M a r . 1973] Department of Applied Mathematics and Theoretical Physics, The Queen's University of Belfast, Belfast BT7 INN Northern Ireland Key Words and Phrases: automatic integration, numerical integration, automatic quadrature, numerical quadrature CR Categories: 5.16 Language: Fortran Editor's note: Algorithm 468 described here & available on magnetic tape from the Department o f Computer Science, University o f Colorado, Boulder, CO 80302. The cost .for the tape is $16.00 (U.S. and Canada) or $18.00 (elsewhere). I f the user sends a small tape (wt. less than 1 lb.) thealgorithm will be copied on it and returned to him at a charge o f $10.O0 (U.S. only). All orders are to be prepaid with checks payable to A C M Algorithms. The algorithm is recorded as one file of BCD 80 character card images at 556 B.P.I., even parity, on seven track tape. We will supply algorithm at a density of 800 B.P.I. i f requested. Cards .for algorithms are sequenced starting at 10 and brcremented by 10. The sequence number is right justified in column 80. Although we will make every attempt to insure that the algorithm cotffbrms to the description printed here, we cannot guarantee it, nor can we guarantee that the algorithm is correct.--L.D.F. and A.K.C. 694 Description Purpose. The algorithm attempts to calculate automatically the integral of F(x) over the finite interval [A, B] with relative error not exceeding a specified value e. Method. The method uses a basic integration algorithm applied under the control of algorithms which invoke, if necessary, adaptive or nonadaptive subdivision of the range of integration. The basic algorithm is sufficiently powerful that the subdivision processes will normally only be required on very difficult integrals and might be regarded as a rescue operation. The Basic Algorithm. The basic algorithm, Q UAD, uses a family of interlacing whole-interval, common-point, quadrature formulas. The construction of the family is described in detail in [1]. Beginning with the 3-point Gauss rule, a new 7-point rule is derived, with three of the abscissae coinciding with the original Gauss abscissae; the remaining four are chosen so as to give the greatest possible increase in polynomial integrating degree; the resulting 7-point rule has degree 11. The procedure is repeated, adding eight new abscissae to the 7-point rule to produce a 15-point rule of degree 23. Continuing, rules using 31, 63, 127, and 255 points of respective degree 47, 95, 191, and 383 are derived. The 255-point rule has not previously been published. In addition, a 1-point rule (abscissa at the mid-point of the interval of integration) is included in the family to make eight members in all. The 3-point Gauss rule is in fact formally the extension of this 1-point rule. The successive application of these rules, until the two most recent results differ relatively by E or better, is the basis of the method. Due to their interlacing form, no integral evaluations need to be wasted in passing from one rule to the next. The algorithm has been used for some time on practical problems and has been found to generally perform reliably and efficiently. Its domain of applicability generally coincides with that of the Gauss formula, which is much wider than commonly supposed [2]. It will perform best on " s m o o t h " functions, but the degree of deterioration of performance when applied to functions with various types of eccentricities depends more on the harshness of these eccentricities than on their presence as such. Integrands with large peaks or even singularities at the ends of the interval of integration are handled reasonably well. It may be noted that none of the rules actually uses the end points of the interval as abscissae. Peaks in the integrand at the center of the interval and discontinuities in the integrand are less easily dealt with. Although it is recommended that the algorithm be applied using the control algorithms described later, if desired it can be used directly as follows. The algorithm is entered by the statement: C A L L QUAD (A, B, RESULT, K, EPSIL, NPTS, ICHECK, F) The user supplies: A lower limit of integration. B upper limit of integration. E P S I L required relative error. F F(X) is a user written function to calculate the integrand. The algorithm returns: R E S U L T an array whose successive elements R E S U L T ( l ) , RESULT(2), etc., contain the results of applying the successive members of the family of rules. The number of rules actually applied depends on EPSIL. The array should be declared by the calling program to have at least eight elements. K element, RESULT(K), of array R E S U L T contains the value of the integral to the required relative accuracy. K is determined from the convergence criterion: [ R E S U L T (K) -- R E S U L T (K -- 1) ] <_ EPSIL* I R E S U L T (K) I N P T S number of integrand evaluations. ICHECK this flag will normally be 0 on exiting from the subroutine. However, if the convergence criterion above is not satisfied after exhausting all members of the family of rules, then the flag is set to 1. Communications of the ACM November 1973 Volume 16 Number 11 Table I. Test Integrals and Their Values .,/ x dx 1. = 2 o 2. [0.92 cosh (x) - cos (x)] dx - 0.4794282267 3. d x / ( x 4 d- x ~ -k- 0.9) - 1.582232964 4. fl 5. foI d x / ( l + x ~) - 0.8669729873 6. fo1 dx/(1 + 0.5 sin (31.4159x)) - 1.154700669 2 x~ dx = 5 1 7. o x d x / ( e ~ -- 1) -- 0.7775046341 1 8. sin (314.159x)/(3.14159x) dx - 0.009098645256 fo .1 9. o ° 50 dx/(2500x 2 -t- 1)/3.14159 - 0.4993638029 3,141,5927 10. f cos (cos (x) -k- 3 sin (x) + 2 cos (2x) "10 + 3cos (3x) + 3 sin (2x)) dx "- 0.8386763234 11. ox In (x) dx = - 1 . 0 12. o 4~'2x sin (20~rx) cos (27rx) dx - 1 -0.6346651825 1 13. fo d x / ( l q- (230x - 30) 5) - 0.0013492485650 The control algorithms. Two control algorithms are provided, QSUBA and QSUB, which if necessary invoke subdivision respectively in either an adaptive or a nonadaptive manner. QSUBA is generally more efficient than QSUB, but since there are reasons for believing [2] that adaptive subdivision is intrinsically less reliable than the nonadaptive form, an alternative is provided. The adaptive algorithm QSUBA. QUAD is first applied to the whole interval. If a converged result is not obtained (that is, the convergence criterion is not satisfied), the following adaptive subdivision strategy is invoked. At each stage of the process an interval is presented for subdivision (initially the whole interval (A, B)). The interval is halved, and QUAD applied to each subinterval. If QUAD fails to converge on the first subinterval, the subinterval is stacked for future subdivision and the second subinterval immediately examined. If QUAD fails to converge on the second subinterval, it is immediately subdivided and the whole process repeated. Each time a converged result is obtained it is accumulated as the partial value of the integral. When QUAD converges on both subintervals the interval last stacked is chosen next for subdivision and the process repeated. A subinterval is not examined again once a converged result is obtained for it, so that a spurious convergence is more likely to slip through than for the nonadaptive algorithm QSUB. The convergence criterion is slightly relaxed in that a panel is deemed to have been successfully integrated if either QUAD converges or the estimated absolute error committed on this panel does not exceed ~ times the estimated absolute value of the integral over (A, B). This relaxation is to try to take account of a common situation where one particular panel causes special difficulty, perhaps due to a singularity of some type. In this case, QUAD could 695 obtain nearly exact answers on all other panels, and so the relative error for the total integration would be almost entirely due to the delinquent panel. Without this condition the computation might continue despite the requested relative error being achieved. The risk of underestimating the relative error is increased by this procedure and a warning is provided when it is used. The algorithm is written as a function with value that of the integral. The call takes the form: QSUBA(A, B, EPSIL, N P T S , ICHECK, RELERR, F) and causes F(x) to be integrated over (A, B) with relative error hopefully not exceeding EPSIL. R E L E R R gives a crude estimate of the actual relative error obtained by summing the absolute values of the errors produced by QUAD on each panel (estimated as the differences of the last two iterates of QUAD) and dividing by the calculated value of the integral. The reliability of the algorithm will decrease for large EPSIL. It is recommended that E P S I L should generally be less than about 0.001. F should be declared E X T E R N A L in the calling program. N P T S is the number of integrand evaluations used. The outcome of the integration is indicated by ICHECK: ICHECK = 0. Convergence obtained without invoking subdivision. This corresponds to the direct use of QUAD. I C H E C K = 1. Subdivision invoked and a converged result obtained. ICHECK = 2. Subdivision invoked and a converged result obtained but at some point the relaxed convergence criterion was used. If confidence in the result needs bolstering, E P S I L and R E L E R R may be checked for a serious discrepancy. ICHECK negative. If during the subdivision process the stack of delinquent intervals becomes full a result is obtained, which may be unreliable, by continuing the integration and ignoring convergence failures of QUAD which cannot be accommodated on the stack. This occurrence is noted by returning I C H E C K with negative sign. The nonadaptive algorithm QSUB. QUAD is first applied to the whole interval. If a converged result is not obtained the following nonadaptive subdivision strategy is invoked. Let the interval (A, B) be divided into 2N panels at step N of the subdivision process. QUAD is first applied to the subdivided interval on which it last failed to converge, and if convergence is now achieved, the remaining panels are integrated. Should a convergence failure occur on any panel, the integration at that point is terminated and the procedure repeated with N increased by one. The strategy insures that possibly delinquent intervals are examined before work, which later might have to be discarded, is invested on well behaved panels. The process is complete when no convergence failure occurs on any panel, and the sum of the results obtained by QUAD on each panel is taken as the value of the integral. The process is very cautious in that the subdivision of the interval (A, B) is uniform the fineness of which is controlled by the success of QUAD. In this way it is much more difficult for a spurious convergence to slip through than for QSUBA. The convergence criterion is relaxed as described for QSUBA. The algorithm is used in the same way as QSUBA and is called with the same arguments as QSUBA. One of the possible values of I C H E C K has a different interpretation: I C H E C K negative. If during the subdivision process the upper limit on the number of panels which may be generated is reached, a result is obtained, which may be unreliable, by continuing the integration ignoring convergence failures of Q UAD. This occurrence is noted by returning I C H E C K with negative sign. Tests. The algorithms have been found to perform reliably on a large number of practical problems. To give a feeling for the performance, results for a number of contrived examples are given using the adaptive control algorithm, QSUBA. It would be difficult to justify these examples as acid tests of any method, but they have the advantage of having being quoted at various times in the literature. For comparison a number of automatic procedures were used, which include S Q U A N K [3] (adaptive Simpson), as well as the Communications of the ACM November 1973 Volume 16 Number 11 Table II. Relative Error Requested, I0 -~ Integral NCADRE NQSUBA TeADRE/TQsvBa 1 17 15 1.8 2 17 7 2.9 3 33 15 4.4 4 9 7 1.9 5 9 7 2.2 6 175 127 3.2 7 9 7 1.8 8 1137 255 8.5 9 97 127 2.4 10 107 63 2.2 11 137 31 9.9 12 252 63 6.3 13 129 787 .52 N and T with a p p r o p r i a t e subscripts give respectively the number of integrand evaluations and the time taken for the computation. Table III. Relative Error Requested, 10-6 1 2 3 4 5 6 7 8 9 10 11 12 13 33 33 49 129 17 401 9 2633 28 l 193 233 532 305 63 15 3l 31 15 255 7 255 255 63 795 127 1001 .75 2.6 3.0 5.0 2.0 2.9 1.8 18. 2.4 3.8 .74 6.4 .90 Table IV. Relative Error Requested, 10-8 1 2 3 4 5 6 7 8 9 10 11 12 13 65 33 97 545 65 569 17 4001 337 305 297 932 48l 255 15 31 31 31 255 15 255 255 127 2415 127 1017 .36 2.7 4.9 20. 3.6 3.8 1.6 24. 2.8 2.8 .28 10. 1.1 modified Havie integrator [4] and CADRE [5] (both based on the Romberg scheme). The latter algorithm, which attempts to detect certain types of singularities using the Romberg table, was found, on the examples tried, to be the best overall competitor to QSUBA, and only this comparison is quoted. The Havie algorithm was particularly poor and had the disturbing feature of converging spuriously on periodic integrands. Thacher [6] has described the shortcomings of Romberg integration, and Algorithm 400 appears to exhibit them. SQUANK was found, to be quite good when used at low accuracy, but the performance deteriorated as the demand for accuracy increased. It also gave trouble on some of the more awkward integrals such as 8 and 11. SQUANK also computes the integral in the context of absolute error, and since this is meaningless unless an estimate of the order of magnitude of the integral is known, the algorithm can hardly be described as automatic. CADRE allows a choice of absolute or relative error. A criticism sometimes levied at relative error is that should the integral turn 696 out to be zero a difficulty will arise. The only advice that can be offered in this respect is that, should a user suspect that this is likely to happen, a constant should be added to the integrand reflecting some appropriate quantity such as the maximum of the integrand. The constant which will be integrated exactly can be removed after the algorithm has done its work. The test integrals are listed in Table I, and the results obtained for various required relative accuracies in Tables II, II1, and IV. Generally QSUBA is superior by a substantial margin. The methods are compared in terms of the number of integrand evaluations needed to obtain the required accuracy and also in terms of the times required. For simple integrands the bookkeeping time of some methods can be significant, and QUAD can obtain a considerable advantage by its relative simplicity. Integrals 11 and 13 are interesting examples of this. The number of integrand evaluations exceeding 255 indicates that QSUBA invoked subdivision to obtain the result. In Tables III and IV QSUBA returned ICHECK = 2 on integral 11, but the requested tolerance was achieved. Integral 8 caused special difficulty to CADRE, and for Tables III and IV a converged result could be obtained only after a relatively large investment of computer time. The feature of CADRE to detect certain singularities should show up in integrals 1 and 11, but the gain does not emerge until high accuracy is requested as in Table IV. For harsher singularities the gain would likely become apparent earlier. References 1. Patterson, T.N.L. The optimum addition of points to quadrature formulae. Math. Comp. 22 (1968), 847-856. 2. Cranley, R., and Patterson, T.N.L. On the automatic numerical evaluation of definite integrals. Comp. J., 14 (1971), 189-198. 3. Lyness, J.N. Algorithm 379, SQUANK. Comm. A C M 13 (Apr. 1970), 260-263. 4. Wallick, G.C. Algorithm 400, Modified Havie integration. Comm. ACM 13 (Oct. 1970), 622-624. 5. de Boor, Carl. CADRE: An algorithm for numerical quadrature. Mathematical Software. J.R. Rice (Ed.) Academic Press, New York, 1971, pp. 417-449. 6. Thacher, H.C. Jr. Remark on Algorithm 60, Comm. A C M (July, 1964), 420-421. Algorithm C C C C C C C C C C C D C C C C C C C C C C C C C C C C C C C SUBROUTINE QUAD(A, O, RESULT, Ks EPSIL, NPTS, ICHECK, F) DIMENSION FUNCTII27), P(381)* RESULT(B) THIS SUBROUTINEATTEMPTS TO CALCULATE THE INTEGRAL OF F(X) OVER THE INTERVAL * A * TO * B * WITH RELATIVE ERROR NDT EXCEEDING * E P S I L * . THE RESULT IS OBTAINED USING A SEQUENCE O F 1 , 3 , 7 , 1 5 s 3 1 , 6 3 , 127, AND 255 PDINT INTERLACING FORMULAE(NG INTEGRAND EVALUATIONS ARE WASTED) OF RESPECTIVE DEGREE 1 , 5 , 1 1 , 2 3 , 47,95,191 AND 3 8 3 . THE FORMULAE ARE BASED ON THE OPTIMAL EXTENSION OF THE 3 - P O I N T GAUSS FD~MULA. DERAILS OF THE FORMULAE ARE GIVEN IN "THE ~PTIMUM ADDITION OF POINT5 TO QUADRATURE FORMULAE' BY T . N . L . PATTERSON,MATHS.COMP. V0L 2 2 , 8 4 7 - B 5 6 , 1 9 6 8 . *** INPUT * * * A L O W E R L I M I T OF INTEGRATION. B UPPER LIMIT OF INTEGRATION. EPSIL R E L A T I V E ACCURACY REQUIRED. WHEN THE RELATIVE DIFFERENCE OF TWO SUCCESSIVE FORMULAE DOES NOT EXCEED *EPSIL* THE LAST FORMULA CDMPUTED IS TAKEN AS THE RESULT, F F(X) I S THE INTEGRAND. *** OUTPUT * * * RESULT T H I S ARRAY, WHICH SHOULD BE DECLARED TO HAVE AT LEAST B ELEMENTS, HOLDS THE RESULTS OBTAINED BY THE 1 , 3 , 7 , E T C - , POINT FORMULAE. THE NUMBER OF FORMULAE COMPUTED DEPENDS ON * E P S I L * . K R E S U L T ( K ) HOLDS THE VALUE OF THE INTEGRAl. TO THE SPECIFIED RELATIVE ACCURACY. NPTS NUMBERINTEGRAND EVALUATIONSICHECK DN E X I T NORMALLY ICHECK=O. HDWEVER I F CBNVERGENCE T O THE ACCURACY REQUESTED I S N o r ACHIEVED I C H E C K = I O N EXIT. ABSCISSAE AND WEIGHTS OF QUADRATURERULES ARE STACKED IN ARRAY * P * I N THE ORDER I N WHICH THEY ARE NEEDED. DATA * P( I ) , P ( 2 ) , P ( 3 ) , P ( 4 ) , P ( 5 ) , P ( 6 ) , P ( 7 ) , * P( B I , P ( 9),P(IO),P(II),P(12)*P(13),P(14), P(15I,P(16I,PIITI,PIIBI,P(19J,P(BO)*P(BI), * P(22),P(B3),P(24),P(2S),P(26I,P(27),P(2B)/ * O.77459666924146337704E 00,0.55555555555555555556E 00, * O.8BBGBBBBBOBOBBOBsSBgE O O 3 0 - 2 6 8 4 G B O B 9 0 6 B 3 3 3 4 4 0 7 3 E O0, * 0.96049126870802028342E 00,0.I0465622602646726519E 00, * 0.43424374934680255BOOE O0, O . 4 D I 3 9 7 4 1 4 7 7 5 9 6 2 2 ~ Q g l E 00, O . 4 5 0 9 1 6 5 3 B 6 5 8 4 7 4 | 4 2 3 5 E O O , O . I S 4 4 1 5 2 5 5 2 4 3 7 8 4 2 2 0 3 6 E DO, * 0.51603eB2997079739697E-OI,O.20062852937698902103E 00, * 0.99303196321275502221E O0, O . 1 7 0 0 1 7 1 9 6 2 9 9 4 0 Q B O 3 3 9 E - 0 1 * Communications of the ACM November 1973 Volume 16 Number 11 * * * * 0.88845923287225699889E 00,0.92927195315124537686E-01~ 0.62ll0294673722640294E OO, O . 1 7 1 5 1 1 9 0 9 1 3 6 3 9 1 3 8 0 7 9 E DO, 0.22338668642896688163E 00,0.21915685840158749640E OO, 0.22551049979820668739E 00,0*67207754295990703540E-01, * O°25807598096176653565E-Ot,O°1003142786|I79557877E 00~ * O.84345657393211062463E-O2~Oo4646289326|757986541E-01~ * O°85755920049990351154E-OI,O.|O9578421D5592463824E OD/ DATA P(29),Pt3O),P(31),P(32),P(33),P(34),P(35), * P(36),P(37)~P(38)*P(39),P(4O),P(4|),P(42), * P(43),P(44),P(45),P(46),P(47),P(48),P(49)* P(50),P(51),P(52)*P(53)*P(S4),P(S5),P(56)/ * 0.9990961~496766759766E 00,0°~54478079|5618744154E-0~, * 0°98153114955374010687E 00,0.16446049854387810934E'01, 0.92965485742974005667E 00,0°359571033071~9322097E-0l, * 0.63672593816886873550E 00,0.56979509494l~33574|2E-0|, 0.70249620649152707861E 00,0.768796~0499003531043E-01, * 0.5313197436443756~397E 00,0o936S7109981264473617E-0|, * 0o33113539325797683309E 0 0 , 0 o 1 0 5 6 6 9 6 9 3 5 8 0 2 3 4 8 0 9 7 4 E DO, * * * * * * * 0.11248894313318662575E 00,0.11195687302095345688E 0.11275525672076869161E O0,O.336D3877146207730542E-Ol, O.12903800100351265626E-OI,0.SOISTI39305899537414E-Ol* O0, O*42176304415588548391E-O2*O.23231446639910269443E-Ol* 0.42877960025007734493E-OI,0.54789210527962865032E-01, O°12651565562300680114E'02,0.8223OO79572359296693E-D2, O.17978551568128270333E-OI,O°284897547458335466|3E-Ol/ DATA P(57),P(58),P(59),P(6O),P(6I),P(62),P(63), P(64),P(65),P(66)*P(67),P(68)*P(69),P(70), * P(7|),P(7~),P(73)~P(74)*P(75)*P(76),P(77), * P(78),P(79),P(SO),P(SI),P(82),P(83),P(84)/ * * * * * O.38439810249455532039E-OI*O.46813554990628OIe403E-O|, 0°528349467901|6519862E-0|,O°559784365104763i9408E-01* 0.99987288812035761|94E 00,0°36322148184553065969E-03, * 0.99720625937222195908E 00,0°25790497946856682724E-02~ * 0.98868475754742947994E 00,0.61155068~21172463397E-02, * 0.972|8287474858179658E 00~0.104982469096213~1898E-01, 0o946342858373402905|5E 00~0.1540675046655949780~E-0l~ * * * * * * * * * * * * * * * * * * * * 0.91037115695700429250E 0o86390793819369047715E 0°80694053195021761186E O0,O.2059423391591eTil149E-OI, 00,0*258696793272|47469|1E-01, 00,0.31073551111687964880E-01~ 0.73975604435269475866E 00,0.36064432780762572640E-01, 0.66290966002476059546E 00,0.407155101169443|8934E'01, 0.5771957|005204581484E 00,0.44914531653632197414E-01, 0°48361802694S84i02756E 00,0.485643304066731987]6E-01/ DATA P( 85),P( 86)~P( 87),P( 88)*P( 89),P( 9O),P( 91)* P( 92),P( 93),P( 94)*P( 95),P( 96)*P( 97),P( 98), P( 9 9 ) , P ( I O O ) , P ( I O I ) , P ( I O 2 ) , P ( I O 3 ) , P ( I O 4 ) , P ( | 0 5 ) , P(IO6),P(IO7),P(ID8),P(IOg),P(|ID)*P(III),P(ll2)/ 0.38335932419873034692E 00,0°5|5~325395~048458777E-01~ 0.27774982202182431507E 00,0*53905499335266063927E-0l, 0.16823525155220746498E 00,0.55481404356559363988E-01, 0.5634431304659278997~E-0|*Oo562776998312543OI273E-0|, O°56377628360384717388E-OI,O.]68Oi938574|O3865271E-01, O.64519000501757369228E-O2,O.25078569652949768707E-O|, Oo21088|524572663~8793E-O~,O.I1615723319955134727E-01, O*21438980012503867~46E-Ol,Oo~739460526398|432516E-Ol, 0.6326073|9369633544~E-O3, O.4|II5039786546930472E-02~ O*~9892757840641357~33E-O2,0oI4244877372916774306E-Ol, O°|9219905%247277660|gE-OI,O°2340677749531400620|E-O|~ O.~6417473395058259931E-OI,O.27989~IS~55238159704E-Ol, O.18073956444538835782E-O3,0.12895240826104173921E-02, * * * * * O.30577534101755311361E-O2,OoS249123454808859|~51E-02/ * DATA * * * P(II3),P(II4),P(IIS),P(116),P(117),P(118)*P(119), P(120),P(121),P(122),P(123),P(124),P(125),P(t26), * * P(127),P(128),P(129),P(|30),P(131),P(132),P(|33), O*98351865757863272876E 00*Oolg197129710135724125E-02, 0o97940628167086268381E 00~0°21944069253638388388E-02, 0.97473445975240266776E 00,0.24789582266575679307E'02/ DATA P(281),P(282),Pt283),P(284),Pt285),Pt286)*F(287)~ p(288),P(289),P(290),P(291)*P(292),P(293),P(294), p(295),P(296),P(297)~P(298)*P(299),Pt300),P(301)* P(302),P(303),P(304),P(305),P(306),P(307),P(308)/ 0°96948465950245923177E 00,0.27721957645934509940E-02, 0.96364062156981213252E 00,0.30730184347025783234E-02, 0.9 5718821610986096274E 00,0*33803979910869203823E-02, * P¢I34),P(135),P(136)~P(137),P(138),P(139)~P(140)/ * Oo7703375233279741~482E-O2~OolO2971|6957956355524E-01~ Oo|~934839663607373455E-O|,Oo|5536775555843982440E-O~, * * * * * * * * * * * * 0.95011529752129487656E * * 0.9424|156519108305981E 0.93406843615772578800E * 0.92507893290707565236E * * * 0.76611781930376009072E 0.7486962936169366028~E 0.73066452124218126133E 00,0.85565435613076896192E-02, * 0.71203315536225203459E 00,0.94636899938300652943E-0~, * * * * * * 0.692813769779||470289E 0.180322|639039|2863~OE-O|,0*20357755D58472159467E-0|, 0o224572658268|6098707E-01,O.24282165203336597358E-01, 0.25791626976D24229388E-Ol*0o2695274966763303|963E-01, 0o27740702178279681994E-01,0.28138849915627150636£-01, 0.99998243035489|59858E 00,0.50536095207862517625E-04, 0.99959879967|9|068325E 00,0o37774664632698466027E-03, O.99B31663531840739253E 00,0.93836984854238150079E-03~ 0.995724|04698407|8851E 00,0o|68||428654214699063E-0~, 0 o 9 7 | 4 9 5 7 ~ | | 7 8 1 0 6 1 3 ~ 4 0 E DO, O o 2 5 6 8 7 6 4 9 4 3 7 9 4 0 ~ O 3 7 3 | E - O ~ , * 0o98537|49959852037|1|E 00,0.357~89~7835172996474E-0~, * 0o977|4|5|46397057|416E 00,0.46710503721|43~|7474E-0~ 0.96663785|5584|656709E 00,0o584344987583563~5076E-02/ DATA * P(|41),P(142),P(|43)*P(144),P(145)*P(|46),P(147), * P(148),P(149),P(150),P(151),P(152),P(153),P(154), * p(|55)~P(156)*P(157)~P(158),P(159),P(|60)~P(|61), * p(162)*P(163),P(164),P(|65),P(|66),P(167)*P(168)/ + 0.953730006425761|364|E 00,0o70724899954335554680E-09, * 0o93832039777959~88365E 00,0o83428387539681577056E-02, 0o9203400~54700|24~073E 00~0.9641|7772970~5366753E-02, 0.89974489977694003664E 00,0.|095573338783790|648E-0|, * 0o8765|34|4484705~6974E 00,0.12275830560082770087E-01, * 0o85064449476835027976E 00,0.1359|57|009765546790E-01, 0o822156254364980~0737E 00,0.|489~6416648|5|82035E-0|, * 0o79106493379984836|43E 00,0o|6|732|87295777|994SE-0|, 0.7574839663805|363793E 00,0.|742|930|59464|73747E-0|, * * * * * * * * * * * * * * * 697 O.84454040083710883710E-O|*O.28076455793817246607E-Ot= 0.28184648949745694339E'O|pO.281763|9033016602131E'OIJ 0.28|88814180192358694E-Ol,O.84009672870519326354E-O2P O*32259500250878684614E-O2*O.]2539284826474884353E-O1, * O.I0544076228633167722E-O2,0.58078616599775673635E-D2, O.I0719490006251933623E-O|,OoI3697302631990?I6258E-O|/ DATA P(197)*P(t98),P(t99),P¢2OO),P(201)*P(202)*P(203)* * P(204),P(205),P(206),P(20?),PC208)JP(209),P(210)~ * P(211),P(212),P(213),P(2[4),P(215)sP(216),P(217), * Pt218),P(219),P(22O),P(221),P(222),P(223)*P(224)/ * O.31630366082226447669E-O3,0*20557519893273465236E-02* : Oo44946378920320678616E-O2,0.71224386864583871532E-02, 0.96099525623638830097E-O2~O.IITO338874765?OO31OlE-Ol* O.13208736697529129966E-OlsO.13994609127619079852E-01, * O.90372734658751149261E-O4nO.64476204|30572477933E-03, OoI52887670508T7655684E-O2*O.262456|7274D44295626E-O2a O*38516876166398709241E-O2,O.51485584789761777618E-02* O.64674198318036867274E-O2,0.776838777792199122OOE'02, O.90161081951756431600E-O2,0.lOi78877529236079733E-Ol* * 0.l]228632913408049354E'01,0.12141082601668299679['01, * Oo12895813488012114694E-OI,OoI34763748338165|5982E-01* Oo13870351089139840997E-O|*O.14069424957813575318E-Ol, O.25157870384280661489E'O4*O*18887326450650491366E-O3s O.46918492424785040975E-O3~O.84057143271O72246365E-03/ DATA * P(225)aPt226),P(227),Pt228),P(229)pP(230)*P(23t)s * P(232),P(233),P(234),P(235),P(236),P(237),P(238), * p(239),P(24O),P(24|),P(242)*P(243),P(244)*P(245), P(246),Pt247),P(248)*P(249)JP(25O),P(251)mP(252)/ * O.1284382471~9701OI766E-O2JO.I7864463917586498247E-D2s * O.23355251860571608737E-O2~Oo29217249379178197538E'02* 0.35362449977167777340E-O2,0.417t4193769840788528E-02, * 0.48205888648512683476E-O2,0.5477866693918950824OE-02, Oo61379|5280D413850435E-O2*Oo67757855048827733948E-02~ * O°74468208324075910174E-O2,0*8086609364788859971OE-02, * 0.87109650797320868736E-D2, O.9315924|280693950932E-02* * O°98977475240487497440E-O2*O°ID452925722906OII926E-01, * OolO97818315265891247OE-OlPO°l147048211469387438OE-01* * O.l1927026053019270040E-OlsO.12345262372243838455E-Ol, * O*12722884982732382906E-Ol,O.1305783668835304884OE-Ol* * O.13348311463725179953E-OI,O.I359275661481239591OE-01, * O.13789874783240936517E-Ol*O°13938625738306850804E-OI, O.14038227896908623303E-Ol,Oo14088159516508301065E-Ol/ DATA * p(253)*P(254)~P(255)sP(256)*P(257),P(258),P(259), * P(26O),P(261P, PC262)*P(263),P(264),P(265)~P(266), * P(267),PC268)eP(269)*P(27O),P(271),Pg272),P(273)* * P(274),P(275),P(g76),P(277),P(278),P(279),P(280)/ * 0.99999759637974846462E 00,0o69379364324108267170E-05, 0.99994399620705437576E 00,0.53275293669780613125E-04, 0.99976049092443204733E 00,0o|3575491094922871973E-03~ 0.99938033802502358193E 00,0.24921240048299729402E-03, * 0.99874561446809511470E 00,0.38974528447328229322E-03, O.9978D535449595727456E 00,0.5542953149303747|492E-03, 0.9965141459|489027385E 00,0.74028280424450333046E-03, O.994E31502800621OOO52E O0*0.94536151685852538246E-D3, * 0.99272134428278661533E 00,0°11674841174299594077E-02, Oo990|5137D400770|5918E 00~0.14049079956551446427E-02, O*987092527954D34067|9E OO*O.1656112728|544526D52E-02, * 0.72t4230853700989|548E 0.68298743109107922809E 00,0.|8631848256|38790186E-0|, 00,0o19795495048097499488E-0|, 0o6422766425097595|377E 00~0.209058514458|2023852E-0|, 0.59940393024224289297E 00,0.21956366305317824939E-01, 0.55449513263193254887E 00,0.2294096422938774876|E-01/ DATA P(169),P(170)*P(17|),P(|72),P(|73),P(|74),P(175), P(176),P(177),P(178),P(|79),P(|80),P(|81),P(182), P(183),p(|84)*P(185),P(186),P(|67),P(188)~P(189)~ P(|9D)*P(|9|),P(|9~)*P(|93),P(194),P(]95),P(|96)/ 0o5076877575337|660215E 00,0o23854052|06038540060E-0|, 0.45913001198983233287E 00,0o24690524744487676909E-01, 0.40897982122988867241E 00,0.254457699654647658|3E-01, OO, O.~6|1567337670609768OE-01, 00,0.26696622927450359906E-0|, 00,0o27185513229624791819E-01, 0o35740383783153215236E 0o3045764415567|404334E 0.25067873030348317661E 0o1958975027|110015392E O0*Oo2757974956646|873035E-Ot, 0.14042423315256017459E 00,0.27877251476613701609E-0|, 00~0.36933779170256508183E-02~ 00,0.40110687240750233989E-02, 00,0o43326409680929828545E-02, 00,0°46573172997568547773E-02, 0o9|543758715576504064E 00,0.49843645647655386012E-02, * 0°9051403568|326|59519E 00,0.53|3086605|870565663E-02, * 0.89418456833555902286E 00,0.56428|810|3844441585E-02* * 0.88256884024734190684E 00,0o59729195655081658049E-02~ 0o870293055548||390585E 00,0.63027734490857587172E-02, * 0.85735831088623215653E 00,0.663|7812429018878941E-02, * 0*84376688267~70860104E 00,0*69593614093904229394E-02/ DATA * P(309),P(31O),P(31t),P(312),Pt313),P(314)*P(315), * P(316),P(317),P(318),P(319),P(32D),P(321),P(322), * p(323),P(324),P(325),P(326)~P(327)~P(328),P(329)~ P(330),P(33i)*P(332),P(333),P(334),P(335),P(336)/ * 0.82952219463740140018E 00,0.72849479805538070639E-02, * O.814628787655|374t344E 00,0.76079896657|90565832E-02, * 0°79909229096084|40|80E 00,0.7927949334294849||03E-02, D.78291939411828301639E DO,Oo82443037630328680306E-02, * * * * 0.67301883023041847920E 0*6526616654100174961DE 0o63175643771119423041E Oo6103181t37151864DO16E 0.58836243444766254143E DATA 00,0.88641732094824942641E-02, 00,O°91667111635607884067E-02, 00,0.97546565363174114611E-02, 00,0.10039172044056840798E-01, 00,0.|0316812330947621682E-01, OO,OolO587167904885|97931E-Dl* 00~0.10849844089337314099E-01~ O0,Ootl104461134006926537E-O|/ P(337),P(338),P(339)~P(340),P(341),P(34~)~P(343), P(344),P(345),P(346),P(347)*P(348),P(349)~P(350), P(351),Pt352),P(353),P(354),P(355),P(356),P(357), P(358),P(359),P(36O),P(36I),P(362),P(363),P(364)/ 0o56590588542365442262E 00,0°113506543|5980596602E-01, * * * 0.54296566649831149049E 0.51955966153745702199E D.49570640791876146017E * 0.471425065~7165887693E 00~0.12~44424981611985899E-01, 0.4467353876620~847374E 00,0.12443560190714035263E-01, * * * * 00,0.11588074033043952566E-01, O0,0.tISt6385890830235763E-0t, 0D,Oo1203527D785279562630E-01, 0o42165768662616330006E Do39621280605761593918E 00,0°12632403643542078765E-01~ 00,0.128|0698163877361967E-01, 0.37042208795007823014E 00~0.1~978202239537399286E-01, Communications of the ACM November 1973 Volume 16 Number 1l * O. 34430734159943802278E DO, O . 1 3 1 3 4 6 9 0 0 9 1 9 6 0 1 5 2 8 3 6 E - O 1, O. 3 I T8908 IBO6847668318E DO* O . 1 3 2 7 9 9 S 1 7 4 3 9 3 0 5 3 0 6 5 0 E - 0 1 . 0.29119514851824668196E DO, 0 . 1 3 4 l 3 7 9 3 0 8 5 1 1 0 0 9 8 5 1 3 E - O k, e O . B 6 4 B A S S T B A I O g B 6 7 6 1 9 4 E DO* 0 . 1 3 5 3 6 0 3 5 9 3 4 9 5 6 B 13614E-O I * * O.B3705884558982972721E GO, O . 1 3 6 4 6 5 1 8 1 0 2 5 7 1 2 9 1 4 2 8 E ' 0 1 / DATA * P(365)*P¢366),P(367),P(368)*P(369),P(STO),P(371 )* * P(372), P(373), P(374)* P(375), P(Q76), P(3?T), P(3TO), .1= P ( S 7 9 ) , P ( 3 8 G ) , P ( S 8 1 )/ 0.209665238243181194TTE OO, O . 1 3 7 4 5 0 9 3 4 4 3 0 0 1 8 9 6 6 3 2 E - 0 1 * * O. 182086496759252 | 9BB5E OO~O.13831631909SO6428676E-Ol, * O. 1 5 4 3 4 6 8 1 1 4 8 1 3 7 8 I O O 6 9 E GO, O . 1 3 9 0 6 0 1 9 6 0 1 3 2 5 4 6 1 2 6 4 E - O I, * O • 12647058437230196685E OO, O. 1 3 9 6 8 1 5 8 8 0 6 5 1 6 9 3 8 5 1 6 E - O I* * o.g8482396598119202090E-OI*O.14OI?96GO39456608GIOE°OI, * 0.'/04069760428 SS179063E-01 * O • 1405538BO7264996427TE-01* * 0 • 42269164765363603212E-O 1, O . 1 4 0 G O S 5 1 9 6 2 5 5 3 6 6 1 3 2 5 E - 0 1 , * O*I4093886410782462614E-OI~O.IAO92OASO691604OBSSSE-01* * O • I 409A40"/0900961T934TE-O I/ ICHECK = 0 CHECK FOR TRIVIAL C A S E . IF (A.EQ.B) GO T O TO. SCALE FACTORS. SUM = ¢ B + A ) / 2 . O DIFF= (B-A)/B.O I-POINT GAUSS FZER0 = F(SUM) RESULT( l ) = B.O*FZEROeDI FF I = O IOLD = 0 INEW = 1 R = 2 ACUM == O . O GO TO 30 I0 IF (K.EG.G) GO TO SO K = K ÷ I ACUM = O.O C O N T R I B U T I O N FROM F U N C T I O N V A L U E S A L R E A D Y COMPUTED. DO BO J = i p I O L D I = I + I ACUM = ACUM + P ( 1 ) * F U N C T ( J ) 20 CONTINUE CONTRIBUTION FROM NEW F U N C T I O N V A L U E S . 3 0 IOLD = IOLD ÷ INEW D0 40 J = I N E W , 1 0 L D I = I + I X = P(1)*DIFF FUNCT(J) = F(SUM+X) + F(SUM-X) I = I + I ACUM = ACUM + P ( I ) * F U N C T ( J ) 40 CONTINUE INEW = IBLD + I I ffi I + I RESULT(K) = (ACUM+P(I)*FZERO)~tDIFF CHECK FOR C O N V E R G E N C E . IF (ABS(RESULT(K)-RESULTCK-I))-EPSIL*ABS(RESULT(K))) 60* * 60. 1O CONVERGENCE NOT A C H I E V E D . 50 I C H E C K = 1 NORMAL TERMINATION. 60 NPTS = I N E W + I O L D RETURN TRIVIAL CASE 70 K = 2 RESULT(I) = O.O RESULT(B) = O.O NPTS = O RETURN END * * C C C C C C C C C F U N C T I O N G S U B ( A * B* E P S I L s N P T S * I C H E C K * RELERR* F) C THIS FUNCTION ROUTINE PERFORMS AUTOMATIC INTEGRATION C OVER A F I N I T E INTERVAL USING THE BASIC INTEGRATION C ALGORITHM QUAD* TOGETHER WITHJ I F NECESSARY* A NBNC A D A P T I V E SUBDIVISION PROCESS. C THE CALL TAKES THE FORM C G S U B ( A * B* E P S I L t N P T S , I CHECK~ RELERR~ F ) C AND CAUSES FOX) TO BE INTEGRATED OVER ¢A*B) WITH R E L A T I V E C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C ERROR H O P E F U L L Y NOT E X C E E D I N G E P S I L . SHOULD QUAD CONVERGE (ICHECK=O) THEN QSUB W I L L RETURN T H E V A L U E O B T A I N E D BY I T OTHERWISE SUBDIVISION WILL BE INVOKED AS A RESCUE OPERATION I N A NON-ADAPTIVE MANNER. THE ARGUMENT RELERR G I V E S A CRUDE E S T I M A T E OF T H E A C T U A L R E L A T I V E ERROR OBTAINED. THE SUBDIVISION STRATEGY I S AS FOLLOWS LET THE INTERVAL CA, B) BE DIVIDED INTO B * * N PANELS AT STEP PROCESS. OUAD IS A P P L I E D F I R S T T O N O F THE S U B D I V I S I O N THE SUBDIVIDED INTERVAL ON WHICH QUAD LAST FAILED TO C O N V E R G E A N D IF C O N V E R G E N C E IS N O W A C H I E V E D T H E R E M A I N I N G P A N E L S ARE I N T E G R A T E D . SHOULD A CONVERGENCE F A I L U R E OCCUR O N ANY PANEL T H E INTEGRATION AT T H A T POINT IS T E R M I N A T E D AND THE PROCEDURE R E P E A T E D W I T H N I N C R E A S E D BY 1o THE STRATEGY INSURES THAT POSSIBLY DELINQUENT INTERVALS ARE E X A M I N E D BEFORE WORK, WHICH L A T E R M I G H T H A V E TO BE D I S C A R D E D * I 5 I N V E S T E D ON WELL B E H A V E D P A N E L S . T H E PROCESS I S COMPLETE WHEN NO CONVERGENCE F A I L U R E OCCURS ON ANY PANEL AND THE SUM OF THE RESULTS OBTAINED BY QUAD ON EACH P A N E L I S T A K E N AS THE V A L U E OF THE I N T E G R A L . THE PROCESS I S VERY CAUTIOUS I N T H A T THE SUBDIVISION OF T H E I N T E R V A L (A~B) IS UNIFORM, T H E F I N E N E S S O F WHICH IS CONTROLLED BY T H E SUCCESS 0 F QUAD. IN THIS WAY I T I 5 RATHER DIFFICULT F0R A SPURIOUS CONVERGENCE TO SLIP THROUGH, T H E CONVERGENCE C R I T E R I O N OF QUAD I S S L I G H T L Y R E L A X E D I N T H A T A P A N E L I S DEEMED TO H A V E BEEN S U C C E S S F U L L Y INTEGRATED I F EITHER QUAD CONVERGES OR THE ESTIMATED A B S O L U T E E~RROR C O M M I T T E D ON T H I S P A N E L DOES NOT EXCEED E P S I L T I M E S THE E S T I M A T E D A B S O L U T E V A L U E OF T H E I N T E G R A L OVER ( A . B ) . THIS RELAXATION IS TO TRY TO TAKE ACCOUNT 0 F A COMMON SITUATION WHERE ONE PARTICULAR PANEL CAUSES SPECIAL DIFFICULTY, PERHAPS DUE T O A S I N G U L A R I T Y O F S O M E TYPE. I N T H I S CASE QUAD COULD O B T A I N NEARLY EXACT ANSWERS ON ALL OTHER PANELS AND 50 THE RELATIVE ERROR FOR THE TOTAL INTEGRATION WOULD BE ALMOST ENTIRELy DUE T0 THE DELINQUENT PANEL° WITHOUT THIS CONDITION THE COMPUTATION MIGHT CONTINUE DESPITE THE REQUESTED RELATIVE ERROR BEING ACHIEVED. 698 C C C C C C C C C C C C C C C C C C C C C C C C C C THE I N T E G R A T I O N I S I N D I C A T E D BY I C H E C K , CONVERGENCE O B T A I N E D W I T H G U T 1 N V O H I N G SUBDIVISION. T H I S CORRESPONDS TO THE D I R E C T USE OF Q U A D . ICHECK=I RESULT OBTAINED AFTER INVOKING SUBDIVISION. ICHECK=B AS FOR I C H E C K = 1 B U T A T SOME P O I N T THE R E L A X E D CONVERGENCE C R I T E R I O N WAS U S E D . THE R I S K O F U N D E R E S T I M A T I N G THE RELATIVE ERROR W I L L BE I N C R E A S E D . IF NECESSARY, C O N F I D E N C E MAY BE RESTORED BY C H E C K I N G E P S I L AND RELERR FOR A S E R I O U S D I S C R E P A N C Y . ICHECK NEGATIVE I F D U R I N G THE S U B D I V I S I O N PROCESS T H E ALLOWED UPPER L I M I T ON T H E NUMBER OF P A N E L S THAT MAY B E G E N E R A T E D (PRESENTLY 4096) IS REACHED A R E S U L T I S OBTAINED WHICH MAY B E U N R E L I A B L E BY C O N T I N U I N G T H E I N T E G R A T I O N W I T H O U T FURTHER S U B D I V I S I O N IGNORING CONVERGENCE F A I L U R E S . T H I S OCCURRENCE I S FLAGGED BY R E T U R N I N G I C H E C K WITH N E G A T I V E SIGN. THE R E L I A B I L I T Y OF T H E ALGORITHM WILL DECREASE FOR L A R G E V A L U E S OF E P S I L . I T IS RECOMMENDED T H A T E P S I L SHOLa-D G E N E R A L L Y BE LESS T H A N ABOUT O . O O l . DIMENSION RESULT(O) I N T E G E R B A D , OUT L O G I C A L RHS EXTERNAL F DATA N M A X / 4 0 9 6 / C A L L G U A D ( A , B , R E S U L T s K* E P S I L , NPTS, ICHECK, F) QSUB = R E S U L T ( K ) RELERR = O.O IF (QSUB.NE.O.O) R E L E R R = * ABS((RESULT(K)-RESULT(K-I))/QSUB) CHECK I F S U B D I V I S I O N IS NEEDED. IF IICHECK.EG.O) RETURN SUBDIVIDE ESTIM = ABS(QSUB*EPSIL) IC = I RHS = . F A L S E . THE OUTCOME 0 R ICHECK=O - N = I H = B - A BAD = 1 GSUB = O . O RELERR = O . O H = H*O.5 N=N+N INTERVAL ( A , B ) DIVIDED INTO N EQUhJ- SUBINTERVALS. INTEGRATE OVER SUBINTERVALS BAD TO (BAD+I) WHERE TROUBLE HAS OCCURRED. H I = BAD MB = BAD ÷ 1 OUT = 1 GO T O SO I N T E G R A T E OVER S U B I N T E R V A L S I TO ( B A D - I ) 2 0 MI = I M2 = BAD - 1 RHS = . F A L S E . OUT= 2 GO TO 5 0 INTEGRATE OVER SUBINTERVALS IBAD+B) TO N. 30 MI = BAD + 2 M2 = N OUT = 3 GO TO 5 0 SUBDIVISION R E S U L T 40 I C H E C K = I C RELERR = RELERR/ABS(OSUG) RETURN INTEGRATE OVER SUBINTERVALS MI TO M2. 50 I F I M I . G T . M 2 ) GO T 0 9 0 DO BO J J = M 1 s M B J = JJ E X A M I N E F I R S T THE L E F T OR R I G H T H A L F OF THE S U B D I V I D E D TROUBLESOME I N T E R V A L D E P E N D I N G ON T H E OBSERVED T R E N D . I F ( R H S ) J = MO + H I - J J ALPHA = A + H*¢J-I) BETA = ALPHA ÷ H CALL QUAD(ALPHAs B E T A . R E S U L T s Ms E P S I L J NF, I C H E C K , COMP = ABS(RESULT(M)-RESULT(M-I)) NPTS = NPTS + NF IF (ICHECK.NE.I) GO T0 7 0 I F (COMP.LE.ESTIM) G 0 T O lOG S U B I N T E R V A L J H A S CAUSED T R O U B L E , CHECK I F FURTHER S U B D I V I S I O N SHOULD BE C A R R I E D O U T . I F (N.EQ.NMAX) GO TO 60 BAD = B * d l RHS = . F A L S E . IF ((J-2*(J/2)).EQ.O) RHS = . T R U E G O TO 10 60 IC = -IABSCIC) TO QSUB = GSUB + RESULT(M) BO C O N T I N U E RELERB = RELERR + COMP 9 0 GO T O ( 2 0 , 3 0 , 4 0 ) , OUT R E L A X E D CONVERGENCE lOG I C = I S I G N ( 2 , 1 C ) GO TO TO END I0 C C C C C C C C C C C - C FUNCTION DSUBA(AJ B, EPSILs NPTS* ICHECK* RELERR* F) C THIS FUNCTION ROUTINE PERFORMS AUTOMATIC INTEGRATION C OVER A F I N I T E INTERVAL USING THE BASIC INTEGRATION C ALGORITHM GUAD TOGETHER WITH, ~F NECESSARY AN ADAPTIVE C SUBDIVISION PROCESS. I T I S GENERALLY MORE E F F I C I E N T THAN C C C C C C C C C C THE N O N - A D A P T I V E A L G O R I T H M QSUB BUT I S L I K E L Y TO BE L E S S COMP.J.*I4*IB9*1971). THE CPJ-L TAKES THE FORM QSUBA(A*B*EPSIL*NPTS*ICHECK,HELERRsF) AND CAUSES F ( X ) TO BE INTEGRATED OVER CA, B) WITH RELATIVE ERROR H O P E F U L L Y NOT E X C E E D I N G E P S I L . SHOULD QUAD CONVERGE (ICHECK=O) THEN QSUBA WILL RETURN THE VALUE OBTAINED BY I T OTHERWISE SUBDIVISION WILL BE INVOKED AS A RESCUE OPERATION I N AN ADAPTIVE MANNER. THE ARGUMENT RELERR GIVES A CRUDE E S T I M A T E OF T H E A C T U A L R E L A T I V E ERROR O B T A I N E D . RELIABLE(SEE Communications of the A C M November 1973 Volume 16 Number 11 F) C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C THE S U B D I V I S I O N STRATEGY I S AS FOLLOWS AT EACH STAGE OF THE PROCESS AN INTERVAL I S PRESENTED FOR SUBDIVISION ( I N I T I A L L Y THIS WILL BE THE WHOLE INTERVAL (A, B I ) . THE INTERVAL IS HALVED AND QUAD APPLIED TO EACH SUBINTERVAL. SHOULDBUAD FAIL ON THE FIRST SUBINTERVAL THE SUBINTERVAL I S STACKED FOR FUTURE S U B D I V I S I O N AND THE SECOND SUBINTERVAL IMMEDIATELY EXAMINED° SHOULD QUAD F A I L ON THE SECOND SUBINTERVAL THE SUBINTERVAL I S IMMEDIATELY S U B D I V I D E D AND THE WHOLE PROCESS REPEATED. EACH TIME A CONVERGED RESULT I S OBTAINED I T I S ACCUMULATED AS THE P A R T I A L VALUE OF THE I N T E G R A L . WHEN OUAD CONVERGES ON BOTH SUBINTERVALS THE INTERVAL LAST STACKED I S CHOSEN NEXT FOR S U B D I V I S I O N AND THE PROCESS REPEATED. A SUBINTERVAL I S NOT EXAMINED A G A I N ONCE A CONVERGED RESULT I S OBTAINED FOR I T 50 THAT A SPURIOUS CONVERGENCE I S MORE L I K E L Y TO S L I P THROUGH THAN FOR THE N O N - A D A P T I V E ALGORITHM QSUB. THE CONVERGENCE C R I T E R I O N OF QUAD I S S L I G H T L Y RELAXED I N THAT A pANEL I S DEEMED TO HAVE BEEN SUCCESSFULLY INTEGRATED I F EITHER DUAD CONVERGES OR THE ESTIMATED ABSOLUTE ERROR COMMITTED ON T H I S PANEL DOES NOT EXCEED E P S I L TIMES THE ESTIMATED ABSOLUTE VALUE OF THE INTEGRAL OVER ( A ~ B ) . THIS RELAXATION IS TB TRY TO TAKE ACCOUNT OF A COMMON S I T U A T I O N WHERE ONE PARTICULAR PANEL CAUSES SPECIAL D I F F I C U L T Y , PERHAPS DUE TO A S I N G U L A R I T Y OF SOME TYPE. IN THIS CASE QUAD COULD OBTAIN NEARLY EXACT A N S W E R S ON A L L O T H E R P A N E L S A N D SO T H E R E L A T I V E ERROR FOR THE TOTAL INTEGRATION WOULD BE ALMOST ENTIRELY D U E TO T H E DELINQUENTPANEL. WITHOUT THIS CONDITION THE COMPUTATION MIGHT CONTINUE DESPITE THE REQUESTED R E L A T I V E ERROR BEING ACHIEVED. THE OUTCOME OF THE INTEGRATION IS I N D I C A T E D BY I C H E C K . C ICHECK=O CONVERGENCEOBTAINED WITHOUT INVOKING SUBC DIVISION. THIS WOULD C O R R E S P O N D TO T H E C DIRECT USE OF QUAD* C ICHECK=I RESULT OBTAINED AFTER INVOKING SUBDIVISION, C ICHECK=2 AS FOR ICHECKBI BUT AT SOME P O I N T THE C RELAXED CONVERGENCE CRITERION WAS USED. C THE RISK OF UNDERESTIMATING THE R E L A T I V E C E R R O R W I L L B E INCREASED. IF N E C E S S A R Y s C CONFIDENCE MAY BE RESTORED BY CHECKING C E P $ I L AND RELERR FOR. A SERIOUS DISCREPANCY. C ICHECK NEGATIVE C I F DURING THE S U B D I V I S I O N PROCESS THE STACK C OF DELINQUENT INTERVALS BECOMES FULL ( I T I S C PRESENTLY SET TO HOLD AT MOST lOB NUMBERS) C A RESULT I S OBTAINED BY CONTINUING THE C INTEGRATION IGNORING CONVERGENCE FAILURES C WHICH CANNOTBE ACCOMMODATED ON T H E STACK. C T H I S OCCURRENCE I S FLAGGED BY RETURNING C ICHECK WITH NEGATIVE SIGN, C THE RELIABILITY OF THE ALGORITHM WILL DECREASE FOR LARGE C VALUES OF E P S I L . I T I S RECOMMENDED THAT EPSIL SHOULD C GENERALLY BE LESS THAN ABOUT O , O O l , DIMENSION RESULT(O)~ STACK(IOG) EXTERNAL F DATA ISMAX/IOO/ CALL OUAD(A, BJ RESULTs K~ EPSILJ NPTSs ICHECK* E) QSUBA = R E S U L T ( K ) RELERR = O.O I F (DSUBA*NE.O°O) * RELERR = ABS((RESULT(K)-RESULT(K-I))/QSUBA) C CHECK IF SUBDIVISION I S NEEDED IF (ICHECK.EQ.O) RETURN C SUBDIVIDE ESTIM = A B S ( Q S U B A * E P S I L ) RELERR = O.O GSUBA = O.O IS= 1 IC = 1 SUBI = A SUB3 = B 10 SUBB = ( S U B I + S U B 3 ) * O , 5 CALL QUAD(SUBIs SUBBt RESULTJ K , E P S I L ~ NF, ICHECK~ F) NPTS = N P T S ÷ NF COMP = A B S f R E S U L T ( H ) - R E S U L T t K - I ) ) IF (ICHECK,EQ.O) GO TO 30 I F I C O M P * L E . E S T I M ) GO TO 7 0 IF (IS.GE*ISMAX) GO TO 20 C STACK SUBINTERVAL I S U B I ~ S U B B ) FOR FUTURE EXAMINATION S T A C K ( I S ) = SUBI IS = IS + l STACK(IS) = SUBB IS = IS + I GO TO 40 20 IC = - I A B S ( I C ) 30 QSUBA = QSUBA + R E S U L T ( K ) RELERR = RELERR + COMP 40 CALL QUADtSUBBn SUB3s R E S U L T ~ Ks E P S I L s NF~ ICHECK~ F) NPTS = NPTS + NF COMP = A B S ( R E S U L T ( K I - R E S U L T ( K - I ) I I F (ICHECK.ED.O) GO TO 50 IF (COMP.LE,ESTIM) G0 TO 80 C SUBDIVIDE INTERVAL ( S U B B , S U B 3 ) SUB1 = SUB2 GO TO 1O SO QSUBA = DSUBA + R E S U L T ( K ) RELKRR = RELERR + CBMP IF (IS,EQ,I) GO TO 60 C SUBDIVIDE THE DELINGUENT INTERVAL LAST STACKED IS = IS - I SUB3 = STACK(IS) IS = IS - I SUB1 = S T A C K ( I S ) GO TO t O C SUBDIVISION RESULT 6 0 ICHECK = ID R E L E R R = RELERR/ABS(BSUBA) RETURN C RELAXED C O N V E R G E N C E 70 I C = ISIGN(2,1C) GO TO 30 80 IC = I S I G N ( 2 ~ I C ) GO TO 50 END 699 Algorithm 469 Arithmetic Over a Finite Field [A1] C . L a m * a n d J. M c K a y t [ R e c d . 8 O c t . 1971] * Department of Mathematics, Caltech University, Pasadena, CA 91101 t School of Computer University, P.O. Box 6070, Montreal Science, McGill 101, P . Q . C a n a d a Key Words and Phrases: algebra; CR Categories: 5.19 Language: Algol Description The rational operations of arithmetic over the finite field Fq, of q = p~(n >_ 1) elements, may be performed with this algorithm. On entry a[i] contains al E F~ with 0 _< ai < P, i = 0,. • •, n -- 1, and x E Fq satisfies the primitive irreducible polynomial n-1 k P(x) = x ~ -t- Y~k=oakx . f q p r o d u c e s e l i n e [ i ] , i = - - I , . . . , q -- 2, where 1 @ x i = x ei with the convention that -- 1 represents • and x* = O. During execution the range o f the a~ is altered to - p < ai ~ 0, i = 0 . . . n -- 1. The storage used is 2q + n -F 6locations including the final array e. With appropriate conventions for *, multiplication and division are trivial, and addition and subtraction are given by x ~ + x b = xa(1 + x b-a) f o r a < b a n d x a - x b = x a + xt(~-I) x b w h e n p # 2. F o r small values of q, it is suggested that addition and multiplication tables be generated by this algorithm. A description of the method and its generalization to a multi-step process when n is composite is in [2]. A list o f primitive irreducible polynomials is given in [1 ]. Further useful information (especially for p = 2) is to be found in [3]. References 1. Alanen, A.J., and Knuth, D.E. Tables of finite fields. Sankhy6- (A) 26 (1964), 305-328. 2. Cannon, J.J. Ph.D. Th., 1967 U. of Sydney, Sydney, Australia. 3. C o n w a y , J.H., and Guy, M.J.T. Information on finite fields. In Computers in Mathematical Research. North-Holland Pub. Co., Amsterdam, 1967. Algorithm procedure fq(p, n, a, e); integer p, n; integer array a, e begin integer array c[0 :n-- 1 ], f[O :p T n - 1 ]; integer q, i, j, d, s, w; q:=P T n; for i : = 0 step 1 until n - 1 do ira[i] # 0 then a[i] : = a[i] -- p; for i : = 1 step 1 u n t i l n -- 1 doe[i] : = O; c[O] : = 1; f [ l ] : = 0; f[0] : = --1: fori:= lstepluntilq-- 2do begin d : = e [ n - - 1]; s : = O ; f o r j : = n -- 1 step --1 until 1 do begin w:=c[i--1]--dXa[j]; w:= w-- w+pXp; c [ j ] : = w; s : = p X s + w end; w:= --dXa[0]; w : = w - - w + p X p ; c [ 0 ] : = w; f [ p X s + w] := i end; for i : = q step - - p until p do begin e[f[i-- 1]] : = f[i--p]; f o r j : = i -- p step 1 until i -- 2 do e[f[j]] : = f [ j + l ] end end Communications of the A C M November 1973 Volume 16 N u m b e r 11