with(linalg); lllgram:=proc(A,n) local G,M,T,k,r,i,j,h; G:=array(1..n); # Phi(b_i^*,b_i^*) T:=evalm(array(identity,1..n,1..n)); #Transformationsmatrix M:=evalm(array(identity,1..n,1..n)); #die mu_{ij} k:=2; G[1]:=A[1,1]; while (k <= n) do print(k); if (k = 1) then k:=k+1; G[1]:=A[1,1]; else for j from 1 to k-1 do M[k,j]:=scapr(M,A,k,j,G); od; r:=round(M[k,k-1]); for j from 1 to n do A[k,j]:=A[k,j]-r*A[k-1,j]; T[k,j]:=T[k,j]-r*T[k-1,j]; od; for j from 1 to n do A[j,k]:=A[j,k]-r*A[j,k-1]; od; M[k,k-1]:=M[k,k-1]-r; for j from 1 to k-2 do M[k,j]:=M[k,j]-r*M[k-1,j]; od; G[k]:=A[k,k]; for j from 1 to k-2 do G[k]:= G[k] - M[k,j]^2*G[j]; od; #Laenge der Projektion von b_k auf ^{\perp } if (G[k] < 3/4 * G[k-1]) then for j from 1 to n do h:=T[k-1,j]; T[k-1,j]:=T[k,j]; T[k,j]:=h; h:=A[k-1,j]; A[k-1,j]:=A[k,j]; A[k,j]:=h; od; for j from 1 to n do h:=A[j,k-1]; A[j,k-1]:=A[j,k]; A[j,k]:=h; od; G[k-1]:=G[k]; k:=k-1; else G[k]:=G[k]-M[k,k-1]^2*G[k-1]; #Phi(b_k^*,b_k^*) for i from k-2 by -1 to 1 do r:=round(M[k,i]); for j from 1 to n do A[k,j]:=A[k,j]-r*A[i,j]; T[k,j]:=T[k,j]-r*T[i,j]; od; for j from 1 to n do A[j,k]:=A[j,k]-r*A[j,i]; od; M[k,i]:=M[k,i]-r; for j from 1 to i-1 do M[k,j]:=M[k,j]-r*M[i,j]; od; od; k:=k+1; fi; fi; #k>=2 od; RETURN(A,T,M,G); end; scapr:=proc(M,A,k,j,G) local m,i; m:=A[k,j]; for i from 1 to j-1 do m:= m - G[i]*M[k,i]*M[j,i]; od; m:=m/G[j]; RETURN(m); end;