Lehman factoring Algol implementation from 1974 research paper

Lehman_factor.alg

begin
   comment  v0.1, only factors sometimes as of now

   comment  runs with this compiler: https://github.com/JvanKatwijk/algol-60-compiler
      $ jff-algol Lehman_factor.alg 
      $ ./Lehman_factor
      32639 3 
      127 
      $
   ;

   comment  Lehman factoring algorithm, Algol code from:
   https://www.ams.org/journals/mcom/1974-28-126/S0025-5718-1974-0340163-2/S0025-5718-1974-0340163-2.pdf#page=8;

   integer procedure Mod (a, b); value a, b; integer a, b; Mod := a - (a / b) * b;

   integer procedure gcd(a, b); value a, b; integer a, b;
   begin integer i;
      if a < b then begin i := a; a := b; b := i end;
   l: i := Mod(a, b); a := b; b := i;
      if i != 0 then goto l; gcd := a
   end gcd;

   integer procedure isqrt(n, u); value n; integer n, u;
   begin integer j,j1,j2;
      j := if n = 0 then 1 else entier (sqrt(n))+1;
      j1 :=j*j - n;
   f: if j1 < 0 then
      begin j1 := j1+2*j+1; j:=j+1; goto f end;
   l: j2:= j1-2*j+1;
      if j2 >= 0 then
      begin j1 := j2; j:= j-1; goto l end;
      isqrt := j; u:= j1
   end isqrt;

   procedure factor(n, r, f); value n, r; integer n, r, f;
   begin integer i, j, p;
      integer array c[1 : 8];
      boolean array qr[0 : 728];
      procedure large(m, mO); value m, mO; integer m, mO;
      begin integer i, i1, j, jump, k, s, t, u, x, y; boolean odd;
         s := 1; k := mO;
      start:
         k := k+c[s]; s := if s=m then 1 else s+1;
         if k<=r then
         begin
            x:= isqrt(4*k*n, u); j := (isqrt(n/k,t) - 1) / (4*(r+1));
            if Mod(x+k, 2) = 0 then
            begin i1 := 1; u := u+2*x+1; x := x+1 end else i1 :=0;
            odd := Mod(k, 2) = 1; jump := if odd then 4 else 2;
            if odd then
            begin
               if Mod(k+n, 4) = Mod(x, 4) then
               begin i1 := i1+2; u := u+4*(x+1); x := x+2 end
            end;
            for i := i1 step jump until j+1 do
            begin
               if qr[Mod (u, 729)] then
               begin
                  y := isqrt(u, t);
                  if t = 0 then
                  begin
                     p := gcd(n, x-y); if p > n/p then p := n/p;
                     goto exit
                  end;
                  comment When a factor p is found, we leave the
                  procedure by going to exit;
               end;
               if odd then begin u := u+8*(x+2); x := x+4 end
               else
               begin u := u+4*(x+1); x := x + 2 end
            end;
            goto start
         end
      end large;
      for i := 0 step 1 until 728 do qr[i] := false;
      for i := 0 step 1 until 364 do
      begin j := Mod(i*i, 729); qr[j] := true end;
      c[1]:= 30;large(1,0);
      c[1]:= 48;c[2]:= c[3]:= c[4]:= 24;large(4, -24);
      c[1]:= c[2]:= c[4]:= 24;c[3]:= 48;large(4,-12);
      c[1]:= c[2]:= c[4]:= 36;c[3]:= 72;large(4, -18);
      c[1]:= c[4] := c[6] := 12; c[2] := c[8] := 36;
      c[3]:= c[5] := c[7] := 24; large(8, -6);
      c[1]:= 4;c[2]:= 2;large(2, -2);
      c[1]:= 2;large(1,-1);
      comment No factor has been found;
      p := 1;
   exit: f := p
   end factor;

   comment  transformed https://rosettacode.org/wiki/Nth_root#ALGOL_W;
   real procedure nthRoot(A, n, precision);
   value A, n, precision; real A; integer n; real precision;
   begin
      real xk, xd;
      integer n1, dummy;
      n1 := n - 1;
      xk := A / n;
      xd := precision + 0.1;
      for dummy := 0 while abs( xd ) > precision do
      begin
         xd := ( ( A / ( xk ** n1 ) ) - xk ) / n;
         xk := xk + xd;
      end;
      nthRoot := xk
   end nthRoot ;


   integer n, r, f;
   n := 257 * 127;
   r := entier(0.1 * nthRoot(n, 3, 1e-5));

   outinteger(1, n); outinteger(1, r); newline (1);
   factor(n, r, f);
   outinteger(1, f); newline (1);
end;