restart;
equA:=[x^2,y^2,z^2,x*y,x*z,y*z]:
with(Groebner):
with(LinearAlgebra):
with(PolynomialTools):

GBA:=Groebner[Basis](equA,plex(x,y,z)):
baseA, rvA:= NormalSet(GBA, plex(x,y,z));
dimA:=nops(baseA):

IsFlat:=proc(morph,eqB1) local phix, phiy, phiz, eqB, GBB, baseB, rvB, dimB, eqB_A, GBB_A, baseB_A, rvB_A, dimB_A, candidat, A,i,j,du,dv,dw;
  phix:=morph[1]:
  phiy:=morph[2]:
  phiz:=morph[3]:
  eqB:=[op(map2(subs,[x=phix,y=phiy,z=phiz],equA)),op(eqB1)]:
  GBB:=Groebner[Basis](eqB,plex(u,v,w)):
  baseB, rvB:= NormalSet(GBB, plex(u,v,w)):
  dimB:=nops(baseB):
  eqB_A:=[op(eqB),phix,phiy,phiz]:
  GBB_A:=Groebner[Basis](eqB_A,plex(u,v,w)):
  baseB_A, rvB_A:= NormalSet(GBB_A, plex(u,v,w)):
  dimB_A:=nops(baseB_A):
  if dimA*dimB_A<>dimB then return(false) fi:
  candidat:=seq(seq(op(j,baseA)*op(i,baseB_A),i=1..dimB_A),j=1..dimA):
  candidat:=map2(subs,[x=phix,y=phiy,z=phiz],[candidat]):
  A:=Matrix(dimB):
  for i from 1 to dimB do
    du[i]:=degree(baseB[i],u):
    dv[i]:=degree(baseB[i],v):
    dw[i]:=degree(baseB[i],w):
  od:
  for j from 1 to dimB do
    for i from 1 to dimB do
      A[i,j]:=coeff(coeff(coeff(candidat[j],u,du[i]),v,dv[i]),w,dw[i]):
    od:
  od:
  if Determinant(A)=0 then false else true fi
end;

C_Basis:=proc(morph,eqB1) local phix, phiy, phiz, eqB, GBB,baseB,rvB;
  phix:=morph[1]:
  phiy:=morph[2]:
  phiz:=morph[3]:
  eqB:=[op(map2(subs,[x=phix,y=phiy,z=phiz],equA)),op(eqB1)]:
  GBB:=Groebner[Basis](eqB,plex(u,v,w)):
  baseB, rvB:= NormalSet(GBB, plex(u,v,w)):
  baseB
end;

A_Basis:=proc(morph,eqB) local eqB_A,GBB_A,baseB_A,rvB_A;
  eqB_A:=[op(eqB),morph[1],morph[2],morph[3]]:
  GBB_A:=Groebner[Basis](eqB_A,plex(u,v,w)):
  baseB_A, rvB_A:= NormalSet(GBB_A, plex(u,v,w)):
  baseB_A
end;

listecoeff:=proc(pol,base,u,v,w) local i,c,l,du,dv,dw;
  du:=[]: dv:=[]: dw:=[]:
  for i from 1 to nops(base) do
    du:=[op(du),degree(base[i],u)]:
    dv:=[op(dv),degree(base[i],v)]:
    dw:=[op(dw),degree(base[i],w)]:
  od:
  l:=[]:
  for i from 1 to nops(base) do
    c:=coeff(coeff(coeff(pol,u,du[i]),v,dv[i]),w,dw[i]):
    if c<>0 then l:=[op(l),c] fi:
  od:
  l
end;

IsKind:=proc(morph,eqB) local phix, phiy, phiz, cas,delta,deltagen,f,d,a,b,c,i,j,equB,mbmaB,GBm,GBmd,basemd,rvmd,cond3,cond2,cond1,listevar,sol;
  phix:=morph[1]:
  phiy:=morph[2]:
  phiz:=morph[3]: 
  equB:=[op(map2(subs,[x=phix,y=phiy,z=phiz],equA)),op(eqB)]:
  mbmaB:=[op(equB),u*phix,u*phiy,u*phiz,v*phix,v*phiy,v*phiz,w*phix,w*phiy,w*phiz]:
  GBm:=Groebner[Basis](mbmaB,plex(u,v,w)):
  cas[1]:=[1,0,0]:
  cas[2]:=[0,1,0]:
  cas[3]:=[0,0,1]:
  cas[4]:=[1,1,0]:
  cas[5]:=[1,0,1]:
  cas[6]:=[0,1,1]:
  cas[7]:=[1,1,1]:
  deltagen:=d[1]*u+d[2]*v+d[3]*w:
  f:=a*phix+b*phiy+c*phiz:
  cond1:=[a,b,c]:
  for i from 1 to 7 do
    delta:=subs(d[1]=cas[i][1]*d[1],d[2]=cas[i][2]*d[2],d[3]=cas[i][3]*d[3],deltagen):
    GBmd:=Groebner[Basis]([op(GBm),delta],plex(u,v,w)):
    basemd, rvmd:=NormalSet(GBmd, plex(u,v,w)):
    cond3:=listecoeff(Reduce(f,GBmd,plex(u,v,w)),basemd,u,v,w):
    cond2:=[]:
    for j from 1 to 3 do
      if cas[i][j]=1 then cond2:=[op(cond2),d[j]] fi:      
    od:
    listevar:=[op(cond2),a,b,c]:
    for j from 1 to 3 do
      sol:=Solve(cond3,listevar,[op(cond2),cond1[j]]):
      if nops(sol)<>0 then return(true) fi:
    od:
  od:
  false
end;


#############------------Examples-------------------------

morph1:=[u^3,v^2,w^3]:
eqB1:=[]:
morph2:=[u^2+v*w,v^2+u*w,w^2+u*v]:
eqB2:=[]:
morph3:=[u^2,v^2,u^2+v^2]:
eqB3:=[w^4,u^2*w^2,v^2*w^2]:
morph4:=[u^3+v*w,v^2+u*w^2,w^3+u^2*v]:
eqB4:=[]:

IsFlat(morph1,eqB1);
IsFlat(morph2,eqB2);
IsFlat(morph3,eqB3);
IsFlat(morph4,eqB4);
C_Basis(morph1,eqB1);
C_Basis(morph2,eqB2);
C_Basis(morph3,eqB3);
C_Basis(morph4,eqB4);
A_Basis(morph1,eqB1);
A_Basis(morph2,eqB2);
A_Basis(morph3,eqB3);
A_Basis(morph4,eqB4);
IsKind(morph1,eqB1);
IsKind(morph2,eqB2);
IsKind(morph3,eqB3);
IsKind(morph4,eqB4);