function L=snepsolver(W,e)
%------------------------------------------
%input: W--polynomial system 
%       e--a small tolerance
%output: 
%       P--the table of dimensions of projection and prolongtion matrices
%       M--the multiplication matrices
%       L--the solutions of polyonmial system
%-------------------------------------------
format long
P=zeros(10);
W=pre_prolongtion(W);%將初給的多項式系統化為齊次，使各個多項式的最高次相等,且為原多項式系統的最高次數
deg=ldeg(W);
n=size(W,2)-2;
A=W;
flag=0;
for j=1:10
    %for語句得到多項式系統A經過若干次prolongtion和projection算子得到的維數表P
    if j==1
        B=A;
    else
        B=prolongtion(A,j-1);
    end
        s=ldeg(B);
        C=convert(B);
        [U,S,V]=svd(C);
        rank=ndim(S,e);
        m=size(V,2);
        t=0;
    for i=1:j+deg
        if i==1
            P(i,j)=size(V,1)-rank;
        else
            t=t+kdeg(s,n);
            D=V(t+1:m,rank+1:m);
            [U,S,F]=svd(D);
            P(i,j)=ndim(S,e);
            s=s-1;
        end
        if i>=2&&j>=2&&P(i-1,j-1)==P(i,j)&&P(i-1,j-1)==P(i,j-1)
            %此時多項式系統對合
            k=j-1;
            flag=1;
        end   
    end
    if flag==1
        break;
    end 
end
P,       %輸出prolongtion和projection算子得到的維數表P
for i=k:-1:1
    if P(i,k)==P(i+1,k)&&P(i,k)==P(i+1,k+1)
        break;
    end
end
l=i;        %在維數表P中選出滿足對合條件的最小的k(prolongtion次數),,,以及在這個k下最大的l(projection次數)
dim=P(i,k);
if k==1
    B=A;
else
    B=prolongtion(A,k-1);
end
C=convert(B);
[U,S,V]=svd(C);
rank=ndim(S,e);
m=size(V,2);
D=V(m-gsum(deg+k-l,n)+1:m,rank+1:m);    %D為對合的approximate spanning set子矩陣
E=V(m-gsum(deg+k-l-1,n)+1:m,rank+1:m);  %E為將D再進行一次projection得到approximate spanning set子矩陣
max=m-rank;
Q=app_basis(E,dim,max,e);
DD=zeros(size(D,1),dim);
EE=zeros(size(E,1),dim);
for i=1:dim
    DD(:,i)=D(:,Q(1,i));        %將D，E利用app_basis隨機數篩選基而得到的滿足對合條件的approximate basis DD，EE
    EE(:,i)=E(:,Q(1,i));
end
[U,S,V]=svd(EE);
P=zeros(dim);
for i=1:dim        %for語句得到各個multiplication matrix M，及multiplication matrix的線性組合T
    P(i,i)=1/S(i,i);
end
G=U(:,1:dim);
X=grlex(deg+k-l,n);
T=zeros(dim);
for i=1:n
    s=1;
    N=zeros(gsum(deg+k-l-1,n),size(V,2));
    for j=1:gsum(deg+k-l,n)
        if X(j,i)>0
            N(s,:)=DD(j,:);
            s=s+1;
        end
    end
    M=G'*N*V*P;
    M,            %輸出各個multiplication matrix M
    T=T+((-1)^i)*M;
end
[V,Y]=eig(T);
R=G*V;
L=zeros(n,dim);
for i=1:n         %利用multiplication matrix的線性組合T的特征向量計算出多項式系統的solutions L,輸出solutions L
    L(i,:)=R(size(R,1)-n+i-1,:)./R(size(R,1),:);
end
%--------------------------------------------
function A=pre_prolongtion(W)
%將初給的多項式系統化為齊次，使各個多項式的最高次相等,且為原多項式系統的最高次數
format long
Q=size(W,1);
n=size(W,2)-2;
m=W(Q,n+1);
F=zeros(m,2);
q=0;
max=0;
s=1;
for i=1:Q    %for語句統計出各個多項式的項數和最高次
    sum=0;
    for j=1:n
        sum=sum+W(i,j);
    end
    if W(i,n+1)==s
        q=q+1;
        if sum>max
            max=sum;
        end
    else
        F(s,1)=q;
        F(s,2)=max;
        q=1;
        max=sum;
        s=s+1;
    end
end
F(s,1)=q;
F(s,2)=max;
for i=1:m
    if F(i,2)>max
        max=F(i,2);
    end
end
sum=0;
for i=1:m
    sum=sum+F(i,1)*kdeg(max-F(i,2),n);
