with(linalg); invmat:=proc(A,l,p,n) local X,N,R,C,B,E; E:=array(1..n,1..n,identity); X:=Inverse(A) mod p; print(X); B:=evalm(l*X); N:=p; do R:= evalm((l*E-B&*A)/N); if (nullmat(R,n)) then RETURN(B); fi; C:=evalm(R&*X); B:=evalm(B+N*C); N:=p*N; for i from 1 to n do for j from 1 to n do B[i,j] := B[i,j] mod N; if (B[i,j] > N/2) then B[i,j] := B[i,j] - N; fi; od; od; od; end; nullmat:=proc(R,n) local i,j; for i from 1 to n do for j from 1 to n do if (R[i,j] <> 0 ) then RETURN(false); fi; od; od; RETURN(true); end; ratapp:=proc(c,N) local n1,n2,s,v1,v2,h,k; v1:=array(1..2,[c,1]); v2:=array(1..2,[N,0]); do n1:=(v1[1]^2+v1[2]^2); n2:=(v2[1]^2+v2[2]^2); s:=(v1[1]*v2[1]+v1[2]*v2[2]); if (n1>n2) then if (abs(s) > n2/2) then k:= round(s/n2); v1:=evalm(v1-k*v2); else if (v2[2]<0) then v2:=evalm(-v2); fi; RETURN(v2); fi; else if (abs(s) > n1/2) then k:= round(s/n1); v2:=evalm(v2-k*v1); else if (v1[2]<0) then v1:=evalm(-v1); fi; RETURN(v1); fi; fi; od; end;