##
##    Title   : Minimal decomposition of a polynomial in 1 variable as
##              a sum of powers.
##    Created :-Tue Jun  7 19:02:30 1994
##   Modified : Tue Mar 11 14:33:58 1997
##    Author  : Bernard Mourrain  <mourrain@korrigan.inria.fr>
##    Use     : > binarydec(x^4+x+1);-
##            : > binarydec(x^2*y+x*y^2);
##            : > binarydec(x^4+y^4);
##            : > ?binarydec
######################################################################
fdiscrim := proc(P,x)
local M,l;
   if P <> 0 then 
     M := linalg[sylvester](expand(P),diff(expand(P),x),x);
     l := convert(evalf(Svd(M)),list);
     l[nops(l)];
   else
     0
   fi;
end:

eker :=proc(M)
local i,v;
  v:= convert(linalg[linsolve](M,[0$linalg[rowdim](M)]),list);
  for i from 1 to linalg[rowdim](M) do v := subs(`_t`[i]=i,v) od;
  eval(v);
end:

lsqrker :=proc(M)
local i,v;
  v:= convert(linalg[linsolve](M,[0$linalg[rowdim](M)]),list);
  for i from 1 to rowdim(M) do v := subs(`_t`[i]=i,v) od;
  eval(v);
end:
######################################################################
# Change these definitions for approximate coefficients. 

DISCRIM := op(discrim): #use fdisrim <-> discrim
Kernel := op(eker):     #use lsqrker <-> eker

isnull :=proc(eps)
  evalb(eps=0); #evalb(abs(eps) <.00001); 
end:
######################################################################
circulant := proc(p,r,d)
local x,l,i,j,m;
  if nops(indets(p)) <>1 then ERROR(`polynomial in one variable please`); fi;
  x := op(1,indets(p)); 
  #  d := degree(p);
  l:=[0$d+1];
  for i from 1 to nops(l) do 
     l :=subsop(i=coeff(p,x,i-1)/binomial(d,i-1),l);
  od;
  m := array(0..d-r,0..r);
  for i from 0 to d-r do
    for j from 0 to r do m[i,j] := l[1+i+j] od; od;
  convert(m,matrix);
end:
######################################################################
binarydec := proc(P)
local x,d,r,i,v,p,m,G,H,q,LAMBDA,l,qq,ds,s;
  d := degree(P); 

  p := subs(y=1,expand(P));
  if nops(indets(p))<>1 then ERROR(`polynomial in x,y please`); fi;
  x := op(1,indets(p)); 
  # recherche du plus petite nombre de puissances.
  r := 1;
  m := circulant(p,r,d);
  v := Kernel(m);
  G :=0;
  for i from 1 to nops(v) do G := G+ x^(i-1)*v[i] od;
  ds :=  DISCRIM(G,x) ;

  while isnull(ds) and r < d  do
      r := r+1; 
      m := circulant(p,r,d);
      v := Kernel(m);
      G :=0;
      for i from 1 to nops(v) do G := G+ x^(i-1)*y^(r-i+1)*v[i] od;
      ds :=  DISCRIM(subs(y=1,G),x) ;
  od;
  # calcul des solutions complexes de G.
  H := G;
  G := subs(y=1,G);
  if(degree(H)=degree(G)) then 
    q := 0; 
    LAMBDA := {};
  else    
    q := u*x^d;
    LAMBDA := {u};
  fi;
  l := [fsolve(G,x,complex)];

  for i from 1 to nops(l) do 
    LAMBDA := {op(LAMBDA),lambda[i]};
    q := q + lambda[i]*(x*l[i] +1)^d
  od;
  qq := expand(q);
  s := {coeffs(qq-p,x)};
  if nops(indets(P))=2 then
    if degree(H)=degree(G) then 
      q :=  sum('lambda[j]*(x*l[j] +y)^d',j=1..nops(l));
    else 
      q := u*x^d + sum('lambda[j]*(x*l[j] +y)^d',j=1..nops(l));
    fi;
  fi;
  subs(linalg[leastsqrs](s,LAMBDA), q);
end:
######################################################################
`help/text/binarydec` := TEXT(` `,
` `,
`FUNCTION: binarydec - compute the minimal decomposition of a binary,`,
`          polynomial as a sum of powers of linear forms `,
`          see  Decomposition of quantics in sums of powers`,
`          P. Comon, and B. Mourrain,Signal Processing, Vol 53, Nb 2, p 
93-107`,
` `,
`CALLING SEQUENCE:`,
`   binarydec(p);`,
` `,
`PARAMETER:`,
`   p - polynomial in x or homogeneous in x and y`,
` `,
`EXAMPLE:`,
` > binarydec(expand((2*x+1)^6+(x+4)^6));`,
`                                      6                                  6`,
`       4096.000001 (.2500000000 x + 1)  + 1.000000000 (2.000000000 x + 1)`,
` `,
` > binarydec(y^3+x^2*y);`,
`                                       3                                  3`,
`     .4999999999 ( - .5773502692 x + y)  + .4999999999 (.5773502692 x + y)`,
` `,
`> binarydec((x+1)^4+x^2);`,
`                                                                              
4`,
`   ( - .01102847626 - .02369918480 I) (( - .7591289466 - 1.212312388 I) x + 
1)`,
`                                                                              
4`,
` + ( - .01102847626 + .02369918480 I) (( - .7591289466 + 1.212312388 I) x + 
1)`,
`                                  4`,
` + 1.022056953 (1.018257893 x + 1)`,
` `):
######################################################################
lprint(`binarydec := proc(p) ... end`);