function [mu,Q,err]=binarydec3(p,eta)
% [mu,Q]=binarydec3(p) ou [mu,Q]=binarydec3(p,eta)
% % P.Comon, 25 Sept 1997
% % D'apres Comon-Mourrain, Signal processing'96
% % Decomposition de quantiques binaires reelles ou complexes
% % Difference avec binarydec2: prise en compte des racines DISTINCTES,
% % et prise de Re quand elles sont COMPLEXES. De plus, err est RELATIVE.
% % eta est un parametre de precision. Il vaut 1.e-4 par defaut
% % binarydec3 choisit le premier r tel que la decomposition
% % soit approchee a eta pres (dans la svd).
% % Exemple: 
% format('compact');
% path(path,'/homes/comon/Matlab/Fonctions')
% p1=[2,1];p2=[1 -4];p3=[1 1]*1.e-7;
% p=conv4(p1)+conv4(p2)+conv4(p3);
% [mu,Q,err]=binarydec3(p); err
% % On a alors: p = sum_k mu(k) conv4(q(k,:)). Soit:
% % err affichee donne l'erreur normalisee de reconstruction
% % Les formes lineaires estimees sont donnees par SOL.
% % Autres exemples (erreur de reconstruction doit tjs etre faible):
% % ORDRES PAIRS (non unicite des formes):
% p=conv4([2 1])+conv4([1 -3])+conv4([0.1 -.1]);
% p=conv4([2 1])+conv4([1 -3])+ones(1,5)*0.01;
% p=convd([2,1],6)+convd([3,-5],6)+convd([-2,1],6)+convd([1,1],6);
% % ORDRES IMPAIRS (unicite des formes a 1 facteur multiplicatif pres):
% p=convd([2,1],5)+convd([3,1],5)+convd([-2,1],5);
% p=convd([2,1],5)+convd([3,-5],5)+convd([-2,1],5);
% % EXECUTION ET AFFICHAGE DU RESULTAT
% [mu,Q]=binarydec3(p) 
% % NB: La reconstruction de p est donnee par:
% % d=length(p)-1;sol=(mu'*convd(Q,d));
% % On peut donc aussi comparer en regardant [sol' p']

cr=isreal(p); % decompositions reelles ou complexes
s=1;r=0;d=length(p)-1;
if nargin==1,eta=1.e-4;end;
fd=facto(d);c=ones(1,d+1);
for i=1:d-1,c(i+1)=fd/facto(i)/facto(d-i);end;
p0=p;p=p./c;v=[];
while s>eta&r<d-r+2,
  r=r+1;
  M=hankel(p(1:d-r+1),p(d-r+1:d+1));
  [U,S,V]=svd(M);
  J=find(diag(S)<eta);
    if length(J)>0,s=S(J,J);J=J(length(J));
    elseif r+1>d-r+1,J=r+1;
    end;
  if length(J)>0, % sortie de boucle si racines distinctes
  v=V(:,J);
  if cr==1,q=real(roots(v));else q=roots(v);end;
  Lq=length(q);tt=q(1:Lq-1)-q(2:Lq);It=find(abs(tt)<eta);
  if length(It)==0,s=0;end; 
  end;
end;
Q=[q,ones(length(q),1)];
mu=convd(Q,d)'\p0';

%% AFFICHAGE %%
sol=(mu'*convd(Q,d));W=diag(ones(1,d+1)./c);
err=sqrt((sol-p0)*W*(sol-p0)')/sqrt(p0*W*p0');
fprintf('nombre de formes = %g\n',r-1);
fprintf('err normalisee de reconstruction = %g\n',err);
%SOL=diag(mu.^(1/d))*Q;
%fprintf('FORMES LINEAIRES ESTIMEES:\n')
%[lig,col]=size(SOL);for j=1:lig,
%fprintf('Forme%g = ',j)
%  for k=1:col, fprintf(' %3.6f  ',SOL(j,k));end;
%  fprintf('\n');
%end;