The Algol string quotes used in the original document have been replaced by the " symbol.
The Library of Algol procedures currently available on Atlas is stored on the Standard Algol Library Tape ALGOLIB. Full details of how to incorporate these library procedures into an Algol program are given in the Atlas Algol System manual. A subset of the more widely used procedures are available on the disc area R502 ALGOLIB. In the following index, library names followed by asterisk indicate that the item is also available on the disc area. The library is divided into three main groups SC, ACM and ICT.
The SC group are procedures specifically designed to run with the GROATS graphic output system. Full details of these can be found in the GROATS manual.
The Algorithms available under the Class Name ACM are those published in the Communications of the ACM. The numbering of the Algorithms is the same as given in the CACM. These Algorithms have been punched at the Atlas Laboratory. Corrections given in Certifications have been added. Trial compilations have been used to give a syntactic check of the Algorithms. No semantic checking of the Algorithms has been attempted. The definitions of these Algorithms can either be found in the relevant copy of the Communications or alternatively a complete set of definitions is available from the ACM. A copy of this is available at the Atlas Laboratory.
The class name ICT has been used for all library items that do not fall into either of the other two classes. Complete details of these are given in the following sections. These procedures should be fully debugged and working.
1 * dump program 2 * release mag tape 3 * dump 4 * print input 5 Array handling procedures 6 Array handling procedures 7 Matrix procedures 8 Matrix procedures 9 Matrix procedures 10 Matrix procedures 11 Matrix procedures 12 * set up directory 13 * break output 14 * random, initialise generator, extract generator 15 * Upper case versions of SIN, COS etc 16 * ham 17 * dump tape positions 18 * whats on my library tape 19 * erase 20 * KDF9 magnetic tape handling procedures 21 22 writekdf9, format 23 char in kdf9, char out kdf9 24 copy text kdf9 25 read array kdf9 26 write array kdf9 27 address, lowbound, range, size 28 advance elliott, buffer elliott 29 * dialect list 30 * initialise dialect : ATLAS CARDS 31 * initialise dialect : ICT 1900 TAPE 32 33 34 35 36 37 * outjob title time date 38 chain 39 * call algol 40 * dump program 41 * lower triangle iliffe vectors 42 * interpret 43 * Fixed block magnetic tape procedures 44 * fast char print routine 45 * chain 46 * elegant output 47 * layout, inlogical, outlogical 48 * elapsed time 49 50 * gresolpiv 51 * tridiagreson 52 * deter 53 * trace 54 * pivotmax 55 * detsym 1, solsym 1, detinvsym 1 56 * detbnd, solbnd 57 * Least square problems 58 * pseudoinv 59 * pseudo inverse 60 * eigenvalues 61 * iter jacobi 62 * jacobi 63 * rutisvalues 64 * hyman values 65 * eigenvec 66 67 68 69 70 * newton complexe 71 * muller 72 * laguerre 73 * bissection 74 * bisdou 75 * rungekutta 76 * rungekutta 77 * rungekutta 78 79 80 * insire 81 * in dourec 82 * int3neville 83 * intcossin 84 85 * approxicon 86 * remez 87 * tchebecha 88 * chebfit 89 * tchebdessous 90 * tchfbing 91 * mcpolysp 92 * methdintepoly 93 * coef spline trois 94 * spline hermite 95 * splderivee 96 * classemarkoff 97 * sous classe cyclique 98 * polweyl 99 * not lag 100* reduction 101* ratqr 102* tred 103* tridi inverse iteration 104* back transformation 105* bisect 106* Character manipulation procedures 107* Lineprinter graphical procedures 108* random 109* chisqaure 110* timer
procedure dump program (k,m,n);
This procedure is designed to dump a fully compiled program in binary on magnetic tape. It is intended to serve the dual purpose of:-
This procedure when called will dump all the blocks in use at that point and all the B-lines onto magnetic tape (logical number k) in blocks m to n of the tape. If the dumping area on the magnetic tape is not sufficient, the diagnostic:-
DUMP SPACE EXCEEDED
will appear and the program will end. However if there is sufficient space the following line will appear:-
DUMP TO TAPE k (Title of tape from Job Description) IN THE RANGE (m,n)
Once the dumping has taken place the program will carry on executing from the point where the procedure dump program was called. There will now be a compiled version of the program on the magnetic tape and the following steering tape, substituting for k and n the same values as above, should be used for subsequent runs:-
JOB ---------- TAPE k --- ---------- COMPILER ABL 1001, k, 0, n 1002, k, 0, 2: 121, 127, 0, 2: EA ***Z
This will bring the program down from magnetic tape and re-enter the program immediately after the call of dump program. An example of the use of this procedure is as follows:-
JOB 10000 HOPGOOD DUMP AFTER READING FIXED DATA
OUTPUT 0 LINEPRINTER 1000 LINES
INPUT 1 CONSTANTDATA
INPUT 2 VARIABLEDATA
TAPE 3 N0009DMP*PERMIT
DISC 99 R502 ALGOLIB
COMPUTING 30 SECONDS
COMPILER ALGOL
INPUT INTERNAL ICT WITH ICT I/O PROCEDURES;
begin
library ICT1;
select input (1);
comment Read in all data which is constant for all runs
of the program and set up all initial conditions;
dump program (3,1,50);
select input (2);
comment Read in variable data and carry on with program;
end
***Z
Note that by inserting data reading instructions in place of the first comment any large amount of fixed data could be stored on magnetic tape with the program. Subsequent runs would use:-
JOB I0000 HOPGOOD READ PROGRAM FROM TAPE OUTPUT 0 LINEPRINTER 1000 LINES INPUT 2 VARIABLEDATA TAPE 3 N0009DMP*INHIBIT COMPUTING 30 SECONDS COMPILER ABL 1001, 3, 0, 50 1002, 3, 0, 2: 121, 127, 0, 2: EA ***Z
See ICT2 in connection with this procedure.
procedure release mag tape (k);
Once a magnetic tape has been finished with it is a help to the operating staff and the efficient running of the computer for the tape to be dismounted by the user. The deck and tape are then available for other users. This procedure could be used in connection with the procedure ICT1 and also the library tape ALGOLIB.
The procedure will release the magnetic tape with logical tape number k.
procedure dump (str);
string str;
This procedure will dump either in octal or in machine order code format certain sections of the Algol program during program execution. It is of use mainly to users having a good knowledge of the internal workings of the Algol Compiler.
The string argument should consist of a set of words separated by commas. The possible words are:-
directory program stack constants i/o
For example:-
dump ('directory, stack');
The areas of store dumped are as follows:-
Each line has its store address in octal printed before it. If several lines are identical the first only is printed. A blank line is inserted to show where lines have not been printed.
procedure print input;
This procedure, when called, will read and print all the remaining characters on the selected input stream. The input is read using the same code conversion as data on that stream. It is designed to give a listing of program or data tapes. For example:-
JOB I0000 HOPGOOD LIST DATA TAPE OUTPUT 0 LINEPRINTER 500 LINES COMPUTING 30 SECONDS TAPE 99 ALGOLIB*INHIBIT COMPILER ALGOL INPUT INTERNAL ICT WITH ICT I/O PROCEDURES; begin library ICT4; print input; end ***T
followed by a tape ending on ***Z would cause the tape to be printed. Some effort is made to keep paperthrows at positions where additional newlines have been inserted.
The following set of procedures can be used for general Algol arrays. The notation used is as follows:-
A,B,C are real arrays, B and C (where used) must have the same number of dimensions as A and corresponding subscripts must have the same number of possible values. For example:-
A [1:3, 10:13, 19:25], B [2 :4, -7:-4, 6:12], C [0:2, 3:6, -7:-1]
is correct, but
A[1:2, 1:2], B[1:2, 1:2], C[1:2, 1:3] and A [1:2, 1:1], B[1:2, 1:1], C[1:2 ]
are incorrect.
If the arrays substituted for A, B, C do not satisfy the above restrictions, there will be a run-time fault. Any or all of the arrays substituted for A, B, C may be the same. Declarations for A, B, C are omitted throughout, in this document, but are of course included in the actual procedures.
The procedures have code bodies.
integer procedure dim(A);
Sets dim = dimension of the array A.
integer procedure bound (A,n); value n;integer n;
Sets bound = nth, bound of the array A. This bound will be the upper or lower bound of the mth subscript according as n = 2m -1,2m . If the nth bound does not exist, there will be a run-time fault.
procedure readarray(A);
Reads data and stores it by rows in A.
procedure printarray (A,m,n); value m,n; integer m,n;
Outputs the a-dimensional array A by rows in style m,n, inserting (2i+1) blank lines when the (a-i)th subscript is incremented.
procedure xarray(A,x); value x; real x;
Sets each element of A equal to x
procedure addxarray (A, B, x, C) ;value x; real x;
Sets A= B + x*C. Efficient code is obeyed for the special cases:
x = 0 giving A = B x = 1 giving A = B + C C ≠ B, x = -1 giving A = B - C C = B, x = -1 giving A = null array C = B, x = -2 giving A = -B C = B giving A = (1 + x) B
A set of additional procedures which call procedures in the ICT5 set:-
procedure nullarray (A);
Sets A = null array = xarray (A,0);
procedure equarray (A,B);
Sets A = B = addxarray (A,B,0,B);
procedure negarray (A,B);
Sets A = -B = addxarray (A,B,-2,B);
procedure addarray (A,B,C);
Sets A=B+C = addxarray (A,B,1,C);
procedure subarray (A,B,C) ;
Sets A = B-C = addxarray (A,B,-1,C);
procedure outarray (A) ;
Outputs the array A = printarray (A,0,10);
ICT7 to ICT11 define a set of procedures which operate on matrices. The notation used is as follows:-
A, B, C are matrices of sizes (a1 × a2), (b1 × b2), (c1 × c2) respectively. D is a diagonal matrix of size (d × d) stored as a vector of length d, I is the unit matrix (with size determined by context).
If A, B, C are not matrices, D is not a vector, or any of the restrictions specified for each procedure are not satisfied, then there will be a run-time fault. Any or all of the matrices substituted for A, B, C may be the same unless stated otherwise. Declarations for A, B, C, D are omitted throughout.
procedure transmx(A);
Restrictions: a1 = a2
Sets A = transpose (A).
procedure equtransmx(A,B);
Restrictions: a1 = b2, a2 = b1
Sets B = transpose (A)
procedure xunitmx(A,x); value x; real x;
Restrictions: a1 = a2
Sets A = x*I.
procedure addxunitmx(A,B,x); value x; real x;
Restrictions: a1 = a2 = b1 = b2
Sets A = B + x * I.
Efficient code is obeyed in the special case A = B giving A = A + x * I
procedure unitmx(A);
Restrictions: a1 = a2
Sets A = I = xunitmx(A,1);
procedure addunitmx(A,B);
Restrictions: a1 = a2 = b1 = b2
Sets A = B + I = addxunitmx(A,B,1);
procedure subunitmx(A,B);
Restrictions: a1 = a2 = b1 = b2
Sets A= B - I = addxunitmx(A,B,-1);
procedure gendiagmx(A,B,D,x, type);
value x, type; real x; integer type;
Restrictions: a1 = a2 = b1 = b2 = d
If type = 0, sets D = x * (leading diagonal of A)
B is not used and should be set as A.
If type ≠ 0, sets A = B + x * D.
Efficient code is obeyed for the special cases:-
type ≠ 0, A = B, x = 1 giving A = A + D type ≠ 0, A = B, x = -1 giving A = A - D type ≠ 0, A = B giving A = A + x*D type = 0, x = 1 giving D = (l. d. of A) type = 0, x = -1 giving D = -(l. d. of A)
procedure addiagmx(A,B,D);
Restrictions: a1 = a2 = b1 = b2 = d
Sets A = B + D = gendiagmx(A,B,D,1,1);
procedure subdiagmx(A,B,D);
Restrictions: a1 = a2 = b1 = b2 = d
Sets A = B - D = gendiamx(A,B,D,-1,1);
procedure equidiagmx(D,A);
Restrictions: a1 = a2 = d
Sets D = (leading diagonal of A) = gendiagmx(A,A,D,1,0);
procedure negdiagmx(D,A);
Restrictions: a1 = a2 = d
Sets D = (leading diagonal of A) = gendiagmx(A,A,D,-1,0);
procedure addxdiagmx(A,B,x,D); real x;
Restrictions: a1 = a2 = b1 = b2 =d
Sets A = B + x D = gendiagmx(A,B,D,x,1);
procedure xdiagmx(D,x,A); real x;
Restrictions: a1 = a2 = b1 = b2 =d
Sets D=x * (leading diagonal of A) = gendlagmx(A,A,D,x,0);
procedure multmx(A,B,C);
Restrictions: A, B, C different a1=b1 a2=c2 b2=c1
Sets A = B*C = genmultmx(A,B,C,"ex");
procedure multtransmx(A,B,C);
Restrictions: A, B, C different a1=b1 a2=c1 b2=x2
Sets A = B* transpose(C) = genmultmx(A,B,C,"ext");
procedure transmultmx(A,B,C);
Restrictions: A, B, C different a1=b2 a2=c2 1=c1
Sets A = transpose (B)* C = genmultmx(A,B,C, "etx");
procedure genmultmx(A,B,C,s); string s;
Restrictions: A must be different from both B and C, but B and C may be the same
The operation is determined by the string s which must be one of the following:
ex giving A = B*C Restriction: a1=b1 a2=c2 b2=c1 ex giving A = -B*C Restriction: a1=b1 a2=c2 b2=c1 px giving A = A + B*C Restriction: a1=b1 a2=c2 b2=c1 mx glvlng A = A - B*C Restriction: a1=b1 a2=c2 b2=c1 ext giving A = B*transpose(C) Restriction: a1=b1 a2=c1 b2=c2 ext giving A = -B*transpose(C) Restriction: a1=b1 a2=c1 b2=c2 pxt giving A = A + B*transpose(C) Restriction: a1=b1 a2=c1 b2=c2 mxt giving A = A- B*transpose( C) Restriction: a1=b1 a2=c1 b2=c2 etx giving A = transpose(B)*C Restriction: a1=b2 a2=c2 b1=c1 etx giving A = -transpose (B)*C Restriction: a1=b2 a2=c2 b1=c1 ptx giving A = A + transpose(B)*C Restriction: a1=b2 a2=c2 b1=c1 mtx giving A = A - transpose(B)*C Restriction: a1=b2 a2=c2 b1=c1
if s is not any of the above, then the action is undefined.
A set of procedures for doing matrix inversion, division and finding determinant.
real procedure invert (A,eps);
value eps; real eps;
Restrictions: a1 = a2
Inverts the matrix A on itself by Gauss-Jordan elimination. If, during elimination, the absolute value of any pivot is eps, there will be a run-time fault. If the inversion is successful, invert will contain the determinant of the original matrix.
real procedure invert part(A,p,q,n,eps);
value p,q,n,eps;
integer p,q,n; real eps;
Restrictions: a1≤p a1≤q a2≥p+n-1 a2≥q+n-1
This procedure acts in the same way as invert except that it inverts the submatrix of A given by A[i,j], p≤i≤p+m-1, q≤j≤q+n-1.
procedure invmx;
This procedure is required for use with invert and/or invertpart.
real procedure matdiv(A,B,eps);
value eps; real eps;
Restrictions: a1 = a2 = b1
Uses a Gauss elimination method to set B = inverse(A)*B.
If, during elimination, the absolute value of any pivot is eps, there will be a run-time fault. Otherwise matdiv will be set to the determinant of A. B may be a matrix or vector. A is destroyed.
real procedure det(A,eps);
value eps; real eps;
Uses a Gauss elimination method to set det to the determinant of A. If, during elimination, the absolute value of any pivot is ≤ eps, then det = 0. A is destroyed.
procedure divmx;
This procedure is required for use with matdiv and/or det.
procedure unfault invmx;
Alters procedure invmx so that if invertor invertpart attempts to invert a singular matrix (i.e. the absolute value of some pivot is ≤ eps), then instead of giving a run-time fault, the procedure exits normally and invert or invertpart is set = 0.
procedure unfault divmx;
Alters procedure divmx so that if matdiv attempts to divide by a singular matrix (i.e. the absolute value of some pivot is ≤ eps), then instead of giving a run-time fault, the procedure exits normally and 'matdiv' is set = 0.
If invmx or divmx is altered by the appropriate above procedure, it cannot be restored to its original condition.
procedure set up directory (str);
string str;
This procedure has been provided for users requiring to set up directories on private library tapes. All details on how to use this procedure are provided in Chapter 13 of the Atlas Algol System Manual. The form of the directory is as follows:-
WORD
0 First free block on tape : Tape number when last used
1 ≠0 for being used when in store : Position of next class name
relative to this one.
2 Length of this class name : First letter of class name.
3
4
.
.
. Last letter of class name : Library item 0 in this class
Library item 1 in this class etc.
procedure break output (n);
value n; integer n;
This procedure when called will break the output on stream n. If large quantities of output on one stream are required by a user it is important that this procedure is used. No output is actually sent to a peripheral until either a complete output stream has been produced or a break in the output stream is made. Consequently, for large amounts of output, if break output is not used then the output disc area will have to carry a large number of blocks of the uncompleted output causing inefficiencies in the system. Also, if a machine failure occurs, the whole of the output may be lost unless break output is used.
The procedure break output should be used about every 1000 lines.
procedure procedure random;
procedure initialise generator (left, right);
value left, right; integer left, right;
procedure extract generator (left,right);
integer left, right;
These three procedures have been implemented to produce, from successive calls of random a sequence of Pseudo Random Numbers chosen from a population having a rectangular distribution and the range 0 to 1. The procedures ensure that different samples are available for re-runs of any given job and also for repeating any particular sequence of Pseudo Random Numbers.
A generator word G for defining the start of the sequence of Pseudo Random Numbers is initialised by a call of initialise generator. The two arguments left and right may be any integer in the range 0 to 524287. They provide respectively the left and right 19 bits of the central 38 bits of the generator word G.
The generator word is reset as follows:
G(new) = ([G(old)*8111**3] modulo 2**39)* 2 ** -39
Each call of random will reset the generator word G and extract the next Pseudo Random Number from G.
If the components of the generator word G are required to be saved at any time then a call of extract generator will store the 38 relevant bits in left and right. Using these values as arguments to initialise generator will cause the sequence to be continued from the point left off.
An example of using these procedures would be:-
JOB
I0000 HOPGOOD RANDOM
OUTPUT 0 LINEPRINTER 200 LINES
COMPUTING 30 SECONDS
DISC 99 R502 ALGOLIB
COMPILER ALGOL
INPUT INTERNAL ICT WITH ICT I/O PROCEDURES;
begin
library ICT14;
integer i, j ,k;
initialise generator (349525, 74762);
for i := 1 step 1 until 100 do
begin
print (random, 1, 8);
newline (1);
end;
extract generator (j,k);
print (j, 5, 0); print (k, 6, 0);
end;
***Z
This contains the definitions or the standard Algol functions with abs replaced by ABS, sin replaced by SIN etc. These procedures are used with programs written for entry 20 or the ICL 1900 compiler. More details are given in the section (ix) or Chapter 4.3 or the Atlas Algol System manual.
procedure ham (DIRECTORY BLOCK, NAME, NUMBER, TITLE);
value DIRECTORY BLOCK, NUMBER;
integer DIRECTORY BLOCK, NUMBER;
string NAME, TITLE;
This outputs the library item stored on tape 99. The directory block is at block DIRECTORY BLOCK and the name or the library item is NAME. Its position is NUMBER. The library item is output on cards on stream 5. Each Algol symbol takes up one column between columns 1 and 72. Then follows the first six letters or TITLE in columns 73 to 78. The last two columns are used to serialise the cards from 0 upwards.
Full details or the way the Algol symbols are punched are given in Chapter 16 of the Atlas Algol System manual.
procedure dump tape positions;
This procedure can be used with the procedure dump program. It prints out the positions or the magnetic tapes when the procedure is called. The information for each tape is printed out a line at a time in ascending order or tape number. For each tape, the block position on the tape is output followed by the word position. Both are printed in octal.
procedure whats on my library tape (DIRECTORY, STREAM, PACKUP,
DIRECTORYNEW, STREAMNEW);
value DIRECTORY, STREAM, DIRECTORYNEW, STREAMNEW;
integer DIRECTORY, STREAM, DIRECTORYNEW, STREAMNEW;
Boolean PACKUP;
The design of the Algol library tape system is extremely simple. If a library item increases in size then, on redefinition, the whole item is moved to the next free blocks and a gap is left in the position where it was originally. If this happens several times then it could cause much longer searches of the magnetic tape than are necessary. The aim of the above procedure is to give the user a map showing where his library items occur and to remove the gaps if he so desires.
To list the positions and size of library items on the library only the first 3 arguments are relevant. The two remaining arguments can be set arbitrarily. The Boolean argument PACKUP must be set to false indicating that only a listing is required. The procedure will list all library items on the tape numbered STREAM which has its directory in block DIRECTORY (The directory block is normally the first block on the magnetic tape unless it has been repositioned by the user).
For example:
JOB
I0000 BLOGGS LIST LIBRARY ITEMS
DISC 99 R502 ALGOLIB
TAPE 1 N0000MYTAPE*INHIBIT
COMPILER ALGOL
INPUT INTERNAL ICT WITH ICT I/O PROCEDURES;
begin
library ICT18;
whats on my library tape(1,1,false,O,O);
end
***Z
The output would be something like:
NEXT AVAILABLE BLOCK IS 370
NAME NUMBER SIZE WHERE
30
ABC 1 173 100
5 90 3
27
DEF 19 20 350
GHI
GAPS ARE AS FOLLOWS
WHERE SIZE
93 7
273 77
YOUR PERCENTAGE OF GAPS IS 22
This indicates a library tape which has library names defined as ABC, DEF, GHI. The possible library names are ABC0 to ABC30, DEF0 to DEF27, GHI0 onwards. Of these only ABC1, ABC5, DEF19 have been defined. In the example, one large gap containing no current information occurs at block 273 for 77 blocks and a smaller gap at block 93 for 7 blocks. The unused blocks amount to 22% of the total space.
These gaps are causing much longer searches down the magnetic tape than are necessary and should probably be removed. The long gaps will cause the EXECUTION time associated with the job to increase and the job to run inefficiently.
This same procedure, when PACKUP is true, can be used to pack up the library tape by removing the gaps. In this case the new library tape is defined as being the tape numbered STREAMNEW and the directory will be put in block DIRECTORYNEW. It is recommended that a different tape from the old library tape is used for the new library tape. If the old library tape is used (STREAM = NEWSTREAM) then any computer failure could cause the total loss of all information from the old tape. The only way this could be safeguarded against while using the same tape would be to define the new directory position further down the tape than the next available block on the old library tape. If a packed up library tape is required then, in the Job Description, a COMMON tape must be defined as a temporary working space for the procedure and be called tape 90.
The library tape given in the first example could be packed up by following program:
JOB
I000 BLOGGS PACK UP LIBRARY TAPE
DISC 99 R502 ALGOLIB
TAPE 1 N0000MYTAPE*INHIBIT
TAPE 2 N0001MYOTHERTAPE*PERMIT
TAPE COMMON 90
COMPILER ALGOL
INPUT INTERNAL ICT WITH ICT I/O PROCEDURES;
begin
library ICT18;
whats on my library tape (1,1,true,1,2);
end
***Z
This defines the new packed up version of the library tape on tape 2 with directory in block 1.
All printing done by this procedure uses the ICT I/O procedures so that these are obligatory in the PROCESSOR command. Two serendipital properties of this procedure are:-
procedure erase(STREAM,DIRECTORY,NAME,NUMBER);
value STREAM,DIRECTORY,NUMBER;
integer STREAM, DIRECTORY, NUMBER;
string NAME;
The removal of library items from a library tape when they become obsolete will result in efficient use of library tapes and is to be encouraged. This is particularly important in the case of library tapes on the disc where space is at a premium.
The above procedure erase, will change the directory entry for the library item having name, NAME, and position, NUMBER, to free. The blocks, that the library item was using, will now be designated as free and will constitute a gap on the library tape (see ICT18). This gap could subsequently be removed from the library tape by using ICT18.
In the first example in the description of ICT18, the library item ABC5 could be deleted by calling:-
erase(1,1, 'ABC' ,5);
If erase was called before packing up the library tape (see second example of ICT18) then the area occupied by ABC5 would also be removed.
Possible diagnostics output by the procedure are:-
This library item contains the set of KDF9 procedures for outputting Algol arrays to magnetic tape and reading them back into store. Each array stored on magnetic tape is given a name by the user and this can be used to find the array when required. If all arrays stored on the tape have unique names then the array name can be presented to the read procedure and this will search the tape until a match is found with the name of an array on the magnetic tape. Alternatively if information is to be stored and retrieved in a serial fashion then the names can be ignored (all arrays can be given the same name) and the position of the tape will define which array is to be read into store. The procedures are as follows:-
procedure write binary(D, A, S);
value D;
integer D;
real array A;
string S;
All elements of array A are output to the magnetic tape numbered D in the Job Description and the array is given the name S on the magnetic tape. If this is the first use of the magnetic tape then it will be rewound before writing takes place. All information further down the magnetic tape than the last write binary order is assumed unobtainable. Thus it is not possible to overwrite information in the middle of the currently used region. All additions to the magnetic tape should therefore be at the end of the currently useful information on the magnetic tape.
procedure read binary(D, A, S);
value D;
integer D;
real array A;
string S;
This procedure will examine the next array past the current position of the magnetic tape D and check to see if its name coincides with S. If it does not agree, the magnetic tape is rewound and each array is examined in turn starting from the front of the magnetic tape until a match is found. Once a match is found, that array is read down into the array A. It is wise to keep the dimensions and limits of arrays the same in the read and write orders.
procedure rewind(D);
value D;
integer D;
The magnetic tape D is rewound to its starting position.
procedure skip(D,N);
value D, N;
integer D, N;
The magnetic tape D will be repositioned by moving the tape over the next N arrays. A negative value of N is allowed to move the tape backwards.
procedure data skip(D);
value D;
integer D;
The magnetic tape D is repositioned at the end of the currently available information on the tape. The tape would then be positioned ready for adding new information to the tape.
Boolean procedure BTC(D);
value D;
integer D;
BTC is true if tape D is at its initial position, otherwise false.
A more flexible format for the output of data can be obtained by using the KDF9 write procedure than any available in the standard ICT set of I/O procedures. This procedure has therefore been defined as ICT22.
procedure write kdf9 (F, Q);
value F,Q;
integer F;
real Q;
This procedure outputs the quantity Q in a form which depends on the format F. The value of F must be set to the result of calling the function procedure:-
integer procedure format(S);
string S;
The string S is set up as a pattern showing the form that the number should take on the output media.
<initial>spaces ::= <empty> | <unsigned integer> s| |s |ss |
sss | ssss | sssss | ........ | sssssssssssssss
<sign> ::= <empty>|-|+|≠
<exponent sign> ::= -|+|≠
<exponent layout> ::= 10<exponent sign> nd
<terminator> ::= <empty> | ;
<editing symbols> ::= <empty> | c | cc | ccc | p
<zeros> ::= 0 | <zeros> 0 | <zeros> s | s <zeros>
<positions> ::= d | <positions> d | <positions> s | s <positions>
<fractional layout> ::= .<positions>
<floating layout> ::= d <fractional layout>
<digit layout> ::= <positions> | <fractional layout>|
<positions> <zeros> | <positions> . <zeros> |
<positions> <fractional layout> <zeros> |
<positions> <fractional layout> |
<positions> <zeros>. <zeros>
<decimal layout> ::= <digit layout> | n <digit layout> | s <decimal layout>
<layout tail> ::= <decimal layout> <terminator> <editing symbol> |
<floating layout> <exponent layout>
<terminator> <editing symbol>
<layout> ::= <initial spaces> <sign> <layout tail>
<format identifier> ::= <identifier>
<format statement> ::= <format identifier> := format ( '<layout>')
A different integer identifier is used for each format required. This identifier is then used as a parameter to a write kdf9 statement.
Tfle symbols of the layout give a picture of the digits, spaces and symbols as they will appear on the printed sheet. The finally printed number will have exactly the same number of printed characters as is present in the layout (except in the case of alarm printing and the characters c, p and the special case if initial spaces). The various symbols of the layout have the following significance:
Alarm printing will occur when the number is too large for the layout or a negative number occurs with a layout having an empty sign position. Alarm printing will also occur if the format is syntactically incorrect. The number will be printed in standard floating format marked with an asterisk as follows:
*≠d.ddddsddddsddds ≠nddc
The alarm printed number will be on a new line followed by another new line with the appropriate number of spaces, so that the overall layout is not disturbed.
For example:-
i := format ('≠ndd.d');
write kdf9 (i, 37.28);
or alternatively
write kdf9 (format ('≠ndd.d'), 37.28);
would output
+37.3
ICT22 comprises the two procedures write kdf9 and format.
This consists of the two procedures:-
integer procedure char in kdf9;
procedure char out kdf9(C);
value C;
integer C;
The procedure char in kdf9 will read the next character of the selected data stream and the Atlas internal code value of this character will be the result of the procedure. For example reading the digit 1 gives the result 17 and A gives the result 33. The case shift causes 64 to be added to the internal code value. Reading a, for example, has the resulting value 97.
The procedure char out kdf9 will output to the selected output stream the character having internal code number C.
These two procedures are particularly useful for testing for the appearance of certain characters on the input and the output of unusual forms.
procedure copy text kdf9(S);
string S;
This procedure causes characters to be copied from the selected input to the selected output. If the string S consists of a single character then copying will take place from the current position of the input stream and will continue until the character in S is encountered. If the string S consists of two characters then all characters on the input preceding the appearance of the first character are ignored. Copying then takes place until the appearance of the second character in the string. The characters in S which indicate the start and finish of the text to be copied will not themselves be copied.
procedure read array kdf9(A, S);
real array A;
string S;
This procedure will read from the currently selected input an array A punched in the following form:
S; w; d; m1; m2; m3;........md; element; element; ........;element;
The symbols have the following significance:
If w = 1, the order in which the elements are expected from the tape is that given by first taking the lower bound values of the first (d-1) subscripts, and taking the dth subscript through its range from lower bound to upper bound; then stepping on the (d-1 )th subscript and taking the dth again through its range; and proceeding thus through the subscripts from right to left. In the case of a two-dimensional array regarded as a matrix, this is reading the elements by rows.
If w = -1, the subscripts are worked through from left to right. For a matrix this is reading the elements by columns. This is regarded as the usual form, rather than reading by rows.
d is the number of subscripts. (d = 2 for a matrix.)
mr is the range of the rth subscript.
procedure write array kdf9(F, D, A, S);
value F, D;
integer F, D;
real array A;
string S;
The array name S is output on the currently selected output followed by the initial parameters necessary for reading the array back in again using ICT25. Following this, the elements of the array are printed in columns.
The library item ICT27 consists of the four Elliott procedures:- address, lowbound, range and size. These four procedures allow the user access to all the information in the array's dope vector. They can be used to generate efficient array handling in code procedures.
integer procedure address(A);
array A;
The result of this procedure is the address of the first element of the array A.
integer procedure lowbound(A, I);
value I;
array A;
integer I;
The result is set to the lower subscript bound for the Ith subscript position from the left in the array A.
integer procedure range(A, I);
value I;
array A;
integer I;
The result is the number of array elements in the range of the Ith subscript position from the left in the array A.
integer procedure size(A);
array A;
The result is the number of elements In the array A.
For example the address of the array element A[I,J] is:-
(I _- lowbound (A,1)) * range (A,2) + J - lowbound (A,2) + address (A)
procedure advance elliott;
This procedure advances the currently selected input stream over the next character.
Boolean procedure buffer elliott(S);
string S;
The value of this Boolean function is true if the next character to be read on the currently selected input is the same as the single character defined as the string S; otherwise the result is false.
integer procedure dialect list (D, S, N, NU, I, T);
value D, S, NU, I, T;
integer D, S, NU, T;
Boolean I;
string N;
This procedure is designed to output a library item in any form or dialect that the user requires. It calls a procedure initialise dialect to be defined by the user which will give the representations of the Algol basic symbols on the particular output form. It is hoped eventually to provide a set of procedures initialise dialect for the main dialects used at the Laboratory.
Tne procedure dialect list outputs the library item (on the magnetic tape numbered S which has its directory at position D) with name N and number NU. The output will be on output stream T. The layout of the Algol library item will be the standard one used at the Laboratory. Indenting will be provided if I is true.
For example if ALGOLIB is defined as tape 99 in the Job Description then:-
dialect list (1, 99, 'ICT', 29, true, 3);
will output the library item ICT29 with indenting on output stream 3 (ALGOLIB has its directory at block 1, the standard position).
ICT30 and ICT31 define copies of the procedure initialise dialect which will produce output in the ATLAS CARDS dialect and ICT 8-HOLE TAPE respectively.
The procedure initialise dialect is defined as:
procedure initialise dialect (D, T);
integer arrayD, T;
The array T is used to store the output forms of the various Algol symbols and the array D is a set of pointers, one per Algol symbol, giving the position of each basic symbol in the array T. The dimensions of the arrays D and T are:-
D [ -9: *] and T[O: *]
The upper limits must be sufficient for the dialect in question. For the positive values of the array D, D[I] points to the position in T of the Ith Algol symbol. The numbering is the standard one used in the Atlas Algol System. The negative entries for D should contain the following:-
For each Algol symbol, the entries in the array T consist of a count of the characters in the symbol followed by the individual characters. If the output is in Internal Code then the internal code values of the characters are given. If the output is in binary, then the numeric values of the characters are output. On paper tape, the three least significant bits of the character are those together on one side of the sprocket feed hole. For example, 11.011 has the numeric value 27.
The Algol symbol Z would require the following:-
D[26] = n T[n ] = number of characters to be output for symbol Z. T[n+1] = first character of Z. T[n+2] = second character of Z. , etc.
In addition, the entries D[I] contain bits defining the layout required for the output. The value n above has added to it:
32768 if a space is required to be output before the Algol symbol 65536 if a space is required to be output after the Algol symbol 131072 if a newline is required to be output before the Algol symbol 262144 if a newline is required to be output after the Algol symbol
The standard settings used in the Atlas Algol listings can be obtained by calling
add in newline space (D);
before leaving the procedure initialise dialect.
This version of initialise dialect together with ICT29 allows library items to be output on cards in the ATLAS CARDS dialect.
This version of initialise dialect together with ICT29 allows library items to be output in the ICT 8-hole tape dialect used on 1900 computers. The output stream used must be defined as SEVEN-HOLE TAPE in the Job Description and a note added to the Job Slip saying that 8-track output is required.
procedure outjobtitletimedate(I);
value I;
integer I;
The procedure outputs on the currently selected stream:-
if I = 1 the Job Title if I = 2 the current time if I = 3 the date.
procedure chain (TAPE, MINBLK, MAXBLK, DEFINEONLY);
value TAPE, MINBLK, MAXBLK, DEFINEONLY;
integer TAPE, MINBLK, MAXBLK;
Boolean DEFINEONLY;
procedure callchain (TAPE, MINBLK, MAXBLK);
procedure returnchain (TAPE, MINBLK, MAXBLK);
The above procedures allow the user to split up a large Algol program into several sections (called chains) which are complete Algol programs. It is then possible to call one of these chains as a subroutine of another and later return to the calling chain. This facility would, of course, be of little value if there was no means of communicating data values between the two chains. This is achieved by insisting that the outer block declarations for all chains are identical as far as the number of variables of each type, the size of arrays and the ordering of the declarations are concerned. Variables may however be given different names in two chains. The equivalence of names depends on the position of the declaration in the outer block.
A program is defined as a chain by calling the procedure chain from the outer block level of the program. A convenient but not essential position would be the first executable statement of the program. If this is later called from another chain then it would be entered at the statement following the chain procedure call. The first three arguments to the procedure chain define a magnetic tape or disc area numbered TAPE where the program is dumped. It is dumped starting at block MINBLK and the program must not exceed the block MAXBLK. These arguments are very similar to those for the procedure dump program (ICT1). When defining a subchain, once the procedure chain has been called, it is often the case that the program may be terminated. The argument DEFINEONLY, if set to true, will terminate the program and re-enter the compiler. The user could then follow the sub-chain program by a second processor command and a second Algol program. This way all chains of a program can be defined in a single run. Alternatively if DEFINEONLY is false, the program continues executing once the chain has been defined. This would normally be used after defining the main program.
To enter one chain from another, the procedure call chain is called giving the position where the chain to be called is defined. The original chain can be returned to by calling return to chain. The arguments, in this case, give the position where the original chain is stored.
It is possible for sub-chains to call other sub-chains as long as a chain is not called recursively.
An example of a program defining two chains on tape 1 is as follows:
JOB
I0000 HOPGOOD TEST CHAIN CALLING
TAPE 1 N0000YOURTAPE~~PERMIT
.......
COMPILER ALGOL
INPUT INTERNAL ICT WITH ICT I/O PROCEDURES.
LIST ON LINEPRINTER;
begin
real a,b,c;
integer d,e,f;
chain (1 ,1,40,true);
a:=a+2;
b:=b+3;
d:=1;
return to chain(1,41,80);
end
INPUT INTERNAL ICT WITH ICT I/O PROCEDURES.
LIST ON LINEPRINTER;
begin
real ap,bp,cp;
integer d,e,f;
chain (1,41,80,false);
ap:=3;
bp:=2;
d:=4;
call chain (1,1,40);
print (ap,2,0); print (bp,2,0); print (d,2,0);
end
***Z
This program would print the results
5 5 1
It may be of interest to note that the I/O package used in each chain need not be the same. Complications could arise, however, if the Elliott I/O package was used in only one chain as the declaration:
array location [location : location];
is inserted at the top. It would be necessary to insert a dummy declaration:-
array location[1 : 1];
in the other chain not using the Elliott I/O package.
procedure call algol;
This procedure, when called, will cause the Algol compiler to be re-entered ready to compile another Algol program. The program will be expected on the currently selected input stream which should be in internal code mode and have a processor command on the next line to be read.
procedure dump program(TAPE, MINBLK, MAXBLK);
value TAPE, MINBLK;
integer TAPE, MINBLK, MAXBLK;
This procedure is similar to the standard dump program procedure, ICT1 but more efficient on reading down. Unlike the standard procedure, the directory block is stored at MINBLK rather than MAXBLK and it is followed by the blocks of the program. The variable MAXBLK contains the last block used on exit from dump program. The procedure uses the block *17 to read down the directory. The standard means 0 entering the dumped program is therefore:-
COMPILER ABL 1001, TAPE, 0, MINBLK 1002, TAPE, 0, J17 121, 127, 0, J17 EA
where TAPE is the number allocated in the Job Description for the tape containing the dumped program and MINBLK is the starting block for the dumped program.
The procedure is shorter than the standard procedure. It also remembers which I/O streams were selected when the program was dumped. It does not print anything on program dumping unless the program is longer than the space allocated. In this case the diagnostic
DUMP SPACE EXCEEDED
is issued.
This procedure was provided by J. Bailey. Further details are given in Chapter 15 of the Atlas Algol System manual.
procedure lower triangle iliffe vectors (A);;
array A;
This procedure enables a symmetric matrix of order N to be stored in a total of 0.5N (N + 1) words, while still treating the matrix as a two dimensional array. When referring to an element A[I,J], the condition J ≤ I must be satisfied, and the array A should be declared as A[1:N, 1:0.5 *(N + 1.1)]. The procedure modifies the Iliffe vector accessing the array elements and must be called before any elements the array are used.
This procedure was provided by J. Bailey.
procedure interpret(I);
value I;
integer I;
This procedure is designed for people with a good understanding of Atlas machine code. It switches a program from running normally into a state where the program is interpreted one instruction at a time and relevent information about changes in registers and storage locations is printed. The interpreter is stored on the disc area R072 INTERPRET. The disc area must be defined in the Job Description and the number allocated must be the same as the one given as the argument I to the procedure. Each instruction obeyed generates a line of print consisting of:-
The value of the modified address or immediate operand if applicable.
This is written as N = ........ for an immediate instruction and S' = ......... for an address.
Fuller details of the output are given in the write-up 'FOOTPRINTS' by W F Lunnon, Manchester University (Jan. 1967). The procedure interpret is a simplification of this very comprehensive system. The interpreter uses B51 and care must be taken not to destroy its value. In addition the interpreter is read down into blocks *0773 to *0776 and these cannot be used by the Algol program (no standard Algol program uses these blocks).
integer procedure array offset (A);
array A;
procedure to tape (I,A);
value I;
integer I;
array A;
procedure from tape (I,A);
value I;
integer I;
array A;
procedure move tape (I,J);
value I, J;
integer I, J;
These procedures allow the user to transfer information between the array A and magnetic tape or disc using fixed block 512 word transfer This is the most efficient method of storing and retrieving information from backing devices. As the array A will be stored in the program's stack of variables, it is unlikely that the base of the array A will be situated at a page boundary (all fixed block transfers must occur from and to 512 words starting at a page boundary).
Therefore to define an area of the array A coincident with a page of Atlas store, the user should define array A as
array A [1:1024];
then S: = array offset (A); will define the element of the array A which starts an Atlas page. The elements A[S] to A[S+511] will take part in the fixed block transfer.
The procedure totape will transfer the contents of the locations A[S + 1] to A[S + 512]to the next block on tape or disc given the number I in the Job Description. Similarly the procedure fromtape will read the next block of tape or disc I into the array positions A[S + 1] to A[S + 512].
The procedure move tape will reposition tape I to point to block J. The first block of tape is given the position 1 so that:-
movetape (I, 1)
is equivalent to the rewind operation.
procedure fast char print routine;
When called, this procedure replaces the standard output routine by one which is 2.5 times faster. All printing normally involves adding characters to a buffer. This allows composite characters to be output but does involve a certain amount of inefficiency. The fast routine cannot handle multiple characters or outer set characters but it does permit the user to freely intermix output from the standard procedures and output generated in code procedures by extra code. A further call of the procedure will result in a return to the original character printing routine. This procedure was provided by J. Bailey.
procedure chain (TAPE, MINBLK, MAXBLK, COMPILER);
value TAPE, MINBLK, MAXBLK, COMPILER;
integer TAPE, MINBLK, MAXBLK, COMPILER;
procedure call chain (TAPE, MINBLK, MAXBLK);
procedure return to chain (TAPE, MINBLK, MAXBLK);
procedure claim chain (TAPE, MINBLK, MAXBLK);
These procedures are derived from the original ICT38. They are slightly more restrictive than ICT38 but also more efficient.
The procedure chain defines an Algol program as a chain and stores it on tape TAPE between blocks MINBLK and MAXBLK. If the variable COMPILER is set:-
< 0 then the job terminates = 0 then the chain just defined is executed > 0 then calls compiler with number COMPILER
for example
2 will call down ABL 12 will call down ALGOL 19 will call down SERVICE 29 will call down SPECIAL
The original procedure return to chain would redefine the currenit program on tape before bringing down the new chain. The new procedure assumes that the program itself has not been modified (this could only be done by code procedures) so that the program itself need not be re-stored before entering the new chain. Similarly call chain assumes that the program has not been modified during execution and so only stores local variables and constants before entering the new chain.
The final procedure claim chain has been defined for use when you are initially not in a chain at all or alternatively never wish to return to this chain. The new chain is brought down and entered. No information about the original chain is stored.
These procedures have been thoroughly tested and are recommended in preference to ICT38. They do not have some of the limitations ICT38 which arise when using different I/O packages in different chains. These procedures were provided by J. Bailey.
procedure elegant output;
Boolean procedure near bottom of page;
These two procedures are designed to produce output which fits nicely onto computer output pages. It incorporates the fast print routine, ICT44, and has the same restrictions. After elegant output has been called, the output passes through a modified output routine which:-
Consequently each output page consists of 38 lines of output centred on the page. The procedure uses B66 as a line count and works on the assumption that output is being produced for one output stream at a time.
The procedure near bottom of page gives an answer true if more than 35 lines of output have been generated on the current page. These procedures were provided by J. Bailey.
This contains the procedures layout, inlogical and outlogical. Their definitions are identical to the procedures format, read boolean and write boolean in the KDF9 I/O package. They are required when running EGDON KDF9 dialect programs on Atlas. Full details are given in section (xii) of Chapter 4.3 of the Atlas Algol System manual.
real procedure elapsed time;
The result of this procedure is the time used in seconds since the program started running. This includes time spent in compilation.
ICT50 - 99 are a set of procedures for various numerical operations. The methods employed are among the most reliable discovered so far and the procedures should be fully debugged.
No procedure is ever completely foolproof, and a brief note is made indicating where difficulties may arise. On the whole, the procedures will be most satisfactory for medium-sized problems where considerations of storage and/or computing time are not of over-riding importance.
Items ICT55, 56, and 57 are from Algol 60 Procedures in Numerical Algebra, Part I by T J Dekker (Mathematisch Centrum Amsterdam Tract no. 22).
All other items are selected from Procedures Algol en Analyse Numerique published by the Centre National de la Recherche Scientifique, Paris.
procedure gresolpiv (A,B,X,N,IMPOSSIBLE);
value N;
array A,B,X;
integer N;
label IMPOSSIBLE;
Solution of AX = B by Gaussian elimination with row pivoting. N is the order of A. Exit to IMPOSSIBLE if A is singular. A and B are preserved. Declare A[1 :N, 1:N], B,X[1:N].
procedure tridiagreson (N,A1,B1,C1,F,IMPOS);
value N;
integer N;
real array A1,B1,C1,F;
label IMPOS;
Solves AX = B where A is tridiagonal, with diagonal elements in A1[1:N] and sub- and super-diagonal elements respectively in B1,C1[2:N]. Solution in F[1:N]. Exit to IMPOS if any minor is less than 10-6 in modulus. A1,B1, and C1 are destroyed.
procedure deter (A,N);
value N;
integer N;
real array A;
Determinant of A[1:N,1:N] by Gaussian elimination with row pivoting. A is preserved.
procedure trace (A,T,N,P);
value N,P;
integer N,P;
real array A,T;
Traces of powers of A [1:N,1:N] up to AP are produced in T[1:P]. A is preserved.
procedure pivotmax (A,B,N,IMPOSSIBLE);
value N;
array A,B;
integer N;
label IMPOSSIBLE;
The inverse of A[1:N,1:N], computed by Gaussian elimination with row interchanges, is stored in B[1:N,1:N]. Exit to IMPOSSIBLE if A is singular. A is preserved.
Cholesky decomposition of a positive definite symmetric matrix.
procedure detsym 1(A,N);
value N;
array A;
integer N;
Determinant of the +ve def. symmetric matrix M whose upper triangle is input in successive columns to A[1:N×(N+1)/2]. The Cholesky matrix of M is overwritten on A.
procedure solsym 1(A,N,B);
value N;
array A,B;
integer N;
Called after detsym 1 and solves Mx = B. The solution x is overwritten on B[1:N]. The Cholesky matrix of M is preserved in A so the procedure can be called repeatedly for various B.
procedure detinvsym 1(A,N);
value N;
array A;
integer N;
Determinant of M. Also the upper triangle of M-1 is overwritten on A (if det(M) ≥ 0). The procedure calls detsym 1. Note that in detsym 1 the decomposition is discontinued if M is not positive definite. In that case detsym 1:= -k, where k is the stage number at which a fault occurs.
procedure detbnd(A,N,LW,RW,AUX,M,P);
value N,LW,RW;
integer N,LW,RW;
integer array P;
array A,M,AUX;
Determinant of a band matrix order N with LW non-zero left codiagonal elements and RW right codiagonal elements. The elements of the matrix are input into A [1 : (N-1 )×(LW + RW) + N] successively by rows. A relative tolerance must be specified as AUX[0]. Method is standard LU decomposition. U is delivered in A and L in M[1:(N-2)×LW+1]. Pivotal indices are declared in P [1:N].
procedure solbnd(A,N,LW,RW,M,P,B);
value N,LW,RW;
integer N,LW,RW;
integer array P;
array A,B,M;
Solves the system Rx = B where R is the original band matrix described above. The solution vector x is overwritten on B. The procedure should be called after detbnd (if the determinant is non-zero). Since solbnd leaves A and M unaltered, it can be called successively for different right-hand sides B.
procedure lsqdec(A,N,M,AUX,AID,CI);
value N,M;
integer N,M;
array A,AUX,AID;
integer array CI;
Rank of the N × M matrix A [1:N,1:M]. AUX [0] is set in AUX [0:2].
procedure lsqsol(A,N,M,AID,CI,B);
value N,M;
integer N,M;
array A,AID,B;
integer array CI;
Least squares solution x of Ax = B is written on the first M elements of B. The Householder-triangularised form of A is left in A, except for the diagonal elements which are output in AID [1:M]. The pivotal indices are left in CI[1:M]. The procedure should only be called if Rank = M. A, AID and CI are left intact so the procedure can be repeatedly called for different RHS. B.
procedure lsqdglinv(A,M,AID,CI,DIAG);
value M;
integer M;
array A,AID,DIAG;
integer array CI;
Can be called after lsqdec if Rank = M, and calculates the diagonal of the inverse of A'A in DIAG [1:M].
procedure lsqdecsol(A,N,M,AUX,DIAG,B);
value N,M;
integer N,M;
array A,AUX,DIAG,B;
Combines the 3 preceding procedures to compute the least squares solution of Ax = B. The elements of A are altered.
procedure pseudoinv(B,M,N);
array B;
integer M,N;
Pseudoinverse (transposed) of a rectangular or non-singular square matrix B, based on the relation BI = (B'B)IB'. B is M × N where M ≤ N. Note that if rank(B) = N then BI = (B'B)-1B'. Declare B[1:M,1:N]. B is destroyed.
procedure pseudoinverse(A,M,N,IMPOSSIBLE);
array A;
integer M,N;
label IMPOSSIBLE;
Pseudoinverse (transposed) of a rectangular or (possibly singular) M × N matrix A. The procedure borders A to obtain a non-singular matrix hence we most declare A [1:M+N, 1:M+N]. An exit to IMPOSSIBLE is only possible if there b a machine failure, or the M+N - order matrix happens to be nearly singular. A is destroyed.
procedure eigenvalues(N,A,EV);
value N;
integer N;
real array A,EV;
Eigenvalues of a real symmetric matrix A[1:N,1:N] by Householder's method. They are arranged in descending order in EV[1:N]. Only the lower triangle of A need be input.
procedure iterjacobi(A,N,R,UT,PREC);
integer N;
real PREC;
real array A,R,UT;
Eigenvalues and normalised eigenvectors of a real symmetric matrix A [1:N,1:N] by Jacobi's method. The eigenvalues are in R[1:N] and the eigenvectors in UT[1:N,1:N] such that the jth row of UT is the eigenvector corresponding to R[j]. The iteration is terminated when
(Σi,j | A[i,j] |) / ( Σi A[i,i] ) ≤ PREC.
10-7 is a reasonable value for PREC. Only the upper triangle of A need be input. Remarks: Jacobi's method is considerably slower than Householder's as far as the eigenvalues are concerned but does have the advantage that an orthogonal set of eigenvectors is produced automatically in cases of multiple roots. Not recommended for N > 50.
procedure jacobi(N,AR,AI,LAMBDA,VR,VI,EPSILON);
integer N;
real EPSILON;
real array AR,AI,LAMBDA,VR,VI;
Eigenvalues and normalised eigenvectors of the complex Hermitean matrix AR + iAI by Jacobi's method. Only the upper triangle need be input. The eigenvector corresponding to the eigenvalue LAMBDA[j] is the jth column of VR + iVI. Declare AR,AI,VR,VI [1:N,1:N], LAMBDA [1:N]. The iteration is terminated when
( Σi<j (| AR[i,j] | + | AI[i,j] | ) ) / Σi=1 N | AR[i,i] |
is less than EPSILON ( = 10-7 say ). Not recommended for N > 20.
procedure rutisvalues(N,A,LAMBDAR,LAMBDAI,COUNT,FAULT,PREC);
integer N,COUNT;
array A,LAMBDAR,LAMBADAI;
label FAULT;
real PREC;
Real or complex eigenvalues of a general real N × N matrix A by Wilkinson quasi-triangulation followed by Rutishauser's method. Real and imaginary parts are stored in LAMBDAR,LAMBDAI[1:N]. The procedure exits to FAULT if a precision determined by PREC is not achieved af COUNT iterations. We suggest PREC = 10-10 and COUNT=50. The procedure is quite fast but multiple eigenvalues are found to less precision.
procedure hymanvalues(N,A,LAMBDAR,LAMBDAI,M,EPS);
integer N,M;
array A,LAMBDAR,LAMBADAI;
real EPS;
Real or complex eigenvalues of a general real N × N matrix A by Householder quasi-triangulation followed by Hyman's method. Real and imaginary parts in LAMBDAR,LAMBDAI[1 N]. The procedure finds M ≤ N eigenvalues, depending on what EPS is (it should find all of them as long as EPS is not set too small). EPS = 10-8 is small enough for practical purposes. The method is even more accurate than rutisvalues but 3-4 times slower for the same tolerance.
procedure eigenvec(N,A,LAMBDA,VEC,EPS,IMPOSSIBLE);
integer N;
real EPS;
array A,LAMBDA,VEC;
label IMPOSSIBLE;
Eigenvectors of a general real N × N matrix A with real eigenvalues. The procedure reduces A to tridiagonal form and calculates the eigenvectors by inverse iteration. If the reduction is singular an exit to IMPOSSIBLE is made. The eigenvalues are arranged in decreasing order in LAMBDA [1:N] and the eigenvectors are the corresponding columns of VEC[1:N,1:N]. We suggest EPS = 10-10.
There may be some instability in the reduction to tridiagonal form although this is not very common.
procedure newton complexe(A,B,N,EPS1,EPS2,EPS3,MAXITER,X,Y);
value N,A,B,EPS1,EPS2,EPS3;
integer N,MAXITER;
array A,B,X,Y;
real EPS1,EPS2,EPS3;
Finds the roots of (A[0] + iB[0]) × ZN+.. .+A[N] + iB[N] = 0 by Newton's method. The iteration is terminated when
max(| x(k+1) - x(k) | , |y(k+1) - y(k) |) < EPS1.
If no convergence is achieved after MAXITER iterations, EPS1 is multiplied by 10 and the process restarted (as long as EPS1≤EPS2). A root is considered real if |y(k)/x(k)| <EPS3. Declare A,B,X,Y [0:N]. The initial guess is X[0], Y[0] (which should be quite close to the origin) and the roots are stored in X[1], Y[1], ... , X[N], Y[N].
procedure muller(N,A,B,EPS,MAXITER,X,Y);
value N,A,B;
real EPS;
integer N,MAXITER;
array A,B,X,Y;
Finds the roots of (A[0] + iB[0]) × ZN + ....A[N] + ib[N] = 0 by Muller's method. The iteration terminates when
| Z(k+1) - z(k) | < EPS × Z(k+1) |, or after MAXITER steps. The initial guess is set internally. The roots are stored in X,Y [1:N].
procedure laguerre(N,A,MAXITER,EPS,EPS2,X,IMPOSS);
value N,A,EPS,EPS2;
integer N,MAXITER;
real EPS,EPS2;
array A,X;
label IMPOSS;
Finds the roots of A[0] × xN + ..... + A[N] = 0 by Laguerre's method, assuming they are all real. The iteration is terminated when | x(k+1) - x(k) | < EPS.
If convergence is not achieved after MAXITER iterations, EPS is multiplied by 10 and the process restarted, as long as EPS≤EPS2. We suggest MAXITER = 30 and EPS2 = 0.01. The initial guess is placed in X[0] and the roots are stored in X[1], ..... , X[N].
procedure bissection(F,Y1,Y2,EPS,ETA,X,FAIL);
value Y1,Y2;
real procedure F;
real Y1,Y2,EPS,ETA,X;
label FAIL;
Finds a zero X of F(y) in [Y1,Y2] by the method of bisection. An exit to FAIL is made if an even number of zeros is found. If there are an odd number then only one is found. Calculation terminates when F (X) < EPS or the interval is less than ETA in width.
procedure bisdou(F,G,A1,B1,A2,B2,EPS,X,Y,IMPOSS);
value A1,B1,A2,B2,EPS;
real procedure F,G;
real A1,B1,A2,B2,EPS,X,Y;
label IMPOSS;
Finds a solution of F(X,Y) = G(X,Y) = 0 in the rectangle [A1,B1] × [A2,B2] by the method of bissection. The solution will lie in the square [X,X + EPS] × [Y,Y + EPS]. Exit to IMPOSS if no solution is found.
These two procedures are rather slow and are most useful when used to supply a starting value for an iterative process such as Newton's method.
procedure rungekutta(X0,Y0,H,M,N,Q,TETA,A,X,V,F);
value X0,H,M,N;
real X0,H;
integer M,N,Q;
array Y0,TETA,A,X,V;
procedure F;
procedure integro 1(X0,Y0,H,M,N,X,V);
value X0,H,M,N;
real X0,H;
integer M,N;
array Y0,X,V;
Integrates the first order system y = f(x,y) from x = X0 to x = X0 + M × H where y and f are N-vectors. The method used is 4th-order Runge-Kutta. Values of x are stored in X[0:M] and the value of y at x = X0 + k × H is stored in the kth column of V[1:N,0:M]. rungekutta is called internally by integro 1, so that the only input parameters required are F,X0,Y0,H,M,N. Only integro 1 is called by the user. (The parameters Q,TETA, and A are set in integro 1).
Difficulties may arise if transient solutions are present since these will not be computed accurately unless H is very small. For example, the solution of y' = -λy (λ > 0) is not computed accurately if λH > 2.8.
The procedure F is defined as:-
procedure F(X,Y,G,N);
real X;
integer N;
array Y,G;
The procedure must define G[i] equal to fi(X,Y[1], ... , Y[N]);
procedure rungekutta(X0,Y0,YP0,H,M,N,Q,TETA,A,B,X,V,VP,F);
value X0,H,M,N;
real X0,H;
integer M,N,Q;
array Y0,YP0,TETA,A,B,X,V,VP;
procedure F;
procedure integro 2(X0,Y0,YP0,H,M,N,X,V,VP);
value X0,H,M,N;
real X0,H;
integer M,N;
array Y0,YP0,X,V,VP;
Integrates the 2nd order system y" = f(x,y,y'), where y, f are N-vectors, from x = X0 to x = X0 + M × H by a 4th order Runge-Kutta method. The initial conditions y(X0) and y'(X0) are stored in Y0,YP0[1:N] and the solution y,y' at x = X0 + k × H = X[k] is stored in the kth columns of V,VP[1:N,0:M]. Only integro 2 is called and the input parameters required are F,X0,Y0,YP0,H,M,N.
Tne procedure F is defined as:-
procedure F(X,Y,YP,G,N);
real X;
integer N;
array Y,YP,G;
The procedure calculates G[i] = fi(X,Y[1], ... Y[N],YP[1], ... YP[N])
procedure rungekutta(X0,Y0,YP0,H,M,N,Q,TETA,A,B,X,V,VP,F);
value X0,H,M,N;
real X0,H;
integer M,N,Q;
array Y0,YP0,TETA,A,B,X,V,VP;
procedure F;
procedure integrobis(X0,Y0,YP0,H,M,N,X,V,VP);
value X0,H,M,N;
real X0,H;
integer M,N;
array Y0,YP0,X,V,VP;
Integrates the special 2nd order system y" = f(x,y) (y, f are N-vectors) from x = X0 to x = X0 + M × H by a 4th order Runge-Kutta method. The significance of the parameters is the same as above i.e. initial values in Y0,YP0[1:N] and solution as columns of V,VP[1:N,0:M] Only integrobis is called by the user and input parameters required are F,X0,Y0,YP0,H,M,N. The definitLon of F is the same as in ICT76.
procedure insiro(F,A,B,ORDMAX,PREC,SORT,RES);
value A,B,ORDMAX,PREC;
real procedure F;
Boolean SORT;
real A,B,PREC;
integer ORDMAX;
real array RES;
Integral of F(x) from A to B by Romberg's rule halving the step at each stage. If the relative difference between successive estimates is less than PREC, an exit is made with SORT=true, while if the number of points increases beyond ORDMAX, SORT=false on exit. The array RES [1:ORDMAX + 1] contains successive approximations, the most precise being in RES[1]. Tne maximum number of function evaluations is 2ORDMAX.
procedure indourec(F,A1,B1,A2,B2,N1,N2,ORDMAX,PREC,SORT,RES);
value A1,B1,A2,B2,N1,N2,ORDMAX,PREC;
real procedure F;
Boolean SORT;
real A1,B1,A2,B2,PREC;
integer N1,N2,ORDMAX;
real array RES;
Integral of F(x,y) over [A1,B1] × [A2,B2] by Romberg's rule, using N1 and N2 steps initially. The other parameters have the same significance as for insiro. The maximum number of function evaluations is N1 × N2 × 4(ORDMAX + 1).
procedure int3neville (F,A1,B1,A2,B2,A3,B3,H1,H2,H3,ALPHA,ORDMAX,
PREC,INT,ENTlER,RES);
value A1,B1,A2,B2,A3,B3,H1,H2,H3,ALPHA,ORDMAX,PREC;
real procedure F;
real A1,B1,A2,B2,A3,B3,H1,H2,H3,PREC,ALPHA;
integer ORDMAX,ENTIER;
real array RES,INT;
Integral of F(x,y,z) over [A1,B1] × [A2,B2] × [A3,B3] by Romberg's rule , with the steps hi reduced from the initial values (H1,H2,H3) according to the rule ki = hi/n where n=entier(ALPHA × n)+1, n(1)=1. We must set ALPHA > 1 (say 1.5). The trapezoidal estimates are in RES and the Neville extrapolations in INT [1:ORDMAX]. If the relative difference between extrapolations after ENTlER iterations is less than PREC, the value of the integral is stored in INT[ENTlER]. The number of function evaluations is then 1 + 8 × (ALPHA (3 × ENTlER -1) / (ALPHA3 -1))/ (H1 × H2 × H3).
procedure intcossin (F,A,B,LAM,RESCOS,RESSIN,ERCOS,ERSIN);
value A,B,LAM,N;
integer N;
real procedure F;
real A,B,LAM,RESCOS,RESSIN,ERCOS,ERSIN;
Integrals of F(x) cos(LAM.x) and F(x) sin(LAM.x) over [A,B] by piecewise polynomial interpolation of F(x) at N points, with an error estimate based on the interpolating polynomial. Integrals in RESCOS ,RESSIN; error estimates in ERCOS,ERSIN.
Not recommended for LAM < 0.05.
procedure approxicon (F,VAR,GAM,R,M,N,A,B,NB,EPS,X);
value M,N,NB,EPS;
real procedure F;
procedure VAR;
integer M,N,NB;
real EPS;
array GAM,R,A,B,X;
Calculation of the best uniform approximation to a continuous function F(x) on the union of NB closed disjoint intervals [A[i], B[i]] by a linear combination x[1] × F1(x)+ ....... + X[N] × FN(x) satisfying the M linear constraints:
X[1] × GAM [1,1] + ............... + X[N] × GAM[1,N] = R[1] ............ ........................... ...... X[1] × GAM [M,1] + ............... + X[N] × GAM[M,N] = R[M]
The basic functions F1, ... , FN must be specified by
procedure VAR (x,y);
real x;
array y;
which associates Y[i] with Fi(x). EPS is usually set to 0.001. The method used is the Remes exchange algorithm, which calculates a lower bound mr and an upper bound Mr to the error norm. These are automatically output by the procedure. The process stops at the rth iteration when Mr - mr < EPS × Mr.
As a result of the generality of this procedure, the execution time may be quite long (2½ minutes for a simple example).
procedure remez (FCTF,N,PO,PI,X,EPS);
value N;
real procedure FCTF;
integer N;
real PO,PI,EPS;
real array X;
Best uniform approximation or a continuous function FCTF on [PO,PI] by a polynomial X[1] × YN + .... + X[N+1]. The error norm is stored in X[N+2]. Set EPS between 0.0001 and 0.1. Uses the Remes algorithm.
procedure tchebecha (M,N,Y,B,X);
array Y,B,X;
integer M,N;
Best uniform approximation of a continuous function on a discrete point set by a linear combination X[1] × F1 + .... + X[N] × FN of continuous functions. Y[1:M] must contain the values of the function and B[1:M,1:N] the values of the basis functions. The error norm is stored in X[N+1]. The method (Stiefel-Remes) converges in a finite number of steps.
procedure chebfit (M,N,X,Y,A);
integer M,N;
real array X,Y,A;
Best uniform approximation of a continuous function Y on M points X by the polynomial A[0] + A[1] × x + ..... a[N] × xN. The method is that of Stiefel-Remes.
procedure chebdessous (A,B,N,F,EPS,H,ALPHA,P,E);
value A,B,N,EPS,H;
real A,B,EPS,H;
integer N;
real procedure F;
real array ALPHA,P,E;
Best one-sided uniform approximation to a continuous function F on [A,B] by a linear combination ALPHA[1] × G1 + ...... + ALPHA[N] ×GN such that the error F-ΣI ALPHA[I] × GI is always positive. The values GI(x) must be specified externally by
procedure cal (X,G);
value X;
real X;
array G;
which associates GI with G[I]. The maximum error is output in ALPHA[N+1]. The values of the errors at the critical points P[1:N+1] are contained in E[1:N+1]. H is a parameter of the discretisation (say 0.001) and EPS is put between 0.1 and 0.001.
procedure tchfbing (M,N,K,G,F,C,A,E,EPS,BOOL);
value M,N,K,EPS;
Boolean BOOL;
integer M,N,K;
real EPS;
array G,F,C,A,E;
Discrete Chebyshev approximation with inequality constraints. The vector F[K+1,K+N] is approximated at N points by a linear combination H = A[1 ] × G1 + ... + A[M] × GM where the Gi have K + N components and from the rows of G[1:M,1:K+N]. The first K components of the approximating function H satisfy the K constraints H[j] ≤ C[j], for C[1:K]. The error norm is output in A [M+ 1] and the errors F - G in E [K+ 1:K+N]. Also the constraints are checked by outputting C - G in E [1:K]. Thus E is declared [1:K+N]. EPS = 0.01 or 0.001. If EPS is too small we have BOOL=false.
procedure mcpolysp (K,M,X,Y,SIGMA,S,ALPHA,BETA,ERROR);
integer K,M;
array X,Y,SIGMA,S,ALPHA,BETA,ERROR;
real procedure P(K,S,ALPHA,BETA,X);
integer K;
real X;
array S,ALPHA,BETA;
Least squares approximation of function values Y[i] at M points X [i] by orthogonal polynomials S[0] P0(x) + ...... + S[K] PK(x). The standard deviations of the random variables Y[i] are input to SIGMA[i]. The standard deviations of the coefficients S[i] are in ERROR[i]. The orthogonal polynomials are generated by standard recurrence relations with coefficients ALPHA[1:K] and BETA[0:K-1] . These are generated internally. The values of the approximating polynomial are calculated by P.
procedure methdintepoly (K,N,X,Y,D);
value K,N,X,Y;
integer K,N;
array X,Y,D;
real procedure polinte(K,N,X,D,T);
value K,N,X,D,T;
integer K,N;
array X,D;
real T;
Interpolating spline of order K which takes on the values Y[i] at points X[i], i =1, ..., N. The derivatives of order 0 to K-1 at these points are output in D[1:N,0:K-1]. Method works best for K ≤ 5, N ≤ 60. Values at any point x = T are produced by polinte.
procedure coefsplinetrois (N,X,Y,M);
value N,X,Y;
integer N;
array X,Y,M;
Cubic spline with values Y[i] at X[i], i = 0, ..., N. The second derivatives are output in M[0:N].
real procedure sp(N,X,Y,M,T);
value N,X,M;
integer N;
array X,Y,M;
real T;
Values of the cubic spline obtained above at the point x=T.
procedure splinehermite (N,Y,YP,X,L,M,CC);
value N,Y,YP,X;
integer N;
array Y,YP,X,L,M,CC;
real procedure hermspl(N,X,L,M,CC,T);
value N,X,L,M,CC;
integer N;
array X,L,M,CC;
real T;
The first procedure calculates the coefficients L[1:N-2], M[1:N-1], CC[1:3] of the Hermite spline taking values Y[i] and derivatives YP[i] at (monotonically ordered) points x[i], i = 1, ..., N. (N ≥ 4). The second procedure calculates the value at x=T of this spline.
real procedure splderivee(N,X,Y,T);
value N,X,Y,T;
integer N;
array X,Y;
real T;
Best approximation (in the sense of Sard) to the derivative at x=T of a function taking the values Y[i] at points X[i], i = 0,..., N. Method uses a cubic spline.
procedure classemarkoff(N,P,R,T,F);
integer N;
array P,R,T,F;
Calculation of non-recurring, transient, and persistent states of a Markoff chain from the boolean matrix P associated with the matrix of transition probabilities. N is the number of states, the non-recurring states are in R[1:N] (state i non-recurring if R[i] = 1), T[1:N] contains the transitory states (state i in transitory class k if T[i] = k), F[1:N] contains the final states (i in class k if F[i] = k). The arrays are declared real.
procedure sousclasse cyclique(N,P,L);
integer N;
array P,L;
Computation of cyclic sub-classes of a set of persistent states of a Markoff chain from the boolean matrix P associated with the stochastic matrix of order N. The results are in the rows of L[1:M,1:M] such that the ith subclass is represented by these elements of L such that L[i,j] = 1.
procedure polyweyl(A,X,R,K,M);
integer array A,X;
integer R,M,K;
Random vector X with R components independent, uniformly distributed in [0,1], in integer form with K digits in each component. We have to declare A,X[1:M+R] where M is the smallest integer satisfying M ≥ 109/K. The initial vector A can be formed from the first (M+R)K digits of π.
procedure notlag(T);
real T;
Inverse of the error function:
ICT100 to ICT105 are taken from Numerische Mathematik. These procedures were provided by Dr Natha.
procedure reduction(n,a,b,dl,fail);
value n;
integer n;
array a,b,dl;
label fail;
Reduction of the general symmetric eigenvalue with the symmetric matrix Ax = λBx, with the symmetric matrix A and symmetric positive definite matrix B, to the equivalent standard problem Pz = λz.
Ref: Numerische Mathematik 11,99-110 (1968), procedure reduce 1.
procedure rator(n,m,posdef,dlam,eps,d,b2);
value n,m,posdef,dlam,eps;
integer n,m;
Boolean posdef;
real dlam,eps;
array d,b2;
QR algorithm for the computation of the lowest eigenvalues of a symmetric tridiagonal matrix.
Ref: Numerische Mathematik 11, 264-272 (1968).
procedure tred(n,tol,a,d,e,e2);
value n,tol;
integer n;
real tol;
array a,d,e,e2;
This procedure reduces the given lower triangle of a symmetric matrix, to tridiagonal form using Householder's reduction.
Ref: Numerische Mathematik 11,181-195 (1968), procedure tred 1
procedure tridiinverseiteration(c,b,n,w,norm,m1,macheps,z);
value n,m1,norm,macheps;
integer n,m1;
real norm,macheps;
array c,b,w,z;
This procedure calculates the eigenvectors of a symmetric tridiagonal matrix by inverse iteration, given very accurate eigenvalues.
Ref: Numerische Mathematik 4, 368-376 (1962).
procedure backtransformation(a,b,x,n,m1);
value n,m1;
integer n,m1;
array a,b,z;
This procedure calculates the eigenvectors of the original symmetric matrix A, from those eigenvectors of the tridiagonal matrix which were computed by ICT103.
Ref: Numerische Mathematik 4, 354-361 (1962).
procedure bisect(c,beta,n,m1,m2,eps1,relfeh,eps2,z,x);
value n,m1,m2,eps1,relfeh;
real eps1,eps2,relfeh;
integer n,m1,m2;
array c,b,x,beta;
This procedure calculates the eigenvalues of a symmetric tridiagonal matrix by the method of bisection.
Ref: Numerische Mathematik 9, 386-393 (1967).
These procedures on character manipulation were provided by Mr Farmer, Birkbeck College.
integer procedure readch;
reads the next character on selected input and produces internal code value as integer.
procedure printch(i);
value i;
integer i;
outputs the character having internal code position i.
integer procedure peep;
looks at the next character on selected input and produces internal code value as integer, but does not alter input pointer.
procedure skip;
moves over next character on selected input.
This set of procedures enables the user to treat the lineprinter as a simple graph plotting device. The plotting area or field can be of any proportions. If the width of the field is greater than 120 Characters it is automatically broken into sections printed below one another. These procedures were provided by Mr Farmer.
procedure startgraph(leftx,topy,rightx,bottomy,chars,lines);
value leftx,topy,rightx,bottomy,chars,lines;
real leftx,topy,rightx,bottomy;
integer chars,lines;
A call to this procedure must precede all other plotting calls. This procedure sets up the list structure necessary for storing all the plotting information and defines the scales that the programmer wishes to use by specifying the co-ordinates at the top left hand corner and the bottom right hand corner of the field:
The line/space ratio of the lineprinter is 6/10 and should be taken into account when choosing coordinate scales.
procedure plot(x,y,char);
value x,y;
real x,y;
string char;
This plots the first character in the Algol string char at the position corresponding to (x,y). Compound characters are allowed and a dot is printed if the string is empty. If the point (x,y) lies outside the specified field, nothing is plotted.
procedure caption(x,y,string,horiz);
value x,y,horiz;
real x,y;
string string;
Boolean horiz;
The given Algol string is plotted starting from point (x,y). If horiz is true it is plotted horizontally, otherwise vertically. Any part of the string lying outside the field will not appear.
procedure border(char);
string char;
The first character of the Algol string char is used to border the field. This is useful if a field ends in the middle of a sheet.
procedure finishgraph;
This procedure must be called to produce any graphical output. All the information stored from the last call to startgraph will be output. The list structure is destroyed during output so the procedure should not be called until after a corresponding call to startgraph.
Any amount of information can be plotted at the same point and everything will be overprinted. It is not necessary for leftx < rightx and topy > bottomy although this is the standard orientation for graphs.
real procedure random;
This procedure generates a pseudo-random number in the range (0,1). A different sequence can be generated by changing the octal numbers at (1) and (1) + 1 within the procedure itself.
procedure chisquare(r,n);
value r,n;
integer r,n;
This procedure performs the chi-square test for randomness of a global procedure called random producing a random sequence in [0,1). It tests n observations dividing the range [0,1) into r intervals.
integer procedure timer;
A call to this procedure sets timer equal to the total number of instruction interrupts used since the start of compilation.