|\^/| Maple 2026 (X86 64 LINUX) ._|\| |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2026 \ MAPLE / All rights reserved. Maple is a trademark of <____ ____> Waterloo Maple Inc. | Type ? for help. > 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 := ; > 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); memory used=6.0MB, alloc=40.3MB, time=0.04 [y1 y2 y3] [ ] A := [y2 y1 y2] [ ] [y3 y2 y1] > b := <1,1,1>; [1] [ ] b := [1] [ ] [1] > d,f,g := FFGE(A,b): i=2 #num=2 #B[i,j]=2 i=3 #num=4 #B[i,j]=4 #det=4 #y[3]=4 #N[2]=8 #y[2]=4 #N[1]=4 #y[1]=4 #f[1]=2 #g[1]=3 #f[2]=3 #g[2]=3 #f[3]=2 #g[3]=3 > factor(d); 2 2 (y1 - y3) (y1 + y1 y3 - 2 y2 ) > seq( x[i]=f[i]/g[i], i=1..3 ); y1 - y2 y1 - 2 y2 + y3 x[1] = -------------------, x[2] = -------------------, 2 2 2 2 y1 + y1 y3 - 2 y2 y1 + y1 y3 - 2 y2 y1 - y2 x[3] = ------------------- 2 2 y1 + y1 y3 - 2 y2 > A := ToeplitzMatrix([y1,y2,y3,y4],symmetric): > b := <1,1,1,1>: > d,f,g := FFGE(A,b): i=2 #num=2 #B[i,j]=2 i=3 #num=4 #B[i,j]=4 i=4 #num=19 #B[i,j]=12 #det=12 #y[4]=13 #N[3]=47 #y[3]=16 #N[2]=29 #y[2]=16 #N[1]=13 #y[1]=13 #f[1]=2 #g[1]=7 #f[2]=4 #g[2]=7 #f[3]=4 #g[3]=7 #f[4]=2 #g[4]=7 > factor(d); memory used=23.3MB, alloc=108.3MB, time=0.16 2 2 2 (y1 + y1 y2 + y1 y4 - y2 - 2 y2 y3 + y2 y4 - y3 ) 2 2 2 (y1 - y1 y2 - y1 y4 - y2 + 2 y2 y3 + y2 y4 - y3 ) > for i to 4 do x[i]=f[i]/g[i] od; y1 - y3 x[1] = ------------------------------------------------- 2 2 2 y1 + y1 y2 + y1 y4 - y2 - 2 y2 y3 + y2 y4 - y3 y1 - y2 - y3 + y4 x[2] = ------------------------------------------------- 2 2 2 y1 + y1 y2 + y1 y4 - y2 - 2 y2 y3 + y2 y4 - y3 y1 - y2 - y3 + y4 x[3] = ------------------------------------------------- 2 2 2 y1 + y1 y2 + y1 y4 - y2 - 2 y2 y3 + y2 y4 - y3 y1 - y3 x[4] = ------------------------------------------------- 2 2 2 y1 + y1 y2 + y1 y4 - y2 - 2 y2 y3 + y2 y4 - y3 > quit memory used=23.3MB, alloc=108.3MB, time=0.16