- Email: [email protected]
The use of recursive procedures in ALGOL 60
The Use of Recursive Procedures in ALGOL 60 H. RUTISHAUSER Swiss Federal Institute of Technology, Zurich
INCREASING emphasis has since recently been given to recursive procedures, i.e, procedures that call themselves. The present comments may be considered as a contribution to the debate that arose over this subject.
I. INTRODUCTION
Unfortunately recursive procedures have been advertised by examples such as real procedurefact(n) ; value n ; integer n ; fact: = if n=O then I else n Xfact(n-l) ;
or real procedure hornerl(a,k,z,n) ; value n,z,k ; real z ; integer k.n ; array a ; hornerl: = if k=n then a[n] else(z X hornerl(a,k+l,z,n)+a[kJ)
True, these examples look quite elegant and may serve to exemplify recursive procedures, butifa computing process canjust as well be described by a simple induction loop, then transforming it into a recursive procedure is sheer extravagance and might well discredit the concept of recursivity altogether. Of course, the solution of many problems present themselves directly in recursive form, but this should not be used for the numerical solution unless it can hardly be avoided. The present author sees mainly two classes of useful applications of recursive procedures. II. INDIRECT RECURSIVITY
Let Xbe an already written non-recursive procedure, one parameter of which represents itself a procedure, e.g, 4-3
H. Rutishauser procedure X(Y) ; procedure Y ;
begin
Y( ... ) ;
_ _ call of procedure Y
end X Now the user of procedure X has to provide the actual parameter corresponding to the formal parameter Y, i.e. he has to declare in his program a procedure Z as the actual counterpart ofY. The body of Z may again call for execution of X as follows:
Program: begin
}--
.
other declarations for the program
procedure Z ( ... ) ;
begin } -. procedure Z 1 ( ... ) ; begin
other declarations for the body of procedure Z.
} _ _ body of procedure 21. end 21 ;
X(Zl)
1
of the procedure ~ ~ statements bodyofZ.
j endZ;
X(2)
end ofprogram
1I -<--- statements of the program.
The Use of Recursive Procedure in ALGOL 60
45
In this way procedure X is called recursively though not being itself recursive. A realistic example is provided by the algorithm for the so-called Rombergintegration (Refs. 4, 10): real procedure rombergintegr(fct,lgr,rgr,ord) ; value lgr.rgr.ard ; reallgr, rgr; integer ord ; real procedurefct ;
begin }
For the details of the body of this procedure see Ref. I.
end
This procedure evaluates rsr
f fct(x) dx Igr
with an error term of order 2 X ord+2. However, the same procedure can be used to evaluate also double integrals. As an example the following piece of program describes the approximate computation of
with 289 mesh points: I: begin real proceduref(x) ; value x ; real x ; begin real procedure g(y) ; valuey ; realy ; comment g uses the non-local quantity real x ; 2: g: = sin (x x y) ; 3: f: = rombergintegr (g,0,I-x,4) endf; 4: J: = rombergintegr(j,0,1,4) end I ;
Explanation: During execution of procedure rombergintegr (of statement 4) the procedure f is called in from time to time in order to provide the integrand. The execution off causes (in statement 3) another call of romberg integrwhich computes the value 1-.-
fo
sin (xy)dy
for a prescribed x, and delivers it back to the rombergintegr (execution of statement 4).
H. Rutishauser
46
Similar programs for integrals in more dimensions get still more complicated since these need one extra procedure declaration for every dimension. It may be shown however, that a special arrangement allows the evaluation ofn-dimensional integrals for arbitrary n: If the domain is an n-dimensional parallelepipedon defined by the lower and upper bounds of the co-ordinates:
lb[k] ~ x[k] ~ ub[k] (for k = 1, 2, ..., n), and if the function to be integrated is given asfct(n,x), where x is the array xCI], x[2], .• " x[n], then such a program looks as follows: real procedure g(n,ftt,lb,ub,ord,k) ; value k,n ; integer n,ord, k ; array lb,ub ; real procedurefct ; begin own real array x[ I : n] ; real proceduref(v) ; real v ; ifk=nthen begin
x[n] : = v ; f: =jct(n,x) end else begin
x[k] : = v ; f: =.g(n,ftt,lb,ub,ord,k+I) end;
g: = rombergintegr(f,lb[k],ub[k],ord) endg
However, this construction is rather intricate and its feasibility for publication may be seriously doubted for two reasons. First, it requires that the implementation of own-arrays occurs according to a proposal of Sattley and Ingerman (Ref. 9) which is not specified by the ALGOL-report nor generally accepted. t Second, the function g thus defined cannot be used freely in so far as outside the procedure body ofg a function designator g(n,jct,lb,ub,ord,k) is only meaningful if k has the actual value 1. Indeed g(n,fct,lb,ub,ord,k) represents formally the value IIb[k]
f
lb[k]
lIb[k+ I] ..• "b[n]
I
I
jet(n,x) dx[k] . dx[k+l] ... dx[n]
lb[k+ 1] ••• [b[Il]
t To the present authors opinion, this proposal is better than a counter proposal which considers the own-variables of different recursion levels as different quantities. The latter makes communication between different recursion levels impossible.
The Use of Recursive Procedure in ALGOL 60
47
as function of the remaining variables x[l], x[2], ..., x[k-l], but the latter are undefined outside the procedure body of g. III. DIRECT RECURSIVITY
Realistic examples of procedures which call themselves without an intermediate procedure, seem to be very rare, at least in numerical applications. As an example we want to set up the algorithm for the determination of the smallest eigenvalue of a positive definite symmetric matrix by aid of the LR-transformation (Refs. 5, 6, 7). Let A be a positive definite symmetric matrix, preferably one which has non-zero elements only in the vicinity of the diagonal. Then starting with Ao = A and Zo = Z the following algorithm produces an infinite sequence A l , A2 , Ag, ••• of matrices such that the As zsI are all similar:
+
Choose j, } decompose A, - y/ into RrRst _ 0 1 2 T S , , , ••• Compute As+1 = R sRs compute Zs +1 = Zs Ys It can be enforced that the last diagonal element of As converges to zero and Zs to the smallest eigenvalue of A for s ~ co, but the speed of this convergence depends heavily upon the choice of the v., For every s this value should be chosen as large as possible, but must of course lie below the smallest eigenvalue of All because the Choleski method can be carried out in the real domain only for positive definite matrices. On the other hand, since the smallest eigenvalue of Asis yet unknown, it may occur thatj , gets too large and accordingly the Choleski decomposition will fail at somestage. Such a failure occurs as follows: After k elimination steps (we assume that k ~ 1) the matrix B = As - yslhas been decomposed partially:
+
U
B
=
v
1 - - - 1 -- - 1
w where pTp order k.
= U, PTQ = V, W -
Q into
o
w*
QTQ = W*, and Uis a square matrix of
t This is the so-called Choleski decomposition of a positive definite matrix into two transposed factors. This requires n square roots for the whole decomposition. For the solution of linear systems th ese square roots can be avoided, but not for the Lk-transformation.
48
H. Rutishauser Now the continuation of the process requires the Choleski decomposition
ofW* which begins with taking the square root of the first diagonal element ofW*. If this element is non-positive (zero could be tolerated fork = n - 1) the decomposition of As - y.! cannot be continued and must be repeated with a smaller value ofys. This repetition seems to be a loss of computing time, but it is actually a net gain. Indeed, it can be shown that the lowest eigenvalue ofthe n-k row matrix W* above isa lower bound for the spectrum of B = As - ysI, i.e, the lowest eigenvalue of W* Y.!* (where 1* is the n-k row unit matrix) is a lower bound for the lowest eigenvalue ofAs' Therefore, if we take this eigenvalue x as the new Ys' then the decomposition must succeed. Practical evidence showed that this x is often a quite narrow lower bound which may result in a tremendous acceleration of the convergence. However, the computation of x already requires the algorithm that we are going to devise, that is, this algorithm will be recursive. This recursion will not be a trivial one, since we cannot know in advance when a failure will occur, especially not ofwhat order the matrix W* will be. Furthermore the recursion will be finite, provided we enforce-by choosing j, < aI,I-that at least one elimination step can be carried out, and therefore the order of W* is at most n - 1. Practically however, we have to reject this way of computing the new Ys after a failure of decomposition in certain cases: First, iffailure takes place already at an early stage of the Choleski decomposition, such that the order of W* will be comparable to n, larger than nf3, say. Second, if W* y)* is not positive definite, then the new Ys would be negative; however, this would not be meaningful since already Ys = 0 is a better choice. In both cases the new value ofy must be furnished by other methods, probably again with the risk of'failure. But this is still better than to make the recursion chain too long. On the whole we obtain the following procedure:
+
+
. procedure reclr(n,a,eps,iriform) result: (a, lambda) exit: (indif) value n ; integer n ; real eps, lambda ; array a ; procedure inform ; label indef ; comm.ent reclr computes the lowest eigenvalue lambda of a positive definite symmetric matrix of order n with elements a[i,j] (i,j = 1, 2, ..., n). The elements of the matrix are usualIy changed by the process, except if the matrix is not positive definite, then an exit through label indef occurs. .inform is a procedure intended for interference from outside (see Re£ 8) and eps is a tolerance which determines the accuracy of lambda. The meaning of the procedures decomp, recomb, ww, adjust and tworow is explained below ;
The Use
of Recursive Procedure in ALGOL 60
49
begin
realy,R;,h,plzi ; integer k,s ; arrayr[l:n, 1:n],w[1:1+entier(nj3], 1:1+entier(nj3)J; procedure nil(a,b,c,d,e,f,g,h) result: (d,g) exit: (i) ; ; 1: iftz =l then begin lambda: = a[l,l] ; goto ex end; 2: if n=2 then begin lambda: = tworow(a) ; goto ex end; 3: s: = 0 ; h : = y : = z : = phi: = 0 ; 4: phi: = (1+phi)j2 ; y: = hXphi ; 5: inform(n,a,s,y,z,r,phi,h) result: (y,phi) exit: (ex) ; 6: decomp(n,a,y) result: (r,k) exit: (fail) ; 7: recomb(n,r) result: (a,h) ; 8: s:=s+l; z:= z+y; 9: if a[n,n] < eps then 10: begin lambda: = z + a[n,n] ; goto ex end ; 11: goto 4 ; 12: fail: if n-k > entier(nj3) 1 then goto mod; 13: ww(n,k,a,r) result : (w) ; 14: reelr(n-k,w,epsj4,nil) result: (w,y) exit: (mod) ;
+
15: 16:
y: = y - epsj2 ;
goto 6 ; 17: mod: adjust(n,k,a,y,phi) result: (y,phi) ; 18: goto 6 ; ex: end reelr In this program the procedures tworow, decomp, recomb, ww, adjust, have the following meaning (no further details are given for these procedures, since this is not relevant for the understanding of the recursive process): real procedure tworow(a) computes the lowest eigenvalue ofthe matrix with the elements a[l,l], a[I,2] =a[2,1], a[2,2]. procedure decomp(n,a,y) result: (r,k) exit: (fail) executes the Choleski decomposition of the matrix A, - yJ into RTR" where in the program the elements of the matrices A"R, and the variabley, are denoted by a[i,j], r[i,j] andy respectively. In case offailure of decomposition a jump to fail occurs, after which k denotes the number of successful elimination steps. procedure recomb (n,r) result: (a,h) computes the M atrix AI+ 1 = R,R'[. The elements of the matrix A, + 1 are again denoted by a[i,j]. At the same time also the square of the smallest diagonal element of R, is computed and denoted by h. procedure ww(n,k,a,r) result: (w) computes the elements of the matrix 4-
H. Rutishauser
50
+
w* y.I* of order n - k from the matrices a and r ,as far as they are available at the moment offailure. procedure adjust(n,k,a,y,phi) result: (y,phi) computes new values ofy and phi after failure of decomposition and it the recursion cannot be applied for one of the reasons mentioned in the text. This procedure decreases the' risk factor' phi and decreases alsoy at a rate depending on phi such that the whole computing process 'learns' whether it should try the decomposition with a smaller or larger value ofy (for details of a similar process see Ref. 7). procedure nil is a dummy procedure. IV. IMPLEMENTATION OF RECURSIVE PROCEDURES
The implementation ofrecursive procedures is usually considered difficult, furthermore it has been reported that the object programs produced by compilers which allow recursive procedures. are lessefficient. Ofcourse this is true to some extent, but the sayings are grossly exaggerated. t In any case the implementation problem for recursive procedures has been solved by Dijkstra (Ref. 3) and also by Sattley and Ingerman (Ref. 9), and accordingly ambitious compiler-builders will certainly include recursive procedures. But ifrestrictions have to bemade, then recursive procedures will be one of the first features to be omitted. This, however, willnot even entirely forbid their actual use. Indeed, in all realistic examples that came to the author's attention, the total depth ofrecursivity was rather low, and therefore such examples could still be treated non-recursively by providing a sufficient number of copies of the procedures in question for the translation process. In view of this possibility the decision of the ALCOR t group, not to implement the recursive procedures, is certainly well founded. One thing however should not be done. Ithas been proposed to simplify the implementation problem by declaring recursive procedures as such by an additional declarator recursive. This would be extremely unfortunate since it obviously would disallow indirect recursivity which contribute at least 50 per cent of all useful examples of recursive procedures. V.REFERENCES F. L., 'Algorithm 60'. Comm. A.C.M. 4, No.6, 255 (1961). R., 'ALGOL-Manual der ALCOR.Gruppe'. Elektron. Rechenanl, 3, 206212; 259-265 (1961) (to be continued).
1.
BAUER,
2.
BAUMANN,
t As an additional indication that the difficulties of implementation cannot be overwhelming it may be reported that the library routine system that was set up for the ERMETH in the pre-ALGOL days, automatically allowed recursive subroutines. t For the intentions of the ALCOR group see Ref. 2.
The Use of Recursive Procedure in ALGOL 60
51
3. DXJKSTRA, E. W., 'Recursive Programming'. Num. Math. 2, 312-318 (1960). 4. ROMBERG, W., 'Vereinfachte Numerische Integration'. Det. Kong. Norske Vidensk. Selsk. Forh. 28, No.7, Trondheim, 1955. 5. RUTISHAUSER, H., 'Solution of Eigenvalue Problems with the LR-Transformation.' Nat. Bur. Standards, Applied Mathematics Series, 49,47-81 (Washington, 1958). 6. RUTISHAUSER, H., 'Ueber eine kubisch konvergente Variante der LR-Transformation ', ZAMM, 40, 49-54 (1960). 7. RUTISHAUSER, H. and SCHWARZ, H. R., 'The LR-Transformation for Symmetric Matrices'. To appear in Num. Math. (1962). 8. RUTISHAUSER, H., 'Interference with an ALGOL-Procedure'. Annual Review in Automatic Programming, Vol. 2, pp. 67-75. Pergamon Press, Oxford (1961). 9. SATTLEY, K., and INGERMAN, P. Z., 'The Al1ocation ofStorage for Arrays in ALGOL60'. Internal progress report, Office of Computer Research and Education, University of Pennsylvania, Nov. 1960. 10. STIEFEL, E., 'Altes und Neues tiber Numerische Quadratur'. ZAMM,41 (1961) p. 408-413.
INCREASING emphasis has since recently been given to recursive procedures, i.e, procedures that call themselves. The present comments may be considered as a contribution to the debate that arose over this subject.
I. INTRODUCTION
Unfortunately recursive procedures have been advertised by examples such as real procedurefact(n) ; value n ; integer n ; fact: = if n=O then I else n Xfact(n-l) ;
or real procedure hornerl(a,k,z,n) ; value n,z,k ; real z ; integer k.n ; array a ; hornerl: = if k=n then a[n] else(z X hornerl(a,k+l,z,n)+a[kJ)
True, these examples look quite elegant and may serve to exemplify recursive procedures, butifa computing process canjust as well be described by a simple induction loop, then transforming it into a recursive procedure is sheer extravagance and might well discredit the concept of recursivity altogether. Of course, the solution of many problems present themselves directly in recursive form, but this should not be used for the numerical solution unless it can hardly be avoided. The present author sees mainly two classes of useful applications of recursive procedures. II. INDIRECT RECURSIVITY
Let Xbe an already written non-recursive procedure, one parameter of which represents itself a procedure, e.g, 4-3
H. Rutishauser procedure X(Y) ; procedure Y ;
begin
Y( ... ) ;
_ _ call of procedure Y
end X Now the user of procedure X has to provide the actual parameter corresponding to the formal parameter Y, i.e. he has to declare in his program a procedure Z as the actual counterpart ofY. The body of Z may again call for execution of X as follows:
Program: begin
}--
.
other declarations for the program
procedure Z ( ... ) ;
begin } -. procedure Z 1 ( ... ) ; begin
other declarations for the body of procedure Z.
} _ _ body of procedure 21. end 21 ;
X(Zl)
1
of the procedure ~ ~ statements bodyofZ.
j endZ;
X(2)
end ofprogram
1I -<--- statements of the program.
The Use of Recursive Procedure in ALGOL 60
45
In this way procedure X is called recursively though not being itself recursive. A realistic example is provided by the algorithm for the so-called Rombergintegration (Refs. 4, 10): real procedure rombergintegr(fct,lgr,rgr,ord) ; value lgr.rgr.ard ; reallgr, rgr; integer ord ; real procedurefct ;
begin }
For the details of the body of this procedure see Ref. I.
end
This procedure evaluates rsr
f fct(x) dx Igr
with an error term of order 2 X ord+2. However, the same procedure can be used to evaluate also double integrals. As an example the following piece of program describes the approximate computation of
with 289 mesh points: I: begin real proceduref(x) ; value x ; real x ; begin real procedure g(y) ; valuey ; realy ; comment g uses the non-local quantity real x ; 2: g: = sin (x x y) ; 3: f: = rombergintegr (g,0,I-x,4) endf; 4: J: = rombergintegr(j,0,1,4) end I ;
Explanation: During execution of procedure rombergintegr (of statement 4) the procedure f is called in from time to time in order to provide the integrand. The execution off causes (in statement 3) another call of romberg integrwhich computes the value 1-.-
fo
sin (xy)dy
for a prescribed x, and delivers it back to the rombergintegr (execution of statement 4).
H. Rutishauser
46
Similar programs for integrals in more dimensions get still more complicated since these need one extra procedure declaration for every dimension. It may be shown however, that a special arrangement allows the evaluation ofn-dimensional integrals for arbitrary n: If the domain is an n-dimensional parallelepipedon defined by the lower and upper bounds of the co-ordinates:
lb[k] ~ x[k] ~ ub[k] (for k = 1, 2, ..., n), and if the function to be integrated is given asfct(n,x), where x is the array xCI], x[2], .• " x[n], then such a program looks as follows: real procedure g(n,ftt,lb,ub,ord,k) ; value k,n ; integer n,ord, k ; array lb,ub ; real procedurefct ; begin own real array x[ I : n] ; real proceduref(v) ; real v ; ifk=nthen begin
x[n] : = v ; f: =jct(n,x) end else begin
x[k] : = v ; f: =.g(n,ftt,lb,ub,ord,k+I) end;
g: = rombergintegr(f,lb[k],ub[k],ord) endg
However, this construction is rather intricate and its feasibility for publication may be seriously doubted for two reasons. First, it requires that the implementation of own-arrays occurs according to a proposal of Sattley and Ingerman (Ref. 9) which is not specified by the ALGOL-report nor generally accepted. t Second, the function g thus defined cannot be used freely in so far as outside the procedure body ofg a function designator g(n,jct,lb,ub,ord,k) is only meaningful if k has the actual value 1. Indeed g(n,fct,lb,ub,ord,k) represents formally the value IIb[k]
f
lb[k]
lIb[k+ I] ..• "b[n]
I
I
jet(n,x) dx[k] . dx[k+l] ... dx[n]
lb[k+ 1] ••• [b[Il]
t To the present authors opinion, this proposal is better than a counter proposal which considers the own-variables of different recursion levels as different quantities. The latter makes communication between different recursion levels impossible.
The Use of Recursive Procedure in ALGOL 60
47
as function of the remaining variables x[l], x[2], ..., x[k-l], but the latter are undefined outside the procedure body of g. III. DIRECT RECURSIVITY
Realistic examples of procedures which call themselves without an intermediate procedure, seem to be very rare, at least in numerical applications. As an example we want to set up the algorithm for the determination of the smallest eigenvalue of a positive definite symmetric matrix by aid of the LR-transformation (Refs. 5, 6, 7). Let A be a positive definite symmetric matrix, preferably one which has non-zero elements only in the vicinity of the diagonal. Then starting with Ao = A and Zo = Z the following algorithm produces an infinite sequence A l , A2 , Ag, ••• of matrices such that the As zsI are all similar:
+
Choose j, } decompose A, - y/ into RrRst _ 0 1 2 T S , , , ••• Compute As+1 = R sRs compute Zs +1 = Zs Ys It can be enforced that the last diagonal element of As converges to zero and Zs to the smallest eigenvalue of A for s ~ co, but the speed of this convergence depends heavily upon the choice of the v., For every s this value should be chosen as large as possible, but must of course lie below the smallest eigenvalue of All because the Choleski method can be carried out in the real domain only for positive definite matrices. On the other hand, since the smallest eigenvalue of Asis yet unknown, it may occur thatj , gets too large and accordingly the Choleski decomposition will fail at somestage. Such a failure occurs as follows: After k elimination steps (we assume that k ~ 1) the matrix B = As - yslhas been decomposed partially:
+
U
B
=
v
1 - - - 1 -- - 1
w where pTp order k.
= U, PTQ = V, W -
Q into
o
w*
QTQ = W*, and Uis a square matrix of
t This is the so-called Choleski decomposition of a positive definite matrix into two transposed factors. This requires n square roots for the whole decomposition. For the solution of linear systems th ese square roots can be avoided, but not for the Lk-transformation.
48
H. Rutishauser Now the continuation of the process requires the Choleski decomposition
ofW* which begins with taking the square root of the first diagonal element ofW*. If this element is non-positive (zero could be tolerated fork = n - 1) the decomposition of As - y.! cannot be continued and must be repeated with a smaller value ofys. This repetition seems to be a loss of computing time, but it is actually a net gain. Indeed, it can be shown that the lowest eigenvalue ofthe n-k row matrix W* above isa lower bound for the spectrum of B = As - ysI, i.e, the lowest eigenvalue of W* Y.!* (where 1* is the n-k row unit matrix) is a lower bound for the lowest eigenvalue ofAs' Therefore, if we take this eigenvalue x as the new Ys' then the decomposition must succeed. Practical evidence showed that this x is often a quite narrow lower bound which may result in a tremendous acceleration of the convergence. However, the computation of x already requires the algorithm that we are going to devise, that is, this algorithm will be recursive. This recursion will not be a trivial one, since we cannot know in advance when a failure will occur, especially not ofwhat order the matrix W* will be. Furthermore the recursion will be finite, provided we enforce-by choosing j, < aI,I-that at least one elimination step can be carried out, and therefore the order of W* is at most n - 1. Practically however, we have to reject this way of computing the new Ys after a failure of decomposition in certain cases: First, iffailure takes place already at an early stage of the Choleski decomposition, such that the order of W* will be comparable to n, larger than nf3, say. Second, if W* y)* is not positive definite, then the new Ys would be negative; however, this would not be meaningful since already Ys = 0 is a better choice. In both cases the new value ofy must be furnished by other methods, probably again with the risk of'failure. But this is still better than to make the recursion chain too long. On the whole we obtain the following procedure:
+
+
. procedure reclr(n,a,eps,iriform) result: (a, lambda) exit: (indif) value n ; integer n ; real eps, lambda ; array a ; procedure inform ; label indef ; comm.ent reclr computes the lowest eigenvalue lambda of a positive definite symmetric matrix of order n with elements a[i,j] (i,j = 1, 2, ..., n). The elements of the matrix are usualIy changed by the process, except if the matrix is not positive definite, then an exit through label indef occurs. .inform is a procedure intended for interference from outside (see Re£ 8) and eps is a tolerance which determines the accuracy of lambda. The meaning of the procedures decomp, recomb, ww, adjust and tworow is explained below ;
The Use
of Recursive Procedure in ALGOL 60
49
begin
realy,R;,h,plzi ; integer k,s ; arrayr[l:n, 1:n],w[1:1+entier(nj3], 1:1+entier(nj3)J; procedure nil(a,b,c,d,e,f,g,h) result: (d,g) exit: (i) ; ; 1: iftz =l then begin lambda: = a[l,l] ; goto ex end; 2: if n=2 then begin lambda: = tworow(a) ; goto ex end; 3: s: = 0 ; h : = y : = z : = phi: = 0 ; 4: phi: = (1+phi)j2 ; y: = hXphi ; 5: inform(n,a,s,y,z,r,phi,h) result: (y,phi) exit: (ex) ; 6: decomp(n,a,y) result: (r,k) exit: (fail) ; 7: recomb(n,r) result: (a,h) ; 8: s:=s+l; z:= z+y; 9: if a[n,n] < eps then 10: begin lambda: = z + a[n,n] ; goto ex end ; 11: goto 4 ; 12: fail: if n-k > entier(nj3) 1 then goto mod; 13: ww(n,k,a,r) result : (w) ; 14: reelr(n-k,w,epsj4,nil) result: (w,y) exit: (mod) ;
+
15: 16:
y: = y - epsj2 ;
goto 6 ; 17: mod: adjust(n,k,a,y,phi) result: (y,phi) ; 18: goto 6 ; ex: end reelr In this program the procedures tworow, decomp, recomb, ww, adjust, have the following meaning (no further details are given for these procedures, since this is not relevant for the understanding of the recursive process): real procedure tworow(a) computes the lowest eigenvalue ofthe matrix with the elements a[l,l], a[I,2] =a[2,1], a[2,2]. procedure decomp(n,a,y) result: (r,k) exit: (fail) executes the Choleski decomposition of the matrix A, - yJ into RTR" where in the program the elements of the matrices A"R, and the variabley, are denoted by a[i,j], r[i,j] andy respectively. In case offailure of decomposition a jump to fail occurs, after which k denotes the number of successful elimination steps. procedure recomb (n,r) result: (a,h) computes the M atrix AI+ 1 = R,R'[. The elements of the matrix A, + 1 are again denoted by a[i,j]. At the same time also the square of the smallest diagonal element of R, is computed and denoted by h. procedure ww(n,k,a,r) result: (w) computes the elements of the matrix 4-
H. Rutishauser
50
+
w* y.I* of order n - k from the matrices a and r ,as far as they are available at the moment offailure. procedure adjust(n,k,a,y,phi) result: (y,phi) computes new values ofy and phi after failure of decomposition and it the recursion cannot be applied for one of the reasons mentioned in the text. This procedure decreases the' risk factor' phi and decreases alsoy at a rate depending on phi such that the whole computing process 'learns' whether it should try the decomposition with a smaller or larger value ofy (for details of a similar process see Ref. 7). procedure nil is a dummy procedure. IV. IMPLEMENTATION OF RECURSIVE PROCEDURES
The implementation ofrecursive procedures is usually considered difficult, furthermore it has been reported that the object programs produced by compilers which allow recursive procedures. are lessefficient. Ofcourse this is true to some extent, but the sayings are grossly exaggerated. t In any case the implementation problem for recursive procedures has been solved by Dijkstra (Ref. 3) and also by Sattley and Ingerman (Ref. 9), and accordingly ambitious compiler-builders will certainly include recursive procedures. But ifrestrictions have to bemade, then recursive procedures will be one of the first features to be omitted. This, however, willnot even entirely forbid their actual use. Indeed, in all realistic examples that came to the author's attention, the total depth ofrecursivity was rather low, and therefore such examples could still be treated non-recursively by providing a sufficient number of copies of the procedures in question for the translation process. In view of this possibility the decision of the ALCOR t group, not to implement the recursive procedures, is certainly well founded. One thing however should not be done. Ithas been proposed to simplify the implementation problem by declaring recursive procedures as such by an additional declarator recursive. This would be extremely unfortunate since it obviously would disallow indirect recursivity which contribute at least 50 per cent of all useful examples of recursive procedures. V.REFERENCES F. L., 'Algorithm 60'. Comm. A.C.M. 4, No.6, 255 (1961). R., 'ALGOL-Manual der ALCOR.Gruppe'. Elektron. Rechenanl, 3, 206212; 259-265 (1961) (to be continued).
1.
BAUER,
2.
BAUMANN,
t As an additional indication that the difficulties of implementation cannot be overwhelming it may be reported that the library routine system that was set up for the ERMETH in the pre-ALGOL days, automatically allowed recursive subroutines. t For the intentions of the ALCOR group see Ref. 2.
The Use of Recursive Procedure in ALGOL 60
51
3. DXJKSTRA, E. W., 'Recursive Programming'. Num. Math. 2, 312-318 (1960). 4. ROMBERG, W., 'Vereinfachte Numerische Integration'. Det. Kong. Norske Vidensk. Selsk. Forh. 28, No.7, Trondheim, 1955. 5. RUTISHAUSER, H., 'Solution of Eigenvalue Problems with the LR-Transformation.' Nat. Bur. Standards, Applied Mathematics Series, 49,47-81 (Washington, 1958). 6. RUTISHAUSER, H., 'Ueber eine kubisch konvergente Variante der LR-Transformation ', ZAMM, 40, 49-54 (1960). 7. RUTISHAUSER, H. and SCHWARZ, H. R., 'The LR-Transformation for Symmetric Matrices'. To appear in Num. Math. (1962). 8. RUTISHAUSER, H., 'Interference with an ALGOL-Procedure'. Annual Review in Automatic Programming, Vol. 2, pp. 67-75. Pergamon Press, Oxford (1961). 9. SATTLEY, K., and INGERMAN, P. Z., 'The Al1ocation ofStorage for Arrays in ALGOL60'. Internal progress report, Office of Computer Research and Education, University of Pennsylvania, Nov. 1960. 10. STIEFEL, E., 'Altes und Neues tiber Numerische Quadratur'. ZAMM,41 (1961) p. 408-413.