% Copyright, J. Maurice Rojas; 2/11/09
% (Any use of this code within other software must acknowledge the
% author. Use within commercial software is permitted only with
% the author's consent.)
% This code decides positive feasibility for honest n-variate
% (n+2)-nomials (checking whether the inputs are honest beforehand),
% via Algorithm 3.2 from [Bihan-Rojas-Stella, 2009].
%
% Warning: This implemenation does NOT use gcd-free bases or
% high-brow linear forms in logs estimates: it merely
% checks linear forms in logs via Matlab's arithmetic.
% Nevertheless, this implementation has worked well on
% numerous examples. Comments and questions are welcome:
% rojas@math.tamu.edu
%
% usage:
% cktposfeas(a,coeffs)
% where
% a = n x m support matrix
% coeffs = 1 x m coefficient array
% tol = small positive number for zero threshold...
% ...and a defines a circuit (degenerate allowed). (In essence,
% this means that a is the union of a point with the vertex set
% of a simplex.) The output is a (hopefully true) declaration as
% to whether f(x) := c_1 x^{a_1} + c_2 x^{a_2} + ... + c_m x^{a_m} has a
% positive root or not. The complexity is polynomial in n+log(deg(f)),
% unlike older methods that were polynomial in deg(f)^n.
function cktposfeas(a,coeffs,tol)
% set up the support
n=size(a,1); m=size(a,2);
% throw out supports that are too large...
if m>n+2
disp(sprintf('Your support has too many points!'));
return;
end;
% keep track of coefficient signs...
posc=[]; pc=0;
zerc=[]; zc=0;
negc=[]; ngc=0;
for j=1:m
if coeffs(j)>0
posc=[posc j]; pc=pc+1;
elseif abs(coeffs(j))0
disp(sprintf('Zero coefficients not allowed! Please re-enter your polynomial in the correct format.'));
end;
% finish early if there are no sign alternations...
if (pc==0)|(ngc==0)
disp(sprintf('Your polynomial has NO positive roots since it''s coefficients are all of the same sign. Duh!'));
return;
end;
% unique difference of sign made negative...
if pc==1
coeffs=-coeffs;
ngc=pc;
negc=posc;
end;
% compute unique circuit relation (or see that you do not have a circuit)...
ahat=[ones(1,m); a];
b=null(ahat,'r');
if size(b,2)>1
disp(sprintf('Your support is similar to an (n''+k)-set in R^{n''} for k>2, and is tthus NOT a circuit...'));
return;
end;
% now let's sort and use my fancy algorithm...
% sort exponents AND coefficients using the b_i as keys. in particular,
% we group by positive b_i, negative b_i, then vanishing b_i...
pos=[]; p=0;
zer=[]; z=0;
neg=[]; ng=0;
for j=1:m
if b(j)>0
pos=[pos j]; p=p+1;
elseif abs(b(j)) change indices of negative circuit rel coords...
coeffs=coeffs(ind(:)); % same shuffle for coeffs...
negc=inv(negc(:)); % shuffle => change indices of negative coeffs...
if (ng==1)&(ngc==1)&(neg==negc)
jp=neg; % j' from my algorithm in [BRS09]...
% log(left-hand side of relations of Steps (b) & (c), Algor 3.2, pg 7, [BRS09])
le=log([-b(jp) coeffs(1:(jp-1))])*[-b(jp);b(1:(jp-1))];
% log(right-hand side of relations of Steps (b) & (c), Algor 3.2, pg 7, [BRS09])
ri=log([-coeffs(jp) b(1:(jp-1))'])*[-b(jp);b(1:(jp-1))];
if (jpri
disp(sprintf('Your polynomial has NO positive roots, and stably so.'));
return;
end;
end;
%
disp(sprintf('Your polynomial has at least one positive root!'));