# once we get candidates [u,v], before we calculate for their corresponding parameters [a,b1,b2,c], we divide them into several orbits. say the candidate [u,v] set is # C, then GetOrb(C) will give a set of orbits which involves all elements in C.

TR1 := proc(x,y) 
   if type(x,list) then return procname(op(x)) fi;
   normal( [x/(x-1), y/(y-1)] )
end:

TR2 := proc(x,y) 
   if type(x,list) then return procname(op(x)) fi;
   normal( [x/(x-1),(y-x)/(1-x)] )
end:

TR3 := proc(x,y) 
   if type(x,list) then return procname(op(x)) fi;
   normal( [(y-x)/(y-1),y/(y-1)] )
end:

TR4 := proc(x,y)
  if type(x,list) then return procname(op(x)) fi;
  normal(eval( [1-x,1-y] ))
end:

# orbit1 gives all elements with same orbit of (x,y) under 4 transformations.
orbit1:=proc(x,y)
  local GG, GG_old,g,gg;
  GG:={[x,y]}; GG_old:={};
  while GG<> GG_old and nops(GG)<5000 do
     GG_old:=GG;
     GG:=GG union normal({seq(TR1(i),i=GG), seq(TR2(i),i=GG), seq(TR3(i),i=GG),seq(TR4(i),i=GG)});
  od;
  normal(GG);
end:

# orbit2 will rule out elements whose degree is greater than 1.
orbit2:=proc(G)
  local GG, g, gg; GG:=G;
  for g in GG do
      for gg in g do
          if degree(denom(gg))>1 or degree(numer(gg))>1 or not has(gg,x) then GG:=GG minus {g}; fi;
      od;
  od;
  normal(GG);
end:

# GetOrb will divide all elements in "candiset" into several orbits.
GetOrb:=proc(candiset)
   local R, s, O, i,A;  A:={}; 
   R:=candiset;
   while R <> {} do
         s:= R[1];
         O:= orbit2(orbit1(s[1],s[2])); A:=A union {O intersect R};
         R:=R minus O;  
        
   od; A
end: