###Maple code of the algorithm described in "On the problem of Detecting when Two Implicit Plane Algebraic Curves Are Similar", Juan Gerardo Alcázar, Gema Mª Diaz-Toca and Carlos Hermoso###

 

### the following procedure, semejanza, checks whether o not two implicit algebraic curves are related by a similarity. If they are, the procedure gives you 
###  -- $\omega$ and $r$ in the general case with cos\= 0
###  -- $\mu$ in the general case with cos= 0, and 
###  -- $b$ and $\bar{b}=beta$ in the special case. 


### input: C1:: an implicit curve ( C2 in the paper), C2:: the inverse image of C1 by the similarity (C1 in the paper).
	


semejanza:=proc(C1,C2) local F,F2,Frev,G,L,Fn,Gn,M,Deltaj,S,indicador,n,i,jbold,jj,arev,brev,lambda,lambda2,lambdarev,fn,gn,k,kk,kkrev,alphabold,betabold,alphaj,betaj,a,b,b2,Beta,Beta2,Betarev,Alpha,Alpharev,P,Prev,ffn,ggn,l;
 #1: general case;
 #3: special case;

if C1=C2 then return WARNING("The curves are equal");fi;

if not irreduc(C1) then return WARNING("The curve C1 is not irreducible");
 elif not irreduc(C2) then return WARNING("The curve C2 is not irreducible");
 elif degree(C1)<>degree(C2) then return WARNING("The curves do not have the same degree");
fi;
 


F:=sort(simplify(subs(x=(z+tau)/2,y=(z- tau)/(2*I),C1))); #tau :: \bar{z}
G:=sort(simplify(subs(x=(z+tau)/2,y=(z- tau)/(2*I),C2)));

n:=degree(C1);

L:=[op(F)];Fn:=0;for i in L do if degree(i)=n then Fn:=Fn+i;fi;od; 
L:=[op(G)];Gn:=0;for i in L do if degree(i)=n then Gn:=Gn+i;fi;od; 

coeffs(collect(Fn ,[z,tau],'distributed') ,[z,tau],'tt'); 
coeffs(collect(Gn ,[z,tau],'distributed') ,[z,tau],'ttt'); 
if ttt<>tt then return WARNING( "The curves are not similar because the powers of the coefficients of Fn and Gn are not same"); fi;

 
for i from 0 to n-1 do 
if coeff(Fn,z,n-i)<>0 then 
 jbold:=i;
 alphabold:= coeff(coeff(Fn,z,n-i),tau,i);
 betabold:=coeff(coeff(Gn,z,n-jbold),tau,jbold); break;
fi; 
od;

  
 
jj:=n+1:
for i from 0 to n-1 do 
 if coeff(coeff(Fn,z,n-i),tau,i)<>0 and 
   (n-i)^2*abs(coeff(coeff(Fn,z,n-i),tau,i))^2-(i+1)^2*abs(coeff(coeff(Fn, z,n-i-1 ),tau,i+1))^2<>0  
 then 
 jj:=i; alphaj:= coeff(coeff(Fn,z,n-i),tau,i);break;
 fi;
od;

k:=jj;

if k<n+1 then indicador:=1;  else  indicador:=3; fi;
 
