FFGE := proc(A::Matrix(polynom(integer)), b::Vector(polynom(integer)))
#
# Solve Ax=b for x Lipson's fraction-free Gaussian elimination
# Output det(A), f, g where x[i] = f[i]/g[i] and gcd(f[i],g[i]) = 1
#
# The entries of A and b must be polynomials over Z
# The entries of f and g are polynomials over Z
#
# Author Michael Monagan, 2025
#
local n,m,B,mu,i,j,k,det,x,y,num,r,numterms,f,g,h;
   numterms := proc(f) if f=0 then 0 elif type(f,`+`) then nops(f) else 1 fi end;
   n,m := op(1,A);
   if n<>m then error "Matrix must be square" fi;
   m := op(1,b);
   if m<>n then error "Matrix and vector must have the same dimension" fi;
   B := <A|b>;
   mu := 1:
   det := 1;
   for k to n-1 do
       i := k;
       while i<=n and B[i,k]=0 do i := i+1; od;
       if i>n then return 0 fi;
       if i>k then # interchange row i and k
          for j from k to n+1 do B[i,j],B[k,j] := B[k,j],B[i,j] od;
          det := -det;
       fi;
       for i from k+1 to n do
           for j from k+1 to n+1 do
               num := expand(B[k,k]*B[i,j]-B[i,k]*B[k,j]);
               divide(num,mu,evaln(B[i,j]));
               if i=k+1 and j=k+1 then 
                  printf("i=%d  #num=%d  #B[i,j]=%d\n",i,numterms(num),numterms(B[i,j]))
               fi;
           od;
           B[i,k] := 0;
       od;
       mu := B[k,k];
   od;
   det := det*B[n,n];
   printf("#det=%d\n",numterms(det));
   # Lipson's back substitution
   y := Vector(n);
   y[n] := B[n,n+1];
   printf("#y[%d]=%d\n",n,numterms(B[n,n+1]));
   for i from n-1 by -1 to 1 do
       num := expand(B[i,n+1]*B[n,n]-add(B[i,j]*y[j],j=i+1..n));
       divide(num,B[i,i],evaln(y[i]));
       printf("#N[%d]=%d  #y[%d]=%d\n",i,numterms(num),i,numterms(y[i]));
   od;
   f := Vector(n);
   g := Vector(n);
   for i from 1 to n do
       h := gcd(y[i],B[n,n],evaln(f[i]),evaln(g[i]));
       printf("#f[%d]=%d #g[%d]=%d\n",i,numterms(f[i]),i,numterms(g[i]));
   od;
   return det,f,g;
end:

with(LinearAlgebra):
A := ToeplitzMatrix([y1,y2,y3],symmetric);
b := <1,1,1>;
d,f,g := FFGE(A,b):
factor(d);
seq( x[i]=f[i]/g[i], i=1..3 );

A := ToeplitzMatrix([y1,y2,y3,y4],symmetric):
b := <1,1,1,1>:
d,f,g := FFGE(A,b):
factor(d);
for i to 4 do x[i]=f[i]/g[i] od;