begin

library A0, A1, A5, A15;

comment GAMM in FORTRAN - single precision version, mark 3;
comment National Physical Laboratory Benchmark Gamm F;
comment This program has a single parameter N;
comment Set N := 10000 for about 60s on a machine like System 360/65;
comment Output is by one write statement to stream 30;
comment It takes KDF9 about 480s with unoptimised KAlgol,
                   and about 154s when optimised;

integer FIVE, I, J, N, REP, TEN, THIRTY;
real ACC, ACC1, DIVN, ROOT, X, Y;
real array A[1:30], B[1:30], C[1:30];

N := 10000;
FIVE := 5;
TEN := 10;
THIRTY := 30;
DIVN := 1.0 / N;
X := 0.1;
ACC := 0.0;

comment INITIALISE A AND B;
Y := 1.0;
for I := 1 step 1 until 30 do
begin
  A[I] := I;
  B[I] := - Y;
  Y := - Y;
end 

comment ONE PASS OF THIS LOOP CORRESPONDS TO 300 GAMM UNITS;
for REP := 1 step 1 until N do
     begin
comment FIRST ADDITION/SUBTRACTION LOOP;
   I := 30;
   for J := 1 step 1 until THIRTY do
     begin
     C[I] := A[I] + B[I];
     I := I - 1;
     end 
comment FIRST POLYNOMIAL LOOP;
   Y := 0.0;
   for I := 1 step 1 until TEN do
     begin
     Y := (Y + C[I]) * X;
     end 
   ACC1 := Y * DIVN; comment not I, which is undefined in Algol;
comment FIRST MAXIMUM ELEMENT LOOP;
   Y := C[11];
   for I := 12 step 1 until 20 do
     begin
     if (C[I] > Y) then Y := C[I];
     end 
comment FIRST SOUARE ROOT LOOP;
   ROOT := 1.0;
   for I := 1 step 1 until 5 do
     begin
     ROOT := 0.5 * (ROOT + Y/ROOT);
     end 
   ACC1 := ACC1 + ROOT * DIVN;
comment SECOND ADDITION/SUBTRACTION LOOP;
   for I := 1 step 1 until THIRTY do
     begin
     A[I] := C[I] - B[I];
     end 
comment SECOND POLYNOMIAL LOOP;
   Y := 0.0;
   for I := 1 step 1 until TEN do
     begin
     Y := (Y + A[I]) * X;
     end 
comment SECOND SQUARE ROOT LOOP;
   ROOT := 1.0;
   for I := 1 step 1 until FIVE do
     begin
     ROOT := 0.5 * (ROOT + Y/ROOT);
     end 
   ACC1 := ACC1 + ROOT * DIVN;
comment FIRST MULTIPLICATION LOOP;
   for I := 1 step 1 until THIRTY do
     begin
     C[I] := C[I] * B[I];
     end 
comment SECOND MAXIMUM ELEMENT LOOP;
   Y := C[20];
   for I := 21 step 1 until THIRTY do
     begin
     if (C[I] > Y) then Y := C[I];
     end 
comment THIRD SQUARE ROOT LOOP;
   ROOT := 1.0;
   for I := 1 step 1 until 5 do
     begin
     ROOT := 0.5 * (ROOT + Y/ROOT);
     end 
   ACC1 := ACC1 + ROOT * DIVN;
comment THIRD POLYNOMIAL LOOP;
   Y := 0.0;
   for I := 1 step 1 until TEN do
     begin
     Y := (Y + C[I]) * X;
     end 
   ACC1 := ACC1 + Y * DIVN;
comment THIRD MAXIMUM ELEMENT LOOP;
   Y := C[1];
   for I := 2 step 1 until TEN do
     begin
     if (C[I] > Y) then Y := C[I];
     end 
comment FOURTH SQUARE ROOT LOOP;
   ROOT := 1.0;
   for I := 1 step 1 until FIVE do
     begin
     ROOT := 0.5 * (ROOT + Y/ROOT);
     end 

     ACC1 := ACC1 + ROOT * DIVN;
     ACC := ACC + ACC1;
end 
open(30);
output(30, N);
output(30, ACC1*N);
comment Should print N then  1.6733 4322 4109 0064 7168 4801E1;
comment KDF9 actually prints 1.6733 4322 411E1;
end
|