## ## 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 ## 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`);