Code covered by the BSD License

groebner

Ben Petschel (view profile)

19 Jun 2009 (Updated )

manipulate and solve systems of multivariate polynomial equations by computing the groebner basis

gbasis=groebner(polyset,ord,varnames,tol)
```function gbasis=groebner(polyset,ord,varnames,tol)
% GROEBNER - calculate the reduced Groebner basis of a set of polynomials
%  usage: gbasis = groebner(polyset,ord)
%         gbasis = groebner(polyset,ord,varnames)
%         gbasis = groebner(polyset,ord,varnames,tol)
%
% INPUTS: polyset, ord, varnames, tol
%  polyset is a cell array of multivariate polynomial coefficients
%    the cell elements can be all strings or all arrays.
%    If the polynomials are specified with arrays, each row of the array
%    represents a single term.  The first element is the coefficient and
%    the others are the degrees of each unknown in the term.
%    e.g. [-1,2,1;1,0,1] represents -x1^2*x2+x2
%    If the polynomials are specified as strings, they must use variable
%    names 'x1','x2',... (unless varnames is provided) and be valid inputs
%    for str2poly.m
%  ord is a string specifying the monomial ordering:
%    'lex': lexicographical, order by highest power of most significant
%      indeterminate, e.g. x1>x2^2, x1*x2>x1.
%    'grlex': graded lex, order by total degree then lex,
%      e.g. x1^2*x2 > x1*x2^2 > x1^2
%    'grevlex': graded reverse lex, order by total degree then by lowest
%      power of least significant indeterminate, e.g. x1*x2^2 > x1*x2*x3
%    (for all cases have x1>x2>...>xn)
%    note: revlex is not a well-ordering because 1>x1>x1^2>...
%  varnames (optional) is a cell array of variable names if polyset is a
%    cell array of strings and the variable names are not 'x1', 'x2', etc.
%    If specified, there must always be at least 2 variables.
%  tol (optional) is the zero tolerance for coefficients, otherwise
%    catastrophic cancellation may occur.  Default value is 0
%
% OUTPUTS: gbasis
%  gbasis is the reduced Groebner basis of the set of polynomials, in the
%    same format (strings or matrices) as polyset{1}.
%
% ALGORITHM:
%  A modified Buchberger's algorithm is used.  In Buchberger's algorithm,
%  at each step calculate Sij = (Lij/ai)*Pi - (Lij/aj)*Pj, for all pairs of
%  polynomials Pi, Pj, where ai and aj are the leading terms of Pi and Pj,
%  and Lij is the least common multiple of ai and aj; then reduce Sij with
%  respect to the polynomials and add it to the set (note that Sij always
%  reduces to zero if ai and aj have no variables in common).  The
%  algorithm terminates when no new polynomials can be added to the set.
%  This function reduces the polynomials at each step in order to find the
%  reduced Groebner basis.
%
% EXAMPLE: simplify equations x^2+2xy^2=0, xy+2y^3=1
%  A1=[1,2,0;2,1,2]; % 1*x^2*y^0 + 2*x^1*y^2
%  A2=[-1,0,0;1,1,1;2,0,3]; % -1*x^0*y^0 + 1*x^1*y^1 + 2*x^0*y^3
%  y=groebner({A1,A2},'lex');
%  poly2str(y) % returns {'x2^3-0.5', 'x1'}
%
% EXAMPLE: same, but with equations specified as strings:
%  groebner({'x^2+2*x*y^2','x*y+2*y^3-1'},'lex',{'x','y'})
%  % returns {'y^3-0.5','x'}
%
% EXAMPLE: show that a set of equations is inconsistent
%  groebner({'x + y', 'x^2 - 1', 'y^2 - 2*x'}, 'lex', {'x', 'y'})
%  % returns {'1'}
%
%
% KNOWN ISSUES
%
%  0. GENERAL: while a significant effort has been made to ensure that the
%    program works as described, unidentified bugs may remain and the results
%    should be regarded with caution.  In particular, the suite of tests
%    performed was by no means comprehensive.
%
%  1. ACCURACY: the algorithm by default works in floating-point which
%    ultimately cannot perform exact arithmetic, so unexpected silent errors may
%    occur even when the coefficients are specified as integers.  A workaround
%    is to use other data types but this has its own issues (see below).
%
%  1a. EXAMPLE (roundoff - linearly independent polynomials)
%    eps % returns 2.2204e-16 (32-bit PC windows machine epsilon)
%    groebner({'x1+x2','x1+1.000000000000001*x2'},'lex') % ok, returns {'x2','x1'}
%    groebner({'x1+x2','x1+1.0000000000000001*x2'},'lex') % not ok, returns {'x2+x1'}
%
%  1b. EXAMPLE (roundoff - linearly dependent polynomials)
%    2/sqrt(2)-sqrt(2) % expect 0, returns -2.2204e-016
%    poly2str(groebner({[1,1,0;sqrt(2),0,1],[1/sqrt(2),1,0;1,0,1]},'lex')) % ok, returns {'1.41421*x2+x1'}
%    poly2str(groebner({[1,1,0;sqrt(2),0,1],[sqrt(2),1,0;2,0,1]},'lex')) % not ok, returns {'x2','x1'}
%    poly2str(groebner({[1,1,0;sqrt(2),0,1],[sqrt(2),1,0;2,0,1]},'lex',{},1e-14)) % ok, returns {'1.41421*x2+x1'}
%
%  1c. EXAMPLE (courtesy of Christophe Lauwerys - catastrophic cancellation & tolerance bug)
%    groebner({'t^3+x+y','t^2+0.5*x^2-x-z^2','t^2+y-z^2'},'lex',{'t','x','y','z'}) % ans{1} contains Inf
%
%  1d. EXAMPLE (from Mathematica documentation - catastrophic cancellation because of large intermediate coefficients)
%    % expect 3 polynomials returns, one of degree 8 in z; instead can get '1'
%    groebner({'x^2 + y^2 + z^2 - 1', 'x*y - z + 2', 'z^2 - 2*x + 3*y'}, 'lex', {'x', 'y','z'})
%
%  2. NON-DOUBLE INPUTS: use of non-double coefficients is not supported.  While
%    it may be tempting to introduce symbolic coefficients, it is not clear how
%    they should be consistently treated - for example 'A*x' could reduce to '0'
%    or 'x' or both, depending upon the assumptions on the coefficient A.
%
%  3. EFFICIENCY: The reduction algorithms aren't terribly efficient.  One of
%    the major bottlenecks at present is the calculation of the leading term,
%    which currently involves a call to sortrows each time.  It could be more
%    efficient to maintain an ordered list or even a priority queue (e.g.
%    implemented as a heap).  Also lattice reduction algorithms such as LLL
%    might speed things up.
%
%
%  poly2str, str2poly, polynsolve

% Author: Ben Petschel 19/6/2009
%
% Change history:
%  19/6/2009 - first release (ord='grlex' and ord='grevlex' not
%    implemented yet)
%  22/6/2009 - implemented ord='grlex' and ord='grevlex'
%    reduction algorithm now reduces lower-order terms as well
%  23/6/2009 - removed ord='revlex' because it is not a well-ordering
%  17/7/2009 - fixed bug in handling polynomials with >=3 variables
%  20/3/2010 - changed polynomial representation from multidimensional
%    arrays to 2d arrays (more memory efficient with large number of vars)
%  11/10/2010 - bugfix for handling error tolerance

if (numel(polyset)>0) && ischar(polyset{1}),
charout = true;
if nargin<3,
polyset = str2poly(polyset);
else
polyset = str2poly(polyset,varnames);
end;
else
charout = false;
end;

if nargin<4,
tol = 0;
end;

% determine the number of degrees
nd=max(cellfun(@(x)size(x,2)-1,polyset)); % error if = 0
if nd == 0
error('number of variables must be >= 1');
end;

gbasis = fullreduce(polyset,ord,tol,nd);
oldgbasis = {};

while ~isequal(oldgbasis,gbasis),
oldgbasis = gbasis;
Qset = SPoly(gbasis,ord,tol,nd);
gbasis = fullreduce([gbasis,Qset],ord,tol,nd);
end;

if charout,
if nargin<3,
gbasis = poly2str(gbasis);
else
gbasis = poly2str(gbasis,varnames);
end;
end;

end % main function groebner(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Qset = SPoly(Pset,ord,tol,nd)
% calculates all reduced S-polynomials of Pset (Pset must be reduced)
%   Sij = (Lij/ai)*Pi - (Lij/aj)*Pj,
% for all pairs of polynomials Pi, Pj, where ai and aj are the leading
% terms of Pi and Pj, and Lij is the least common multiple of ai and aj
% then reduces Sij wrt Pset

Qset = {};
for i=1:numel(Pset)-1,
for j=i+1:numel(Pset),
if any(and(d1>0,d2>0)),
% leading terms have a variable in common, so S-poly is nontrivial
L = max(d1,d2); % LCM of leading terms
S = reduceset(S,Pset,ord,tol,nd);
%if any(S(:)~=0),
Qset{end+1}=S; % add S to the list
end;
end;
end;
end;

end % helper function SPoly(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Qset = fullreduce(Pset,ord,tol,nd)
% reduces all polynomials wrt each other
Qset = Pset;
oldQset = {};

while ~isequal(oldQset,Qset),
oldQset = Qset;
for i=1:numel(Qset)-1,
for j=i+1:numel(Qset),
Qset{i}=reduce(Qset{i},Qset{j},ord,tol,nd);
Qset{j}=reduce(Qset{j},Qset{i},ord,tol,nd);
end;
end;
end;

% keep only the nonzero results (to within zero tolerance tol)

if numel(Qset)==1,
% special case, a single equation can slip through without being reduced
% coefficients occupy first column of the array
Qset{1} = multconst(Qset{1},1/c); % know c~=0 because have already removed zero eqns
end;

end % helper function fullreduce(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q = reduceset(P,Pset,ord,tol,nd)
% reduces P wrt all polynomials in Pset
Q = P;
oldQ = [];

while ~isequal(oldQ,Q),
oldQ = Q;
for i=1:numel(Pset),
Q=reduce(Q,Pset{i},ord,tol,nd);
end;
end;

end % helper function reduceset(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q1=reduce(P1,P2,ord,tol,nd)
% reduces a polynomial P1 by subtracting multiples of P2 so that leading
% term coefficient is 1 and leading term of P2 does not divide that of P1

Q2=P2;
if abs(c2)<=tol,
% ignore terms less than tol
keepgoing = false;
Q1 = P1;
else
keepgoing = true;
Q2 = multconst(Q2,1/c2);
Q1 = 0; % will successively add terms to Q
% remove the leading terms, in order to reduce wrt lower terms
remain = P1;
end;

while keepgoing,
if abs(c1)<=tol,
% ignore terms less than tol
keepgoing = false;
else
%reduce remain by Q2, if possible, or remove leading term
if d1>=d2,
% can reduce Q1 by Q2, so subtract a multiple of Q2 from Q1
else
% cannot reduce any further wrt leading term, but try lower terms
end;
end;
end;

% reduce leading coefficient, if possible
if abs(c1)>tol,
Q1 = multconst(Q1,1/c1);
end;

end % helper function reduce(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=multconst(P,c)
% returns c*P where c is a constant

Q = P;
Q(:,1)=Q(:,1)*c;

end % helper function multconst(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s1 = size(P1,2);
s2 = size(P2,2);

if s1>s2
Q = [P1;P2,zeros(size(P2,1),s1-s2)];
else
Q = [P1,zeros(size(P1,1),s2-s1);P2];
end;

Q = sortrows(Q,2:size(Q,2));
i=1;
while i<size(Q,1)
% merge rows representing terms with the same exponents
if Q(i,1)==0,
Q(i,:)=[];
elseif all(Q(i,2:end)==Q(i+1,2:end))
Q(i,1)=Q(i,1)+Q(i+1,1);
Q(i+1,:)=[];
if Q(i,1)==0
% remove terms that add to zero
Q(i,:)=[];
end;
else
i=i+1;
end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=multiplyterm(P,d)
% multiplies a polynomial by a term (no coefficient)

% padding a matrix at the end is same as multiplying by x
d=d(:).';
sP = size(P,2);
sd = size(d,2)+1; % d doesn't include coefficient
if sP<sd,
Q=[P,zeros(size(P,1),sd-sP)];
else
Q=P;
d=[(d(:)).',zeros(1,sP-sd)];
end;
Q=Q+repmat([0,d],size(Q,1),1);

end % helper function multiplyterm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% returns the leading term of P
% deg is a row vector of length size(P,2)-1
% coeff is the leading coefficient, deg(i) is the degree of xi
% ignores terms with coefficients <= tol
% nd is the minimum number of degrees (in order to handle cases where deg=0)

d=size(P,2);

if all(abs(P(:,1))<=tol),
% P=0
coeff=0;
deg=zeros(1,d-1);
else
switch ord
case 'lex'
P=sortrows(P,2:d);

case 'revlex'
error('revlex is not a well-ordering');

case 'grlex'
P = [P,sum(P(:,2:d),2)];
P = sortrows(P,[d+1,2:d]);
P = P(:,1:d);

case 'grevlex'
P = [P,sum(P(:,2:d),2)];
P = sortrows(P,[d+1,-(d:-1:2)]);
P = P(:,1:d);

otherwise
error('ordering %s not recognized',ord);

end; % switch ord

ind = find(abs(P(:,1))>tol,1,'last');
if isempty(ind),
coeff = 0;
deg = zeros(1,d-1);
else
coeff = P(ind,1);
deg = P(ind,2:end);
end;

end; % if all(P(:)==0)

if nd>d-1,
deg = [deg,zeros(1,nd-d+1)];
end;