end
A=zeros(sum,n+2);
sum=0;
p=0;
flag=0;
for i=1:m                %for語句完成對各個多項式的延拓，使其最高次相等
    D=rlex(max-F(i,2),n);
    for j=1:kdeg(max-F(i,2),n)
        for t=1:F(i,1)
            p=p+1;
            A(p,1:n)=W(sum+t,1:n)+D(j,:);
            A(p,n+1)=flag+j;
            A(p,n+2)=W(sum+t,n+2);
        end
    end
    sum=sum+F(i,1);
    flag=flag+kdeg(max-F(i,2),n);
end
%--------------------------------------------------
function B=prolongtion(A,k)
%對多項式系統A進行k次prolongtion算子延拓
format long
flag=0;
Q=size(A,1);
n=size(A,2)-2;
m=A(Q,n+1);
q=Q;
B=zeros(Q*gsum(k,n),n+2);
B(1:Q,:)=A;
for i=1:k
    D=rlex(i,n);
    for j=1:kdeg(i,n)
        flag=flag+1;
        for t=1:Q
            q=q+1;
            B(q,1:n)=A(t,1:n)+D(j,:);
            B(q,n+1)=A(t,n+1)+flag*m;
            B(q,n+2)=A(t,n+2);
        end
    end
end
%-------------------------------------------------
function clum=app_basis(E,dim,max,e)
%將對合的approximate spanning set子矩陣E利用隨機數篩選出一組approximate basis
B=zeros(size(E,1),dim);
clum=sort(ceil(max*rand(1,dim)));
for i=1:dim
    B(:,i)=E(:,clum(1,i));
end
[U,S,V]=svd(B);
while ndim(S,e)~=dim
    clum=sort(ceil(max*rand(1,dim)));
    for i=1:dim
        B(:,i)=E(:,clum(1,i));
    end
    [U,S,V]=svd(B);
end
%-------------------------------------------------
function d=ndim(A,e)
%計算SVD分解的核空間的近似維數
r=size(A,1);
l=size(A,2);
if r>l
    min=l;
else
    min=r;
end
d=0;
for i=1:min
    if A(i,i)>=e
        d=d+1;
    end
end
%---------------------------------------------------
function C=convert(B)
%將延拓后的多項式系統B的指標矩陣轉化為系數矩陣
format long
Q=size(B,1);
n=size(B,2)-2;
deg=ldeg(B);
q=Q;
C=zeros(B(q,n+1),gsum(deg,n));
for i=1:q
    w=B(i,n+1);
    E=grlex(deg,n); 
    for j=1:gsum(deg,n)
        b=E(j,:)-B(i,1:n);
        a=b.*b;
        s=0;
        for t=1:n
            s=s+a(t);
        end
        if s==0
           v=j;
           break;
        end
    end
    C(w,v)=B(i,n+2);
end
C=flipud(C);
%-----------------------------------------------
function l=gsum(k,n)
%計算n個未知數的所有次數不超過k次的單項式的總數
l=0;
for i=0:k
    l=l+kdeg(i,n);
end
%----------------------------------------------
function g=grlex(k,n)
%n個未知數的所有次數不超過k次的單項式分次字典序的指標列表(從高到低)
sum=0;
for i=0:k
    sum=sum+kdeg(i,n);
end
g=zeros(sum,n);
r=1;
for i=k:-1:0
    g(r:r+kdeg(i,n)-1,1:n)=rlex(i,n);
    r=r+kdeg(i,n);
end
%----------------------------------------------
function output=rlex(k,n)
%n個未知數的所有k次單項式分次字典序的指標列表(從高到低)
output=zeros(kdeg(k,n),n);
q=1;
if n==1
   output=k;
else
    for p=k:-1:0
        r=kdeg(k-p,n-1);
        output(q:q+r-1,2:n)=rlex(k-p,n-1);
        output(q:q+r-1,1)=p;
        q=q+r;
    end
end
%-----------------------------------------------
function s=ldeg(A)
%計算多項式系統A的最高次數
Q=size(A,1);
n=size(A,2)-2;
s=0;
for i=1:Q
    sum=0;
    for j=1:n
        sum=sum+A(i,j);
    end
    if sum>s
        s=sum;
    end
end
%-----------------------------------------------
function c=cdeg(a,b)
%計算組合數
if a<=b
    c=1;
for e=1:b
    c=c*e;
end
for f=1:b-a
    c=c/f;
end
for g=1:a
    c=c/g;
end
else
    c=0;
end
%--------------------------------------------
function d=kdeg(k,n)
%計算n個未知數的k次單項式的總數
if k>0
    d=0;
    for i=1:k
        d=d+cdeg(i-1,k-1)*cdeg(i,n);
    end
else
    d=1;
end