module linmodp; % Routines for solving linear equations mod p.
% Authors: A. C. Norman and P. M. A. Moore, 1979.
fluid '(prime!-base);
symbolic procedure lu!-factorize!-mod!-p(a,n);
% A is a matrix of size N*N. Overwrite it with its LU factorization.
begin scalar w;
for i:=1:n do begin
scalar ii,pivot;
if w eq 'singular then return nil;
ii:=i;
while n>=ii and ((pivot:=getm2(a,ii,i))=0
or iremainder(pivot,prime!-base)=0) do ii := ii+1;
if ii>n then return(w := 'singular);
if not(ii=i) then begin
scalar temp;
temp:=getv(a,i);
putv(a,i,getv(a,ii));
putv(a,ii,temp) end;
putm2(a,i,0,ii); % Remember pivoting information;
pivot:=modular!-reciprocal pivot;
putm2(a,i,i,pivot);
for j:=i+1:n do
putm2(a,i,j,modular!-times(pivot,getm2(a,i,j)));
for ii:=i+1:n do begin
scalar multiple;
multiple:=getm2(a,ii,i);
for j:=i+1:n do
putm2(a,ii,j,modular!-difference(getm2(a,ii,j),
modular!-times(multiple,getm2(a,i,j)))) end end;
return w
end;
symbolic procedure back!-substitute(a,v,n);
% A is an N*N matrix as produced by LU-FACTORIZE-MOD-P, and V is a
% vector of length N. Overwrite V with solution to linear equations.
begin
for i:=1:n do
begin scalar ii;
ii:=getm2(a,i,0); % Pivot control;
if ii neq i
then begin scalar temp;
temp:=getv(v,i);
putv(v,i,getv(v,ii));
putv(v,ii,temp)
end
end;
for i:=1:n do begin
putv(v,i,times!-mod!-p(!*n2f getm2(a,i,i),getv(v,i)));
for ii:=i+1:n do
putv(v,ii,difference!-mod!-p(getv(v,ii),
times!-mod!-p(getv(v,i),!*n2f getm2(a,ii,i)))) end;
% Now do the actual back substitution;
for i:=n-1 step -1 until 1 do
for j:=i+1:n do
putv(v,i,difference!-mod!-p(getv(v,i),
times!-mod!-p(!*n2f getm2(a,i,j),getv(v,j))));
return v
end;
endmodule;
end;