if indicador=3 then  
  k:=1; 
  while coeff(coeff(G,z,n-1),tau,0)*coeff(coeff(G,tau,n-1),z,0)=0 do G:=simplify(subs(z=z+k,tau=tau+k,G));k:=k+1;od;  
  if k>1 then L:=[op(G)];Gn:=0;for i in L do if degree(i)=n then Gn:=Gn+i;fi;od; fi;
  k:=simplify( n*alphabold*b+coeff(coeff(Fn, z,n-1 ),tau,1)*Beta + coeff(coeff(F, z,n-1 ),tau,0) );
  kk:=1/(alphabold*coeff(coeff(G,z,n-1),tau,0));
  kkrev:=1/(alphabold*coeff(coeff(G,tau,n-1),z,0));
  a:=simplify ( expand(coeff(coeff(Gn,z,n),tau,0)*k*kk) );
  arev:=simplify ( expand(coeff(coeff(Gn,tau,n),z,0)*k*kkrev) );	
  lambda:=simplify(expand(alphabold*a^n/coeff(coeff(Gn,z,n),tau,0)));
  lambdarev:=simplify(expand(alphabold*arev^n/coeff(coeff(Gn,tau,n),z,0)));
  k:=simplify( n*conjugate(alphabold)*Beta+conjugate(coeff(coeff(Fn, z,n-1 ),tau,1))*b + conjugate(coeff(coeff(F, z,n-1 ),tau,0)) );
  kk:=1/(conjugate(alphabold)*conjugate(coeff(coeff(G,z,n-1),tau,0)));
  kkrev:=1/(conjugate(alphabold)*conjugate(coeff(coeff(G,tau,n-1),z,0)));
  Alpha:=simplify (expand( conjugate(coeff(coeff(Gn,z,n),tau,0))*k*kk ));
  Alpharev:=simplify (expand( conjugate(coeff(coeff(Gn,tau,n),z,0))*k*kkrev ));
  F2:=simplify(expand(subs( z =a*z+b  ,tau=Alpha*tau+Beta,F))); #esto es un pol en z y tau,b, BETA, eq (2)
  S:=[coeffs(collect(F2 - simplify(expand(lambda*G)),[z,tau],'distributed') ,[z,tau])];
   
  WARNING("special case, or. preserving");
  #print(S);
  print(SolveTools[PolynomialSystem](S,[b,Beta],engine= traditional));
  printf("\n a in b,Beta is %a \n", a);
  printf("\n lambda is %a \n", lambda);

 
  Frev:= subs(z=Z,tau=Tau,F) ; 
  F2:=simplify(expand(subs( Z =arev*tau+b ,Tau=Alpharev*z+Beta,Frev)));  
  S:=[coeffs(collect(F2 - simplify(expand(lambdarev*G)),[z,tau],'distributed') ,[z,tau])];
    
  WARNING("special case, or. reversing "); 
 
  print(SolveTools[PolynomialSystem](S,[b,Beta],engine= traditional));
  printf("\n a in b,Beta is %a \n", arev);
  printf("\n lambda is %a \n", lambdarev);

