function m_eig
global n
global m
global ALLSEQ
global ALLSEQ1
m = 4;
n = 6;

global A C

for t=1:1000
  
  A = pnormdim(rand(m,m),1);
  %C = pnormdim(rand(m,1),1);
  C = zeros(m,1); C(end)=1;
  %load tmp
  %A = pnormdim(A+realmin,1);
  
  E = sort(abs(eig(A)));
  h = E(end-1);
  
  i0=1;
  ALLSEQ =  list_all_ind(repmat(m,1,n));
  ALLSEQ1 = list_all_ind(repmat(m,1,n-i0+1));
  XXi = list_all_ind(repmat(m,1,i0));
  YY = [];
  %for j=1:size(XXi,1)
  for j=1
    Xi = XXi(j,:);
    YY(j) = Y(Xi);
    %YY_pr(j) = Y_pr(Xi);
  end
  YY
  %if (totdif(YY,YY_pr)>1e-10)
  %  'dont match up!!!!'
  %  keyboard
  %end
  C = sum(h.^[0:n-1])

  g = all(abs(YY)<=C+eps);
  if ~g
    'eig dont work'
    keyboard
  end
  fprintf('good %d \n',t);
end

return


function y = h(x)
y = (x==1) + 2*(x~=1);
return

function r = F(x)
r = length(find(x~=1));
return

function Ef = E(Xi)
% E[f(X)|Xi]
global ALLSEQ
global A
Ef = 0;
i = length(Xi);
n = size(ALLSEQ,2);
for k=1:size(ALLSEQ,1)
  x = ALLSEQ(k,:);
  p = P(x);
  if (i & (i<n))
    pxi = A(x(i+1),x(i));
    pXi = A(x(i+1),Xi(i));
    p = p / pxi * pXi;
  end
  f = F([Xi x(i+1:end)]);
  Ef = Ef + f*p;
end
return


function y = Y(Xi)
% Y1 = E[f(X)|X_i] - E[f(X)|X_{i-1}]
Ef1 = E(Xi);
Ef0 = E(Xi(1:end-1));
y = Ef1 - Ef0;
return


function y = Y_pr(Xi)
% Y1 = E[f(X)|X_i] - E[f(X)|X_{i-1}]
global ALLSEQ1
global A
i = length(Xi);
%n = size(ALLSEQ1,2);
y = 0;

for j=1:size(ALLSEQ1,1)
  x = ALLSEQ1(j,:);
  x = [x 0];
  Xix = [Xi(1:i) x(2:end)];
  Xi_1x = [Xi(1:i-1) x];
  a = pcp(Xix(i+1),Xix(i));
  b = pcp(Xi_1x(i+1),Xi_1x(i));
  p = PCOND(x,Xi(1:i-1));
  p = p / b;

  f0 = F(Xi_1x);
  f1 = F(Xix);
  y = y + p*(a*f1-b*f0);
  %show(X,p,alph,f1,f0);
  %show(X,p,A(X(2),X1),A(X(2),X(1)),f1,f0,X1);
  %show(x(1:end-1),(a*f1-b*f0));
end
return

function p = P(x)
global A
global C
if isempty(x)
  p = 1;
else
  p = C(x(1));
end
for j=2:length(x)
  p = p * A(x(j),x(j-1)); % markov case
  %p = p * C(x(j)); % independent case
end
return

function p = PCOND(x,Xi)
global C
if isempty(Xi)
  p = C(x(1));
else
  p = pcp(x(1),Xi(end));
end
for j=2:length(x)
  p = p * pcp(x(j),x(j-1)); % markov case
end
return

function p = pcp(i,j)
% primitive cond prob
global A
if i
  p = A(i,j);
else
  p = 1;
end
return

function r = randinseg(a,b)
r = rand*(b-a) + a;
return