fi:

 
if  indicador=1 then  
  lambda:=simplify(alphaj*a^(n-jj)*Alpha^jj/coeff(coeff(Gn,z,n-jj),tau,jj));
 
  Deltaj:=simplify((n-jj)^2*abs(alphaj)^2-(jj+1)^2*abs(coeff(coeff(Fn, z,n-jj-1 ),tau,jj+1))^2);

 

  M:=Matrix(2,2);

 
  M[1,1]:=simplify(a*alphaj*coeff(coeff(G, z,n-jj-1 ),tau,jj)/coeff(coeff(Gn, z,n-jj ),tau,jj) -coeff(coeff(F, z,n-jj-1 ),tau,jj));
 
  M[2,1]:=simplify(Alpha*conjugate(alphaj*coeff(coeff(G, z,n-jj-1 ),tau,jj))/conjugate(coeff(coeff(Gn, z,n-jj ),tau,jj)) -conjugate(coeff(coeff(F, z,n-jj-1 ),tau,jj)));

  M[1,2]:=(jj+1)*coeff(coeff(Fn, z,n-jj-1 ),tau,jj+1);

  M[2,2]:=(n-jj)*conjugate(alphaj);
 
  b:=simplify(LinearAlgebra[Determinant](M)/Deltaj);

 

  M:=Matrix(2,2);

  M[1,2]:=simplify(a*alphaj*coeff(coeff(G, z,n-jj-1 ),tau,jj)/coeff(coeff(Gn, z,n-jj ),tau,jj) -coeff(coeff(F, z,n-jj-1 ),tau,jj));
 
  M[2,2]:=simplify(Alpha*conjugate(alphaj*coeff(coeff(G, z,n-jj-1 ),tau,jj))/conjugate(coeff(coeff(Gn, z,n-jj ),tau,jj)) -conjugate(coeff(coeff(F, z,n-jj-1 ),tau,jj)));

  M[2,1]:=(jj+1)*conjugate(coeff(coeff(Fn, z,n-jj-1 ),tau,jj+1));

  M[1,1]:=(n-jj)*alphaj;

  Beta:=simplify(LinearAlgebra[Determinant](M)/Deltaj);
 
   
  
  lambdarev:=simplify(alphaj*a^(n-jj)*Alpha^jj/coeff(coeff(Gn,tau,n-jj),z,jj));
   
  M:=Matrix(2,2);

  M[1,1]:=simplify(a*alphaj*coeff(coeff(G, tau,n-jj-1 ),z,jj)/coeff(coeff(Gn, tau,n-jj ),z,jj) -coeff(coeff(F, z,n-jj-1 ),tau,jj));
 
  M[2,1]:=simplify(Alpha*conjugate(alphaj*coeff(coeff(G, tau,n-jj-1 ),z,jj))/conjugate(coeff(coeff(Gn, tau,n-jj ),z,jj)) -conjugate(coeff(coeff(F, z,n-jj-1 ),tau,jj)));

  M[1,2]:=(jj+1)*coeff(coeff(Fn, z,n-jj-1 ),tau,jj+1);

  M[2,2]:=(n-jj)*conjugate(alphaj);
 
  brev:=simplify(LinearAlgebra[Determinant](M)/Deltaj);


  M:=Matrix(2,2);

  M[1,2]:=simplify(a*alphaj*coeff(coeff(G, tau,n-jj-1 ),z,jj)/coeff(coeff(Gn, tau,n-jj ),z,jj) -coeff(coeff(F, z,n-jj-1 ),tau,jj));
 
  M[2,2]:=simplify(Alpha*conjugate(alphaj*coeff(coeff(G, tau,n-jj-1 ),z,jj))/conjugate(coeff(coeff(Gn, tau,n-jj ),z,jj)) -conjugate(coeff(coeff(F, z,n-jj-1 ),tau,jj)));

  M[2,1]:=(jj+1)*conjugate(coeff(coeff(Fn, z,n-jj-1 ),tau,jj+1));

  M[1,1]:=(n-jj)*alphaj;

  Betarev:=simplify(LinearAlgebra[Determinant](M)/Deltaj);
 
  L:=[op(C1)];fn:=0;for i in L do if degree(i)=n then fn:=fn+i;fi;od;fn;
  L:=[op(C2)];gn:=0;for i in L do if degree(i)=n then gn:=gn+i;fi;od;fn;

  L:=factors(fn);
 

  
 
  if subs(x=0,fn)=0 and subs(x=0,gn)=0 and nops(L[2])=2 then 
   ffn:=quo(fn,gcd(fn,x^n),x);k:=degree(gcd(fn,x^n));
   L:=factors(ffn)[2][1][1];
   if degree(L)=1 then  
   	ggn:=quo(gn,gcd(gn,x^n),x);l:=degree(gcd(gn,x^n));  
      if k<>l and k<>n-l then return  WARNING("case 2 fails, curves are not similar"); fi;
      L:=factors(ggn)[2][1][1];
      if degree(L)=1 then   
  		if subs(y=0,fn)=0 and subs(y=0,gn)=0 then 
			if k<>n-k then 
				if k=l then P:=t; fi;
			else
				P:=t;  
				Prev:=t;
			fi;
		else
			 
			if k<>l then P:=t+coeff(L,y)/coeff(L,x); Prev:=t-coeff(L,y)/coeff(L,x);   
			elif  (n-k)<>k then P:=t;Prev:=t;
			else P:=t*(t+coeff(L,y)/coeff(L,x));  
			     
			     Prev:=t*(t-coeff(L,y)/coeff(L,x));  
			fi;
		fi;
	fi;
   fi;

  else 
   
   L:=subs(x=1,gn);  
   L:=numer(normal(algsubs(y=(y-t)/(1+y*t),L))); #by experiments...
   P:=resultant(subs(x=1,fn), L,y);
   L:=subs(x=1,gn);
   L:=numer(normal(algsubs(y=(t-y)/(1+y*t),L)));
   Prev:=resultant(subs(x=1,fn), L,y);

  fi;


 

 b2:=subs(a=r*(1+I*omega),Alpha=r*(1-I*omega),b);
 lambda2:=subs(a=r*(1+I*omega),Alpha=r*(1-I*omega),lambda);
 Beta2:=subs(a=r*(1+I*omega),Alpha=r*(1-I*omega),Beta);
 F2:=simplify(subs( z =r*(1+I*omega)*z+b2 ,tau=r*(1-I*omega)*tau+Beta2,F));  
  
 S:=[coeffs(collect(F2 - lambda2*G,[z,tau],'distributed') ,[z,tau]),subs(t=omega,P)]; 
 
 WARNING("general case, cos<>0, or.preserving:");
 
 print(SolveTools[PolynomialSystem](S,[r,omega],engine= traditional));
 printf("\n b in r, omega is %a \n", b2);
 printf("\n lambda is %a \n", lambda2);

 
 b2:=subs(a=I*mu,Alpha=-I*mu,b);
 lambda2:=subs(a=I*mu,Alpha=-I*mu,lambda); 
 Beta2:=subs(a=I*mu,Alpha=-I*mu,Beta);
 F2:=simplify(subs( z =I*mu*z+b2 ,tau=-I*mu*tau+Beta2,F));  
 S:=[coeffs(collect(F2 - lambda2*G,[z,tau],'distributed') ,[z,tau])]; 
 WARNING("general case, cos=0, or.preserving:");
 #print(S);
 print(SolveTools[PolynomialSystem](S,[mu],engine= traditional));
 printf("\n b in mu is %a \n", b2);


 b2:=subs(a=r*(1+I*omega),Alpha=r*(1-I*omega),brev);
 lambda2:=subs(a=r*(1+I*omega),Alpha=r*(1-I*omega),lambdarev);
 Beta2:=subs(a=r*(1+I*omega),Alpha=r*(1-I*omega),Betarev);
 Frev:= subs(z=Z,tau=Tau,F) ;#tau conjugado de z
 F2:=simplify(subs( Z =r*(1+I*omega)*tau+b2 ,Tau=r*(1-I*omega)*z+Beta2,Frev));  
 S:=[coeffs(collect( F2 - lambda2*G,[z,tau],'distributed') ,[z,tau]),simplify(subs(t=omega,  Prev) )]; 
 
 WARNING("general case, cos<>0, or.reversing:");
  
 print(SolveTools[PolynomialSystem](S,[r,omega],engine= traditional));
 printf("\n b in r, omega is %a \n", b2);
 printf("\n lambda is %a \n", lambda2);
 
 

 b2:=subs(a=I*mu,Alpha=-I*mu,brev);
 lambda2:=subs(a=I*mu,Alpha=-I*mu,lambdarev);
 Beta2:=subs(a=I*mu,Alpha=-I*mu,Betarev);
 Frev:= subs(z=Z,tau=Tau,F) ;
 F2:=simplify(subs( Z =I*mu*tau+b2 ,Tau=-I*mu*z+Beta2,Frev));  
 S:=[coeffs(collect(F2 - lambda2*G,[z,tau],'distributed') ,[z,tau])]; 
 
 WARNING("general case, cos=0, or.reversing:");
 
 print(SolveTools[PolynomialSystem](S,[mu],engine= traditional));
 printf("\n b in mu is %a \n", b2)

fi;

end:

######in order to obtain similar curves

pre_curve_ztau:=proc(F,a,b) local k;
k:=sort(simplify(subs(x=(z+tau)/2,y=(z- tau)/(2*I),F)));
return simplify(subs(z=a*z+b,tau=conjugate(a)*tau+conjugate(b),k));
end;

rev_curve_ztau:=proc(F,a,b) local k,kk;
k:=sort(simplify(subs(x=(z+tau)/2,y=(z- tau)/(2*I),F )));
kk:=subs(z=Z,tau=T,k);
return simplify(subs(Z=a*tau+b,T=conjugate(a)*z+conjugate(b),kk));
end;