#debugging function			
  Assert:= function(b,s)
    if b then 
      return;
    else return Error("Assumption: ",s,"failed!");
    fi;
  end;

#some functions that are forwarded
   MultMap:=function(a)
     return a;
   end;

   IsogenyFromSubgroup:=function(a)
      return a;
   end;

   extFrob:=function(a)
      return a;
   end;

   IsogenyFromKernel:=function(a)
       return a;
   end;

   IdealFactorization:=function(a)
      return a;
    end;

   MakeSmooth2:=function(a)
      return a;
   end;

   PicardNumber:=function(a)
      return a;
   end;

   PseudoRandom:=function(a)
      return a;
   end;

   FindIso:=function(a)
      return a;
   end;

#true if the isogeny f is the zero-map
  IsZero:=function(f)
    return "zero" in RecFields(f);
  end;


#isogeny-operations
  _iso_ops:=rec();
  _iso_ops.Print:= function(f)
    if IsZero(f) then 
      if f.domain = f.codomain then
        Print(" ZeroMap on ",f.domain,"\n");
      else Print("ZeroMap from ",f.domain, "to ",f.codomain,"\n"); 
      fi;
    else
      if "fro" in RecFields(f) then 
        Print("Isogeny from ",f.domain," to ",f.codomain,"(x,y)|->(x^",Degree(f),", y^",Degree(f),")\n");
        return;
      fi;
      AssignNames_(CoefficientRing(Parent(f.xpol)),["x"]);
      AssignNames_(Parent(f.ypol),["y"]);
      Print("Isogeny from ",f.domain," to ",f.codomain,"(x,y)|->(",f.xpol,",");
      Print(f.ypol,")\n");
    fi;
  end;


#Frobenius-isogeny of degree q
   FrobeniusToOtherCurve:=function(E,q)
     local E2,j,L,BR,L2,i,_frob,frob,x,y,j2,Funf;
     j:=jInvariant(E);
     L:=Coefficients(E);
     BR:=BaseRing(E);
     x:=PolynomialAlgebra(BR).1;
     Funf:=EllipticFunctionField(E);
     y:=Funf.1;
     if GCD(q,Characteristic(BR)) = 1 then
         return Error;
     fi;
     L2:=[];
     for i in [1..Length(L)] do
         L2[i]:=L[i]^q;
     od;
     E2:=EllipticCurve(BR,L2);
     _frob:=function(P)
        if P = Zero(Parent(P)) then
           return Zero(E2);
        else return MakePoint(E2,[PX(P)^q, PY(P)^q]);
        fi;
     end;
     frob:=Map(E,E2,_frob);
     frob.type:= iso;
     frob.fro:=true;
     frob.operations:= _iso_ops;
     frob.degree:=q;
     frob.xpol:=Coerce(Funf,x^q);
     frob.ypol:=y^q;
     return frob;
   end;


# computes invariants a_i and b_i
   ECInvariants:=function(E)
     local L,Inv;
     L:=Coefficients(E);
     Inv:=rec();
     if Length(L) = 5 then
        Inv.a1:=L[1];
        Inv.a3:=L[2];
        Inv.a2:=L[3];
        Inv.a4:=L[4];
        Inv.a6:=L[5];
     else
        Inv.a1:=0;
        Inv.a3:=0;
        Inv.a2:=0;
        Inv.a4:=L[1];
        Inv.a6:=L[2];
     fi;
     Inv.b2:= Inv.a1^2+4*Inv.a2;
     Inv.b4:=2*Inv.a4+Inv.a1*Inv.a3;
     Inv.b6:=Inv.a3^2+4*Inv.a6;
     Inv.b8:=Inv.a1^2*Inv.a6+4*Inv.a2*Inv.a6-Inv.a1*Inv.a3*Inv.a4+Inv.a2*Inv.a3^2-Inv.a4^2;
     return Inv;
   end;


# the generic point (x,y) of a curve
   GenericPoint:= function(E)
     local i,Funf,y0,x0,EE,x,y,coLi;
     Funf:=_OldFunf(E);
     coLi:=[];
     for i in [1..Length(E.L)] do
        coLi[i]:=E.L[i];
     od;
     EE:= EllipticCurve(Funf,coLi);
     y:= BaseRing(EE).1;
     x:= BaseField(BaseRing(EE)).1;
     return MakePoint(EE,[x,y]);
   end;


# computes an isomorphic curve in short WeierstrassNormalForm
   ShortWNF:=function(E)
     local L,Inv,E2,Iso,xmap,ymap,P,PA,PA2,x,y,f,ff,co;
     L:=Coefficients(E);
     if Length(L) = 2 then
        return rec(base:=E,ext1:=Isomorphism(E,E,1));
     else Inv:=ECInvariants(E);
        E2:=_WeierstrassNormalForm(BaseRing(E),L);
        E2:=EllipticCurve(BaseRing(E),E2[2]);
        xmap:=function(P)
          return 36*PX(P)+3*Inv.b2;
        end;
        ymap:=function(P)
          return (2*PY(P)+Inv.a1*PX(P)+Inv.a3)*108;
        end;
        Iso:=Map(E,E2,
           function(P) 
             if P = Zero(E) then 
               return Zero(E2); 
             fi; 
             return Coerce(E2,[xmap(P),ymap(P)]);
           end 
         );
        P:=GenericPoint(E);
        Iso.type:= iso;
        Iso.operations:= _iso_ops;
        Iso.xpol:= xmap(P);
        Iso.ypol:= ymap(P);
        Iso.degree:= 1;
        return rec(base:=E2,ext1:=Iso);
     fi;
   end;


# returns -f if f is an isogeny
   Negate:=function(f)
     local fmap,Iso,E2,a1,a3,L;
     fmap:=function(P)
       return -f(P);
     end;
     E2:=Codomain(f);
     L:=Coefficients(E2);
     if Length(L) = 2 then
       a1:=0;
       a3:=0;
     else 
       a1:=L[1];
       a3:=L[2];
     fi;
     Iso:=Map(Domain(f),E2,fmap);
     Iso.type:=iso;
     Iso.operations:=_iso_ops;
     Iso.xpol:=f.xpol;
     Iso.ypol:=f.ypol-a1*f.xpol-a3;
     Iso.degree:=f.degree;
     return Iso;
   end;



# returns true if E is of the form y^2+xy=x^3+bx^2+a over GF(2,n) and trace(b) =1 or b =0

   InNormalFormChar2:=function(E) 
     local BR,L,b1,b2,b;
     BR:=BaseRing(E);
     L:=Coefficients(E);
     if not Characteristic(BR) =2 then
        return Error("function only applicable for curves over fields of genus 2");
     fi;
     b1:= (L[1]=1 and L[2] = 0 and L[4] =0);
     b:=L[3];
     if Is(Type(Parent(b)),fld^fun) then
        b2:= (Trace(Coerce(CoefficientRing(CoefficientRing(Parent(b))),b)) <> 0 ) or (L[3] = 0);
     else 
        b2:=(Trace(b) <>0) or (b =0);
     fi;
     return b1 and b2;
   end;


#isomorphic curve of the form y^2+xy=x^3+bx^2+a over GF(2,r), n defines direction of isomorphism

   NormalFormChar2:=function(E,n) 
     local L,f,BR,x,a1,a2,a3,a4,a6,g,b,a,s,E2,Iso,_Iso,xmap,ymap,P,EE,b2;
     if InNormalFormChar2(E) then
         return rec(base:=E,ext1:=Isomorphism(E,E,1));
     fi;
     if n <> 1 and n <> 2 then 
        return Error(" n must be in [1,2] ");
     fi;
     L:=Coefficients(E);
     BR:=BaseRing(E);
     f:=_GetPoly(BR,L);
     x:=Parent(f).1;
     a1:=L[1];
     a3:=L[2];
     a2:=L[3];
     a4:=L[4];
     a6:=L[5];
     b:=((a1^3)*a3+(a1^4)*a2)/(a1^6);
     s:=0;
     if Trace(b) = 0 then
        s:=Roots(PolynomialAlgebra(BR).1^2+PolynomialAlgebra(BR).1+b)[1][1];
        b:=0;
     fi;
     a:=(a1^4*a4^2+a3^4+a1^3*a3^3+a1^4*a2*a3^2+a1^5*a3*a4+a1^6*a6)/(a1^12);
     E2:=EllipticCurve(BR,[1,0,b,0,a]);
     E2.size:=E.size;
     if n = 1 then 
        P:=GenericPoint(E);
        xmap:= function(P)
          if PX(P) = INFTY then 
            return INFTY; 
          fi; 
          return (PX(P)+a3/a1)/a1^2;
        end;
        ymap:= function(P) 
          if PY(P) = INFTY then 
            return INFTY; 
          fi;
          return (PY(P)+(a1^2*a4+a3^2)/a1^3+s*xmap(P))/a1^3;
        end;
        _Iso:= function(P)
            if P = Zero(E) then
              return Zero(E2);
           else return Coerce(E2,[xmap(P),ymap(P)]);
           fi;
        end;
        Iso:=Map(E,E2,_Iso);
     else 
        P:=GenericPoint(E2);
        xmap:= function(P)
          if PX(P) = INFTY then 
              return INFTY; 
          fi; 
          return a1^2*PX(P)+a3/a1;
        end;
        ymap:= function(P) 
          if PY(P) = INFTY then 
              return INFTY; 
          fi;
          return a1^3*PY(P)+(a1^2*a4+a3^2)/a1^3+s*xmap(P);
        end;
        _Iso:= function(P)
           if P = Zero(E) then
              return Zero(E2);
           else return Coerce(E,[xmap(P),ymap(P)]);
           fi;
        end;
        Iso:=Map(E2,E,_Iso);
     fi;
     Iso.type:= iso;
     Iso.operations:= _iso_ops;
     Iso.xpol:= xmap(P);
     Iso.ypol:= ymap(P);
     Iso.degree:= Degree(Iso.xpol)/2;
     AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]);
     AssignNames_(Parent(Iso.ypol),["y"]);
     return rec(base:=E2,ext1:=Iso);
   end;


#recursion for division polynomials
   _ff := function(EC,m) 
     local a1,a2,a3,a4,a6,L,b2,b4,b6,b8,Phi,kx,F,x,y,n;
     L:= Coefficients(EC);
     if Length(L)=2 then
        a1:=0;
        a3:=0;
        a2:=0;
        a4:=L[1];
        a6:=L[2];
     else 
        a1:=L[1];
        a3:=L[2];
        a2:=L[3];
        a4:=L[4];
        a6:=L[5];
     fi;
     b2:=a1^2+4*a2;
     b4:=a1*a3+2*a4;
     b6:=a3^2+4*a6;
     b8:=a1^2*a6+4*a2*a6-a1*a3*a4+a2*a3^2-a4^2;
     kx:= PolynomialAlgebra(BaseRing(EC));
     x:= kx.1;
     F:= EllipticFunctionField(EC);
     y:=F.1;
     if m<=0 then return 0;
     elif m =1 then return 1;
     elif m=2 then return 2*y+a1*x+a3;
     elif m=3 then return 3*x^4+b2*x^3+3*b4*x^2+3*b6*x+b8;
     elif m=4 then return (2*x^6+b2*x^5+5*b4*x^4+10*b6*x^3+10*b8*x^2+(b2*b8-b4*b6)*x+b4*b8-b6^2)*(2*y+a1*x+a3);
     elif IsOdd(Coerce(Integers(),m)) then 
        n:=(m-1)/2;
        return _ff(EC,n+2)*_ff(EC,n)^3-_ff(EC,n-1)*_ff(EC,n+1)^3;
     else 
        n:=m/2;
        return  (_ff(EC,n+2)*_ff(EC,n-1)^2-_ff(EC,n-2)*_ff(EC,n+1)^2)*_ff(EC,n)/(2*y+a1*x+a3);
     fi;
   end;


#recursion for division polys over char 2
   _ff2 := function(EC,m) 
     local a1,a2,a3,a4,a6,L,b2,b4,b6,b8,Phi,kx,F,x,y,n,E,Iso,li;
     L:= Coefficients(EC);
     a2:=L[3];
     a6:=L[5];
     x:=PolynomialAlgebra(BaseRing(EC)).1;
     F:= EllipticFunctionField(EC);
     y:=F.1;
     if m<=0 then return 0;
     elif m =1 then return 1;
     elif m=2 then return x;
     elif m=3 then return 3*x^4+x^3+a6;
     elif m=4 then return x^6+a6*x^2;
     elif IsOdd(Coerce(Integers(),m)) then 
        n:=(m-1)/2;
        return _ff2(EC,n+2)*_ff2(EC,n)^3-_ff2(EC,n-1)*_ff2(EC,n+1)^3;
     else 
        n:=m/2;
        return (_ff2(EC,n+2)*_ff2(EC,n-1)^2-_ff2(EC,n-2)*_ff2(EC,n+1)^2)*_ff2(EC,n)/x;
     fi;
   end;

#Division polynomials for characteristic 2
   _DivisionPolynomial2:=function(EC,m) 
      local PA,L,a2,a6,Phi,Phi2,els,ro,kx,x,F,y,E,Iso;
      L:= Coefficients(EC);
      PA:= PolynomialAlgebra(BaseRing(EC));
      x:= PA.1;
      return Coerce(PA,_ff2(EC,m));
   end;

#returns Division polynomial of degree m
   DivisionPolynomial:=function(EC,m) 
     local PA,L,a1,a3,Phi,Phi2,els,ro,kx,x,F,y;
     if Characteristic(BaseRing(EC)) =2 then  
        return  _DivisionPolynomial2(EC,m); 
     fi;
     L:= Coefficients(EC);
     if Length(L) <>2 then 
        return Error("division polynomial: expecting elliptic curve in short WNF");
     fi;
     a1:=L[1];
     a3:=L[2];
     kx:= PolynomialAlgebra(BaseRing(EC));
     x:= kx.1;
     F:= EllipticFunctionField(EC);
     y:=F.1;
     PA:= PolynomialAlgebra(BaseRing(EC));
     if m mod 2 = 0 then
        Phi:= Coerce(PA,_ff(EC,m)/_ff(EC,2));
        if Characteristic(BaseRing(EC))>3 then
           if Length(L) = 2 then
              return Phi*(x^3+a1*x+a3);
           else return Error("case not implemented yet");
           fi;
        elif Characteristic(BaseRing(EC))=3 then return Error("case not implemented yet");
        else Phi2:= -a1*x-a3;
        fi;
        return Phi*Phi2;
     elif IsOdd(m) then 
        return Coerce(PA,_ff(EC,m));
     else 
        return Coerce(PA,_ff(EC,m)/_ff(EC,2));
     fi;
   end;




#the characteristic polynomial of the frobenius endomorphism
   CharacteristicFrobeniusPolynomial:= function(E)
     local t,size,p,Zx;
     size:= Size(E);
     p:= Size(BaseRing(E));
     Zx:= PolynomialAlgebra(Integers());
     t:= p+1-size;
     return Coerce(Zx,[p,-t,1]);
   end;


#returns true if l is an Elkies prime of EC
   IsElkiesPrime:=function(EC,l)
     local mo,F,PA,ro,nr,i,ch,O;
     ch:= CharacteristicFrobeniusPolynomial(EC);
     O:=EquationOrder(ch);
     if Conductor(O) mod l = 0 then 
        return rec(base:=false); 
     fi;
     PA:= PolynomialAlgebra(GF(l));
     ch:= Coerce(PA,ch);
     ro:=Roots(ch);
     return rec(base:=ro <>Sequence([]),ext1:=ro);
   end;


#returns Elkies primes <=37 of E 
   ElkiesPrimes:= function(E)
     local i,el,bool;
     el:= [];
     for i in [2..Minimum(Size(BaseRing(E))-1,37)] do
        if IsPrime(i) then
           bool:=IsElkiesPrime(E,i);
           if  bool.base then
               el:= Append(el,[i]);
           fi;
        fi;
     od;
     return el;
   end;

#returns Elkies primes <=59 of E 
   ElkiesPrimes59:= function(E)
     local i,el,elli,bool;
     el:= [];
     elli:=[];
     for i in [2..Minimum(Size(BaseRing(E))-1,59)] do
        if IsPrime(i) then
           bool:=IsElkiesPrime(E,i);
           if  bool.base then
              el:= Append(el,[i]);
           fi;
        fi;
     od;
     return el;
   end;




#the torsion subgroup EC[l]/K
   TorsionSubgroup:=function(EC,l)
     local di,i,ro,li,litmp;
     di:= DivisionPolynomial(EC,l);
     ro:= Roots(di);
     if ro = [] then
        return [Zero(EC)];
     else li:=[];
        for i in [1..Length(ro)] do
            litmp:= MakePoints(EC,ro[i][1]);
            Append_(li,litmp);
        od;
        Append_(li,[Zero(EC)]);
        return li;
     fi;
   end;

#recursion for evaluation of Division polynomial at certain point in char 2
   _DivisionPolynomialEval2 := function(EC,m,P)  
     local a1,a2,a3,a4,a6,L,b2,b4,b6,b8,Phi,kx,F,x,y,n,E,Iso1,nor,Iso2,re,P2;
     nor:=InNormalFormChar2(EC);
     if not nor then
        E:=NormalFormChar2(EC,1);
        Iso1:=E.ext1;
        Iso2:=NormalFormChar2(EC,2).ext1;
        E:=Base(E);
     else 
        Iso1:=false;
        E:=EC;
     fi;
     L:= Coefficients(E);
     a2:=L[3];
     a6:=L[5];
     if Iso1 = false then
        x:=PX(P);
        y:=PY(P);
     else 
        x:=PX(Iso1(P));
        y:=PY(Iso1(P));
     fi;
     if m<=0 then return 0;
     elif m =1 then 
        re:=1;
        if Iso1 = false then
           return re;
        else 
           P2:=MakePoints(E,re)[1];
           return PX(Iso2(P2));
        fi;
     elif m=2 then 
        re:=x;
        if Iso1 = false then
           return re;
        else 
           P2:=MakePoints(E,re)[1];
           return PX(Iso2(P2));
        fi;
     elif m=3 then re:= x^4+x^3+a6;
        if Iso1 = false then
           return re;
        else 
           P2:=MakePoints(E,re)[1];
           return PX(Iso2(P2));
        fi;
     elif m=4 then re:= x^6+a6*x^2;
        if Iso1 = false then
           return re;
        else 
           P2:=MakePoints(E,re)[1];
           return PX(Iso2(P2));
        fi;
     elif IsOdd(Coerce(Integers(),m)) then 
        n:=(m-1)/2;
        re:= _DivisionPolynomialEval2(EC,n+2,P)*_DivisionPolynomialEval2(EC,n,P)^3+_DivisionPolynomialEval2(EC,n-1,P)*_DivisionPolynomialEval2(EC,n+1,P)^3;
        if Iso1 = false then
           return re;
        else 
           P2:=MakePoints(E,re)[1];
           return PX(Iso2(P2));
        fi;
     else 
        n:=m/2;
        re:= (_DivisionPolynomialEval2(EC,n+2,P)*_DivisionPolynomialEval2(EC,n-1,P)^2+_DivisionPolynomialEval2(EC,n-2,P)*_DivisionPolynomialEval2(EC,n+1,P)^2)*_DivisionPolynomialEval2(EC,n,P)/x;
        if Iso1 = false then
           return re;
        else 
           P2:=MakePoints(E,re)[1];
           return PX(Iso2(P2));
        fi;
     fi;
   end;




#recursion for evaluation of Division polynomial at certain point
   _DivisionPolynomialEval := function(EC,m,P) 
     local a1,a2,a3,a4,a6,L,b2,b4,b6,b8,Phi,kx,F,x,y,n;
     if Characteristic(BaseRing(EC))=2 then 
        return _DivisionPolynomialEval2(EC,m,P); 
     fi;
     L:= Coefficients(EC);
     if Length(L)=2 then
        a1:=0;
        a3:=0;
        a2:=0;
        a4:=L[1];
        a6:=L[2];
     else 
        a1:=L[1];
        a3:=L[2];
        a2:=L[3];
        a4:=L[4];
        a6:=L[5];
     fi;
     b2:=a1^2+4*a2;
     b4:=a1*a3+2*a4;
     b6:=a3^2+4*a6;
     b8:=a1^2*a6+4*a2*a6-a1*a3*a4+a2*a3^2-a4^2;
     kx:= PolynomialAlgebra(BaseRing(EC));
     x:= PX(P);
     F:= EllipticFunctionField(EC);
     y:= PY(P);
     if m<=0 then return 0;
     elif m =1 then return 1;
     elif m=2 then return 2*y+a1*x+a3;
     elif m=3 then return 3*x^4+b2*x^3+3*b4*x^2+3*b6*x+b8;
     elif m=4 then return (2*x^6+b2*x^5+5*b4*x^4+10*b6*x^3+10*b8*x^2+(b2*b8-b4*b6)*x+b4*b8-b6^2)*(2*y+a1*x+a3);
     elif IsOdd(Coerce(Integers(),m)) then 
         n:=(m-1)/2;
         return _DivisionPolynomialEval(EC,n+2,P)*_DivisionPolynomialEval(EC,n,P)^3-_DivisionPolynomialEval(EC,n-1,P)*_DivisionPolynomialEval(EC,n+1,P)^3;
     else 
         n:=m/2;
         return (_DivisionPolynomialEval(EC,n+2,P)*_DivisionPolynomialEval(EC,n-1,P)^2-_DivisionPolynomialEval(EC,n-2,P)*_DivisionPolynomialEval(EC,n+1,P)^2)*_DivisionPolynomialEval(EC,n,P)/(2*y+a1*x+a3);
     fi;
   end;






#Conductor of an order
   _IntCond:= function(O)
     local Omax;
     Omax:= MaximalOrder(O);
     return Index(Omax,O);
   end;


#the i-th power of the frobenius endomorphism of E
   Frobenius:= function(E,i)
     local ma,p,frob,EE,P,q;
     P:= GenericPoint(E);
     EE:= Parent(P);
     p:= Size(BaseRing(E));
     q:= p^i;
     frob:= function(P1)
       local P2,E2;
       E2:= Parent(P1);
       return MakePoint(E2,[PX(P1)^q, PY(P1)^q]);
     end;
     ma:= Map(E,E,frob);
     ma.type:= iso;
     ma.fro:=true;
     ma.operations:= _iso_ops;
     ma.degree:=q;
     return ma;
   end;


#isogeny of degree 0, i.e. phi(P) = O for all P in E1(K)
   NuMap:=function(E1,E2) 
     local ma,ze;
     ze:=function(P)
       return Zero(E2);
     end;
     ma:=Map(E1,E2,ze);
     ma.type:=iso;
     ma.operations:= _iso_ops;
     ma.degree:=0;
     ma.zero:=true;
     return ma;
   end;

#converts integer to binary list
   _IntToBinList:=function(n) 
     local L;
     L:=[];
     while n>0 do
       Append_(L,[n mod 2]);
       n:= Div(n,2);
     od;
     return L;
   end;
#converts integer to finite field element, T generator of the field, only for characteristic 2
   _IntToFldElt:=function(n,T) 
     local T,K,L,i,el;
     L:=_IntToBinList(n);
     el:=0;
     for i in [1..Length(L)] do
       el:=el+L[i]*T^(i-1);
     od;
     return el;
   end;



#returns the dual isogeny of f
   DualIsogeny:= function(f)  
     local VR,ba,li,M,apoly,bpoly,polys,i,d,isodu,xpoly,xmap,E2,Fu,tmp,x,fu,mux, xpol,ypol,ymap,ypoly,y,yoben,yunten,BR,E,ymupol,Funf,xpol,muma,xmuma,xoben,xunten,PA,xmupol,Iso,f,dy,fact,P,j,dego1,dego2,degu1,degu2,GenP,_Iso,d2,xpart,ypart,xdenom,ydenom,yy,fact,ro,E3,psi,isotest;
     E:= f.domain;
     E2:= f.codomain;
     BR:= BaseRing(E);
     PA:= PolynomialAlgebra(BR);
     polys:= [];
     xpol:= f.xpol;
     d:= Coerce(Integers(),Degree(f));
     if not IsPrime(d) then
        return Error("assuming prime degree for dual isogeny-computation");
     fi;
     muma:=MultMap(E,d);
     for i in [1..d] do
        polys[i]:= xpol^i;
     od;
     M:= [xpol^0];
     for i in [2..d+1] do
        M[i]:= polys[i-1];
     od;
     mux:=-Coerce(Parent(xpol),muma.xpol);
     M:= Append(M,[mux]);
     for i in [2..d] do
        M:= Append(M,[mux*polys[i-1]]);
     od;
     VR:= Relations(M,BR,1);
     ba:= Basis(VR);
     Assert(ba<> Sequence([])," Vectorspace of dimension > 0 ");
     li:= List(ba[1]);
     apoly:= [];
     bpoly:= [];
     for i in [1..d+1] do
        apoly:= Append(apoly,[li[i]]);
     od;
     for i in [d+2..Length(li)] do
        bpoly:= Append(bpoly,[li[i]]);
     od;
     xunten:= Coerce(PA,bpoly);
     ro:=Roots(xunten,SplittingField(xunten))[1][1];
     P:=MakePoints(ChangeBaseRing(E2,Extension(Parent(ro),2)),ro)[1];
     isodu:=IsogenyFromSubgroup(P,E2,d);
     E3:=Codomain(isodu);
     if not E3 = E then
        isodu:= Composition(Isomorphism(E3,E,1),isodu);
     fi;
     isotest:=Composition(isodu,f);
     P:=RandomPoint(E);
     while 2*P = Zero(E) do
         P:=RandomPoint(E);
     od;
     if d*P = isotest(P) then
        return isodu;
     else Assert(d*P = -isotest(P), "dual assumption");
     fi;
     return Negate(isodu);
   end;



#multiplication-by-m-map
   MultMap:= function(E1,m) 
     local muma,mu,mult,P,EE,xmuma,ymuma,xmu,ymu,ymuma2,cha,E,Iso1,Iso2,E2,di,Psi,Psip,Psim,Phi,Psim2,omega,x,y;
     if m = 0 then return NuMap(E1,E1); fi;
     if m = 1 then return Isomorphism(E1,E1,1);fi;
     cha:= Characteristic(BaseRing(E1));
     if cha = 2  then
        E:=NormalFormChar2(E1,1);
        Iso1:=E.ext1;
        E:=Base(E);
        Iso2:=NormalFormChar2(E1,2).ext1;
     else 
        if Length(Coefficients(E1)) = 2 then
           Iso1:=false;
           E:=E1;
        else 
           E:=ShortWNF(E1);
           Iso1:=E.ext1;
           Iso2:=DualIsogeny(Iso1);
           E:=Base(E);
        fi;
     fi;
     P:= GenericPoint(E);
     EE:= Parent(P);
     if cha = 2 then
        Psi:= _ff(E,m);
        Psip:= _ff(E,m+1);
        Psim:= _ff(E,m-1);
        Phi:= Psip*Psim;
        Psim2:=_ff(E,m-2);
        xmuma:= function(P2)
           local psi,phi,BR,PA;
           if Is(Type(PX(P2)),elt-fld^fin) then
              BR:=BaseRing(Parent(P2));
              PA:=PolynomialAlgebra(BR);
              psi:=Coerce(PA,Psi);
              if Evaluate(psi,PX(P2)) = 0 then 
                  return "INFTY POINT"; 
              fi;
              phi:=Coerce(PA,Phi);
              return PX(P2)+Evaluate(phi/psi^2,PX(P2));
           fi;
           if Evaluate(Psi,PX(P2)) = 0 then 
              return "INFTY POINT"; 
           fi;
           return PX(P2)+Evaluate((Phi/Psi^2),PX(P2));
        end;
        ymuma:= function(P2)
           local Phi2,phi2,psi,psim,psip,psim2,BR,PA;
           if Is(Type(PX(P2)),elt-fld^fin) then
              BR:=BaseRing(Parent(P2));
              PA:=PolynomialAlgebra(BR);
              psi:=Coerce(PA,Psi);
              psim:=Coerce(PA,Psim);
              psip:=Coerce(PA,Psip);
              psim2:=Coerce(PA,Psim2);
              phi2:= (PX(P2)^2+PX(P2)+PY(P2))*Evaluate(psim*psi*psip,PX(P2))+Evaluate(psim2*psip^2,PX(P2));
              return PX(P2)+PY(P2)+ phi2/(PX(P2)*Evaluate(psi^3,PX(P2)));
           fi;
           Phi2:= (PX(P2)^2+PX(P2)+PY(P2))*Evaluate(Psim*Psi*Psip,PX(P2))+Evaluate(Psim2*Psip^2,PX(P2));
           return PX(P2)+PY(P2)+ Phi2/(PX(P2)*Evaluate(Psi^3,PX(P2)));
        end;
     else
        Psi:= _ff(E,m);
        Psip:= _ff(E,m+1);
        Psim:= _ff(E,m-1);
        if Is(Type(Psi),elt-fld^fun) then
           x:=CoefficientRing(Parent(Psi)).1;
        else
           x:=PolynomialAlgebra(BaseRing(E)).1;
        fi;
        Phi:= x*Psi^2-Psip*Psim;
        omega:=_ff(E,m+2)*Psim^2-_ff(E,m-2)*Psip^2;
        y:=EllipticFunctionField(E).1;
        xmuma:= function(P2)
           local BR,PA,phi,psi;
           if Is(Type(PX(P2)),elt-fld^fin) then
              BR:=BaseRing(Parent(P2));
              PA:=PolynomialAlgebra(BR);
              psi:=Coerce(PA,Psi);
              if Evaluate(psi,PX(P2)) = 0 then 
                 return "INFTY POINT"; 
              fi;
              phi:=Coerce(PA,Phi);
              return Evaluate(phi/psi^2,PX(P2));
           fi;
           if Evaluate(Psi,PX(P2)) = 0 then 
              return "INFTY POINT"; 
           fi;
           return Evaluate(Phi/Psi^2,PX(P2));
        end;
        ymuma:= function(P2)
          local omega2,psi,BR,PA;
          if Is(Type(PX(P2)),elt-fld^fin) then
             BR:=BaseRing(Parent(P2));
             PA:=PolynomialAlgebra(BR);
             omega2:=Coerce(PA,omega);
             omega2:=Evaluate(omega2,PX(P2))/(4*PY(P2));
             psi:=Coerce(PA,Psi);
             return omega2/(Evaluate(psi,PX(P2))^3);
          fi;
          omega2:= Evaluate(omega,PX(P2))/(4*PY(P2));
          return omega2/(Evaluate(Psi,PX(P2)))^3;
        end;
     fi;
     muma:=function(P1)
        local P2,EE,P;
        if P1 = Zero(E) then 
           return P1;
        fi;
        P:= GenericPoint(E);
        EE:= Parent(P);
        P2:= Coerce(EE,P1);
        xmu:= xmuma(P1);
        ymu:=ymuma(P1);
        if xmu = "INFTY POINT" then 
           return Zero(E);
        fi;
        return MakePoint(E,[xmu,ymu]);
     end;
     mu:= Map(E,E,muma);
     mu.type:=iso;
     mu.operations:= _iso_ops;
     if Characteristic(BaseRing(E))>3 then
        mu.xpol:= Phi/Psi^2;
        mu.ypol:= omega/(4*y*Psi^3);
     else
        x:=PolynomialAlgebra(BaseRing(E)).1; 
        y:=EllipticFunctionField(E).1;
        mu.xpol:=x+Phi/Psi^2;
        mu.ypol:=x+y+(((x^2+x+y)*Psim*Psi*Psip)+Psim2*Psip^2)/(x*Psi^3);
     fi;
     if GCD(cha,Coerce(Integers(),m))<>1 then
        mu.degree:= Degree(mu.xpol)/2;
     else mu.degree:=m^2;
     fi;
     if  Iso1 = false then
        return mu;
     else return 
        Composition(Iso2,Composition(mu,Iso1));
     fi;
   end;

#returns true if Z[pi] maximal order
   EndRingIsMaximal:= function(EC) 
     local q,eq,nop,t,eq,orde,z;
     z:= PolynomialAlgebra(Integers()).1;
     q:= Size(EC.F);
     nop:= Size(EC);
     t:= q+1-nop;
     eq:= EquationOrder(z^2-t*z+q);
     orde:=rec(base:=IsMaximal(eq));
     if orde.base then
        orde.ext1:= eq;
     fi;
     return orde;
   end;







#returns an isomorphism between E1 and E2 over an extension of ground field of degree n
   Isomorphism:= function(E1,E2,n)
     local _Iso,L1,L2,xmap,ymap,P,a,b,A,B,j,delta,u1,u2,Iso,u,i,K,E12,E22,s,a1,a2,x,Iso1,Iso2;
     if not jInvariant(E1) = jInvariant(E2) then
        return Error(E1," and ", E2," are not isomorphic.");
     fi;
     if E1 = E2 then
        Iso:=Map(E1,E2,function(P) return P; end);
        Iso.type:= iso;
        Iso.operations:= _iso_ops;
        P:= GenericPoint(E1);
        Iso.xpol:= PX(P);
        Iso.ypol:= PY(P);
        Iso.degree:= Degree(Iso.xpol)/2;
        AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]);
        AssignNames_(Parent(Iso.ypol),["y"]);
        return Iso;
     fi;
     K:=Base(Extension(BaseRing(E1),n));
     if Characteristic(K)=2 then
        E12:=NormalFormChar2(ChangeBaseRing(E1,K),1);
        Iso1:=E12.ext1;
        E12:=ChangeBaseRing(Base(E12),K);
        E22:=NormalFormChar2(ChangeBaseRing(E2,K),2);
        Iso2:=E22.ext1;
        E22:=ChangeBaseRing(Base(E22),K);
     else 
        E12:=ChangeBaseRing(E1,K);
        E22:=ChangeBaseRing(E2,K);
     fi;
     L1:= Coefficients(E12);
     L2:= Coefficients(E22);
     P:= GenericPoint(E12);
     if Length(L1) =2 then
        a:= L1[1];
        b:= L1[2];
        A:= L2[1];
        B:= L2[2];
        if a = 0 then
           u:=AllRoots(b/B,6);
           if u = Sequence([]) then return Error("no isomorphism defined over BaseRing");
           else
              if u[1] =1 then
                 u := u[2];
              else u:= u[1];
              fi;
           fi;
        else
           u1:=AllRoots(A/a,4);
           if u1 = Sequence([]) then return Error("no isomorphism defined over BaseRing");
           else 
              if b=0 then
                 if u1[1] =1 then
                    u := u1[2];
                 else u:= u1[1];
                 fi;
	      else u2 := AllRoots(B/b,6);
	         if u2 =Sequence([]) then return Error("no isomorphism defined over BaseRing");
	         fi;
                 i:=0;
                 j:=0;
	         while i <Length(u1) do
                     i:=i+1;
	             if not u1[i] =1 then
                        j:=0;
		        while j <Length(u2) do
                            j:=j+1;
		            if u2[j] = u1[i] then
			       u:= u1[i];
			       j:= Length(u2);
			       i:= Length(u1);
			    fi;
		        od;
		     fi;
                 od;
	      fi;
           fi;
        fi;
        xmap:= function(P)
           if PX(P) = INFTY then return INFTY; fi;
           return PX(P)/(u^2);
        end;
        ymap:= function(P)
           if PY(P) = INFTY then return INFTY; fi;
           return PY(P)/(u^3);
        end;
        _Iso:= function(P)
           if P = Zero(E12) then
              return Zero(E22);
           else return Coerce(E22,[xmap(P),ymap(P)]);
           fi;
        end;
        Iso:=Map(E12,E22,_Iso);
        Iso.type:= iso;
        Iso.operations:= _iso_ops;
        Iso.xpol:= xmap(P);
        Iso.ypol:= ymap(P);
        Iso.degree:= Degree(Iso.xpol)/2;
        AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]);
        AssignNames_(Parent(Iso.ypol),["y"]);
        return Iso;
     elif Characteristic(K) =2 then 
        Assert(L1[2]=0 and L1[4]=0," curve in characteristic 2 form ");
        a1:=L1[3];
        a2:=L2[3];
        x:=PolynomialAlgebra(K).1;
        s:=Roots(x^2+x+a1+a2);
        if s = Sequence([]) then 
           return Error("no isomorphism defined over BaseRing");
        else s:=s[1][1];
        fi;
        xmap:= function(P)
           if PX(P) = INFTY then return INFTY; fi;
           return PX(P);
        end;
        ymap:= function(P)
           if PY(P) = INFTY then return INFTY; fi;
           return PY(P)-s*PX(P);
        end;
        _Iso:= function(P)
           if P = Zero(E12) then
              return Zero(E22);
           else return Coerce(E22,[xmap(P),ymap(P)]);
           fi;
        end;
        Iso:=Map(E12,E22,_Iso);
        Iso.type:= iso;
        Iso.operations:= _iso_ops;
        Iso.xpol:= xmap(P);
        Iso.ypol:= ymap(P);
        Iso.degree:= Degree(Iso.xpol)/2;
        AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]);
        AssignNames_(Parent(Iso.ypol),["y"]);
        return Composition(Iso2,Composition(Iso,Iso1));
     else return Error("missing case");
     fi;
   end;


#new base ring is extension of degree n
   ChangeBaseRingN:=function(E,n)
     local K;
     K:=Extension(BaseRing(E),n);
     return ChangeBaseRing(E,K);
   end;

#G subgroup of E(K), computes E/G
   GetCodomain:= function(E,G) 
     local A1,A2,A3,A4,A6,E2,F2,xx,M,a4,a6,a1,a3,a2,b2,b4,b6,nu,i,tq,uq,S,t,w,tQ,E2;
     E2:=Parent(G[1]);
     F2:= [];
     xx:=[];
     M:=[];
     if Length(E.L) = 2 then
       a4:= E.L[1];
       a6:= E.L[2];
       a1:= 0;
       a3:= 0;
       a2:= 0;
     else a1:= E.L[1];
       a3:= E.L[2];
       a2:= E.L[3];
       a4:= E.L[4];
       a6:= E.L[5];
     fi;
     b2:= a1^2+4*a2;
     b4:=a1*a3+2*a4;
     b6:= a3^2+4*a6;
     nu:= Zero(E2);
     for i in [1..Length(G)] do
       if G[i]<>nu then
          if 2*G[i] = nu then
             F2:=Append(F2,[G[i]]);
          elif not G[i].x in xx then
             xx:=Append(xx,[G[i].x]);
             M:=Append(M,[MakePoints(E2,[G[i].x])[1]]);
          fi;
       fi;
     od;
     tq:= function(Q)
        if 2*Q = Zero(Parent(Q)) then
            return (3*Q.x^2+2*a2*Q.x+a4-a1*Q.y);
        else return 6*Q.x^2+b2*Q.x+b4;
        fi;
     end;
     uq:= function(Q)
        return 4*Q.x^3+b2*Q.x^2+2*b4*Q.x+b6;
     end;
     S:= Append(F2,M);
     t:=0;
     w:=0;
     for i in [1..Length(S)] do
        Q:= S[i];
        tQ:= tq(Q);
        t:= t+tQ;
        w:= w+(uq(Q)+PX(Q)*tQ);
     od;
     A1:= a1;
     A2:=a2;
     A3:= a3;
     A4:= a4- 5*t;
     A6:= a6-b2*t-7*w;
     if A1 = 0 and A3 =0 and A2 = 0 then
        E2:= EllipticCurve(BaseRing(E),[A4,A6]);
     else
        E2:= EllipticCurve(BaseRing(E),[A1,A3,A2,A4,A6]);
     fi;
     return E2;
   end;

#the degree of an isogeny
   _IsoDegree:= function(f)
     return f.degree;
   end;

#returns E[l]
   CompleteTorsionSubgroup:=function(E,l)
     local poi,E2,di,ro,nu,i,Sp,PA,F,din,ex,k,li,l,BR;
     di:= DivisionPolynomial(E,l);
     BR:= BaseRing(E);
     Sp:= Extension(SplittingField(di),2);
     di:=Coerce(PolynomialAlgebra(Sp),di);
     ro:= Roots(di);
     E2:= ChangeBaseRing(E,Sp);
     poi:= [];
     for k in [1..Length(ro)] do
        li:= MakePoints(E2,ro[k][1]);
        if (Length(li) <> 2 and not IsCoercible(Integers(),l/2)) then      return Error("something went wrong");
        fi;
        poi:= Append(poi,li);
     od;
     Append_(poi,[Zero(E2)]);
     return poi;
   end;

#return isogeny phi:E ->E/P if P has order 2
   IsoFromSub2:=function(P,E,E2)
     local coeffs,a1,a2,a3,a4,a6,b2,b4,b6,xmap,ymap,_Iso,PP,Zx,EE,polx,poly,Iso,FFchange,f,PA2,xdenom,xnom,xpol,els,f1,f2,ydenom,ynom,y,ypol,PA,P2;
     coeffs:=Coefficients(E);
     if Length(coeffs) = 2 then
        a1:=0;
        a2:=0;
        a3:=0;
        a4:=coeffs[1];
        a6:=coeffs[2];
     else 
        a1:=coeffs[1];
        a2:=coeffs[3];
        a3:=coeffs[2];
        a4:=coeffs[4];
        a6:=coeffs[5];
     fi;
     b2:=a1^2+4*a2;
     b4:=a1*a3+2*a4;
     b6:=a3^2+4*a6;
     PA:=PolynomialAlgebra(BaseRing(E));
     xmap:=function(Q)
        local PR;
        PR:=Coerce(Parent(Q),P);
        return PX(Q)+PX(PR+Q) - PX(PR);
     end;
     ymap:=function(Q)
        local PR;
        PR:=Coerce(Parent(Q),P);
        return PY(Q)+PY(PR+Q) - PY(PR);
     end;
     _Iso:=function(Q)
        if Q = P or Q = Zero(E) then
           return Zero(E2);
        else 
           return MakePoint(E2,[xmap(Q),ymap(Q)]);
        fi;
     end;
     PP:= GenericPoint(Parent(P));
     Zx:= PolynomialAlgebra(Integers());
     EE:= Parent(PP);
     P2:=Coerce(EE,PP);
     polx:= xmap(PP);
     poly:=ymap(PP);
     Iso:= Map(E,E2,_Iso);
     Iso.operations:= _iso_ops;
     Iso.type:= iso;
     FFchange:=true;
     if Parent(polx)= EllipticFunctionField(E) then
        Iso.xpol:= polx;
        Iso.ypol:= poly;
        FFchange:=false;
     else 
        f:=Eltseq(polx)[1];
        PA2:=PolynomialAlgebra(CoefficientRing(Parent(f)));
        xdenom:=Denominator(f);
        xnom:=Coerce(PA,Coerce(PA2,xdenom*f));
        xdenom:=Coerce(PA,xdenom);
        Iso.xnom:=xnom;
        Iso.xdenom:=xdenom;
        xpol:=Coerce(EllipticFunctionField(E),xnom/xdenom);
        els:=Eltseq(poly);
        f1:= els[1];
        f2:= els[2];
        xdenom:=Denominator(f1);
        xnom:=Coerce(PA,Coerce(PA2,xdenom*f1));
        xdenom:=Coerce(PA,xdenom);
        ydenom:=Denominator(f2);
        ynom:=Coerce(PA,Coerce(PA2,ydenom*f2));
        ydenom:=Coerce(PA,ydenom);
        y:=EllipticFunctionField(E).1;
        ypol:=xnom/xdenom+(ynom/ydenom)*y;
        Iso.xpol:=xpol;
        Iso.ypol:=ypol;
     fi;
     Iso.degree:= Degree(Iso.xpol)/2;
     Iso.sub:= [P];
     AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]);
     AssignNames_(Parent(Iso.ypol),["y"]);
     if not FFchange then
        els:= Eltseq(Iso.xpol);
        xdenom:= Denominator(els[1]);
        Iso.xnom:=Coerce(PolynomialAlgebra(BaseRing(Parent(P2))), els[1]*xdenom);
        Iso.xdenom:= xdenom;
     fi;
     return Iso;
   end;


#returns phi:E->E/<P> of degree l
   IsogenyFromSubgroup:= function(P,E,l) 
     local PA,P2,EE,i,Iso,xmap,ymap,E2,_Iso,E,G,L,P,Zx,polx,xdenom,dseq,nseq,xnom,ydenom,ynom,els,BR,poly,f,y,ydenom,ynom,els,f1,f2,FFchange,xpol,ypol,PA2,x,y,mytime, coeffs,a1,a2,a3,a4,a6,b2,b4,b6,R,F2,xcoords,ze,t,w,S,tqtmp,tq,uq,gqy,gqx,xxunten,d,fro;
     fro:=extFrob(E,Degree(BaseRing(Parent(P)),BaseRing(E)));
     if not fro(P) in Subgroup(P) then
        return Error("Subgroup is not galois-invariant");
     fi;
     if P.x = "INFTY POINT" then 
        return Isomorphism(E,E,1);  
     fi;
     if l=2 then
        S:=[P];
     else
        d:=Coerce(Integers(),(l-1)/2);
        S:=[P];
        for i in [2..d] do 
          Append_(S,[i*P]);
        od;
     fi;
     BR:=BaseRing(E);
     PA:=PolynomialAlgebra(BR);
     coeffs:=Coefficients(E);
     if Length(coeffs) = 2 then
        a1:=0;
        a2:=0;
        a3:=0;
        a4:=coeffs[1];
        a6:=coeffs[2];
     else 
        a1:=coeffs[1];
        a2:=coeffs[3];
        a3:=coeffs[2];
        a4:=coeffs[4];
        a6:=coeffs[5];
     fi;
     b2:=a1^2+4*a2;
     b4:=a1*a3+2*a4;
     b6:=a3^2+4*a6;
     gqx:=function(P)
        return 3*PX(P)^2+2*a2*PX(P)+a4-a1*PY(P);
     end;
     gqy:=function(P)
        return -2*PY(P) - a1*PX(P) -a3;
     end;
     uq:=function(P)
        return 4*PX(P)^3+b2*PX(P)^2+2*b4*PX(P)+b6;
     end;
     tq:=function(P)
        if 2*P = Zero(Parent(P)) then
           return gqx(P);
        else 
           return 6*PX(P)^2+b2*PX(P)+b4;
        fi;
     end;
     t:=0;
     w:=0;
     for i in [1..Length(S)] do
        tqtmp:=tq(S[i]);
        t:=t+tqtmp;
        w:=w+uq(S[i])+tqtmp*PX(S[i]);
     od;
     if Length(E.L) =2 then
        E2:= EllipticCurve(BaseRing(E),[a4-5*t,a6-b2*t-7*w]);
     else
        E2:= EllipticCurve(BaseRing(E),[a1,a3,a2,a4-5*t,a6-b2*t-7*w]);
     fi;
     E2.size:=E.size;
     if 2*P = Zero(Parent(P)) then
        return IsoFromSub2(P,E,E2);
     fi;
     P2:=P;
     xmap:= function(P1)
        local sum,i,xx,xxunten;
        xx := PX(P1);
        sum:=xx;
        for i in [1..Length(S)] do
           if sum = INFTY then 
              return sum; 
           fi;
           xxunten:=xx-PX(S[i]);
           if xxunten = 0 then 
              return INFTY;
           fi;
           sum := sum +tq(S[i])/xxunten+uq(S[i])/(xxunten^2);
        od;
        return sum;
     end;
     ymap:= function(P1)
        local sum,i,yy,Q,xx,xxunten;
        xx:=PX(P1);
        yy := PY(P1);
        sum:=yy;
        for i in [1..Length(S)] do
            if sum = INFTY then 
               return sum; 
            fi;
            Q:=S[i];
            xxunten:=xx-PX(Q);
            if xxunten = 0 then 
               return INFTY;
            fi;
            sum := sum -(uq(Q)*((2*yy+a1*xx+a3)/(xxunten)^3) + tq(Q)*((a1*xxunten+yy-PY(Q))/(xxunten)^2)+(a1*uq(Q)-gqx(Q)*gqy(Q))/(xxunten)^2);
        od;
        return sum;
     end;
     _Iso:=function(P1)
        local P2,m;
            if P1 = Zero(Parent(P1)) then  
                return Zero(E2); 
            else 
                if xmap(P1) = INFTY then
                    return Zero(E2);
                fi;
            fi;
            return MakePoint(E2,[xmap(P1),ymap(P1)]);
     end;
     P:= GenericPoint(Parent(P));
     Zx:= PolynomialAlgebra(Integers());
     EE:= Parent(P);
     L:= [];
     for i in [1..Length(S)] do
        S[i]:=Coerce(EE,S[i]);
     od;
     polx:= xmap(P);
     poly:=ymap(P);
     E2.size:=Size(E);
     Iso:= Map(E,E2,_Iso);
     Iso.operations:= _iso_ops;
     Iso.type:= iso;
     FFchange:=true;
     if Parent(polx)= EllipticFunctionField(E) then
         Iso.xpol:= polx;
         Iso.ypol:= poly;
         FFchange:=false;
     else 
         f:=Eltseq(polx)[1];
         PA2:=PolynomialAlgebra(CoefficientRing(Parent(f)));
         xdenom:=Denominator(f);
         xnom:=Coerce(PA,Coerce(PA2,xdenom*f));
         xdenom:=Coerce(PA,xdenom);
         Iso.xnom:=xnom;
         Iso.xdenom:=xdenom;
         xpol:=Coerce(EllipticFunctionField(E),xnom/xdenom);
         els:=Eltseq(poly);
         f1:= els[1];
         f2:= els[2];
         xdenom:=Denominator(f1);
         xnom:=Coerce(PA,Coerce(PA2,xdenom*f1));
         xdenom:=Coerce(PA,xdenom);
         ydenom:=Denominator(f2);
         ynom:=Coerce(PA,Coerce(PA2,ydenom*f2));
         ydenom:=Coerce(PA,ydenom);
         y:=EllipticFunctionField(E).1;
         ypol:=xnom/xdenom+(ynom/ydenom)*y;
         Iso.xpol:=xpol;
         Iso.ypol:=ypol;
     fi;
     Iso.degree:= Degree(Iso.xpol)/2;
     AssignNames_(CoefficientRing(Parent(Iso.xpol)),["x"]);
     AssignNames_(Parent(Iso.ypol),["y"]);
     if not FFchange then
         els:= Eltseq(Iso.xpol);
         xdenom:= Denominator(els[1]);
         Iso.xnom:=Coerce(PolynomialAlgebra(BaseRing(Parent(P2))), els[1]*xdenom);
         Iso.xdenom:= xdenom;
     fi;
     Codomain(Iso).size:=E.size;
     return Iso;
   end;

#frobenius applicable to points of E(GF(q^i))
   extFrob:= function(E,i)
     local E2,BR,F,fro,re;
     BR:= BaseRing(E);
     F:= Extension(BR,i);
     E2:= ChangeBaseRing(E,F);
     fro:= Frobenius(E,1);
     re:= Map(E2,E2,fro.base);
     re.domain:=E2;
     re.codomain:=E2;
     re.type:=iso;
     re.operations:=fro.operations;
     re.fro:=fro.fro;
     re.degree:=fro.degree;
     return re;
   end;


#subfunction for computing isogenies between elliptic curves over fields of  characteristic 2
   B:=function(i,j)
     return (Binomial(i,j) mod 2);
   end;


#returns elliptic curve with j-invariant j and size n
   ECFromJGivenNumber:=function(j,n)
     local E,P,i,O,q,t,l1,l2,E2;
     E:=EllipticCurveFromJInvariant(j);
     O:=Zero(E);
     q:=Size(Parent(j));
     t:=q+1-n;
     P:=RandomPoint(E);
     l1:=n;
     l2:=q+1+t;
     while true do
        if not l1*P = O then
           if not l2*P = O then
              return Error("there is no elliptic curve with j-invariant ",j," and size ",n," over ",Parent(j));
           fi;
           E2:= QuadraticTwist(E);
           E2.size:=n;
           return E2;
        elif not l2*P = O then
           E.size:=n;
           return E;
        else P:=RandomPoint(E);
        fi;
     od;
   end; 


#subfunction for computing isogenies between elliptic curves over fields of  characteristic 2
   Fillp:=function(p,alpha,beta,PAMult,d,x)
     local ok,k,k1,k2,elimvars,sum,sum1,sum2,i,ck,bk,f,coeffs,vars,T,elim,co,varelim,i2;
     ok:=false;
     k:=1;
     k1:=0;
     k2:=1;
     elimvars:=[];
     while k<d-2*k do   
        sum:=0;
        sum2:=0;
        for i in [1..k] do                                
           sum:=sum+(p[i]^2)*(p[d-k+i]^2)*alpha^(2*(i-1));
        od;                                               
        for i in [1..Floor(k/2)] do                       
           sum2:=sum2+p[k-2*i+1]*B(d-k+2*i,i);
        od;                                               
        sum:=AllRoots(alpha,4)[1]*sum+(AllRoots(beta,4)[1]*Sqrt(alpha)^(d+2*k)*sum2);
        ck:=sum/(alpha^(2*k)*AllRoots(alpha,4)[1]);
        bk:=(AllRoots(beta,4)[1]*Sqrt(alpha)^(d+2*k))/(alpha^(2*k)*AllRoots(alpha,4)[1]);
        f:=ck/(bk^2);
        if Is(Type(f),elt-fld^fin) then                   
           f:=MultPol(PAMult,[0],f);
        fi;                                               
        if Is(Type(f),elt-alg^mult) then                  
           coeffs:=Coefficients(f);
           vars:=Variables(f);
           T:=1;
           while T<> 0 do                                 
              T:=0;
              for i in [1..Length(coeffs)] do             
                 if coeffs[i]=0 then                        
                    coeffs[i]:=Coerce(GF(2),0);
                 fi;                                       
                 if Trace(coeffs[i],GF(2)) = 1 then        
                    T:=T+MultPol(PAMult,vars[i],Coerce(GF(2),1));
                 fi;                                       
              od;                                         
              if T<>0 then                                
                 elim :=false;
                 coeffs:=Coefficients(T);
                 vars:=Variables(T);
                 i:=0;
                 while i <Length(vars) do 
                    i:=i+1;                                  
                    if Length(vars[i]) = 1 and vars[i][1]<>0 then 
                       elim :=true;
                       for i2 in [1..Length(vars)] do            
                           if i2 <> i then                       
                              if vars[i][1] in vars[i2] then     
                                  elim:=false;
                              fi;                                
                           fi;                                   
                       od;                                       
                       if elim then                             
                           elimvars:=Append(elimvars,[vars[i][1]]);
                           co:=coeffs[i];
                           varelim:=vars[i][1];
                           Remove_(vars,i);
                           Remove_(coeffs,i);
                           f:=MultPol(PAMult,vars,coeffs)/co;
                           for i2 in [1..Length(p)] do
                               if  (Type(GetEntry(p,i2))=elt-fld^fin or Type(GetEntry(p,i2))=elt-alg^mult) then   
                                  if Is(Type(p[i2]),elt-alg^mult) then      
                                      p[i2]:=Elimination_(p[i2],varelim,f);
                                  fi;    
                               fi;                         
                           od;
                           i:=Length(vars)+1;           
                       fi;
                    fi; 
                 od;       
                 if not elim then
                    return Error("not yet implemented");
                 fi;
                 sum:=0;
                 sum2:=0;
                 for i in [1..k] do
                    sum:=sum+(p[i]^2)*(p[d-k+i]^2)*alpha^(2*(i-1));
                 od;
                 for i in [1..Floor(k/2)] do
                    sum2:=sum2+p[k-2*i+1]*B(d-k+2*i,i);
                 od;
                 sum:=AllRoots(alpha,4)[1]*sum+(AllRoots(beta,4)[1]*Sqrt(alpha)^(d+2*k)*sum2);
                 ck:=sum/(alpha^(2*k)*AllRoots(alpha,4)[1]);
                 bk:=(AllRoots(beta,4)[1]*Sqrt(alpha)^(d+2*k))/(alpha^(2*k)*AllRoots(alpha,4)[1]);
                 f:=ck/(bk^2);
                 coeffs:=Coefficients(f);
                 vars:=Variables(f);
              fi;
           od;           
           sum:=0;
           for i in [1..Length(coeffs)] do
              sum:=sum+MultPol(PAMult,vars[i],Roots(x^2+x+coeffs[i])[1][1]);
           od;  
           sum:=sum*bk;
           p[k+1]:=sum+MultPol(PAMult,[k],bk);
           k:=k+1;
           if not k>d-2*k+1 then
              k1:=k1+1;
              sum1:=0;
              sum2:=0;
              for i in [1..k1] do
                 sum1:=sum1+p[d-2*k1+2*i]*B(d-2*k1-1+2*i,i)*alpha^(2*i);
                 sum2:=sum2+p[d-2*k1+2*(i-1)+1]*B(d-2*k1+2*(i-1),i-1)*alpha^(2*(i-1));
              od;
              sum2:=sum2+p[d+1]*B(d,k1)*alpha^(2*k1);
              p[d-2*k1]:= alpha^(4*k1+1-2*d)*p[k1+1]^4+sum1+alpha*sum2;
              if not k>d-2*k then 
                 k2:=k2+1;
                 sum1:=0;
                 sum2:=0;
                 for i in [1..k2] do 
                    sum1:=sum1+p[d-2*k2+2*i]*B(d+1-2*k2+2*(i-1),(i-1))*alpha^(2*(i-1));
                    sum2:=sum2+p[d-2*k2+2*i+1]*B(d-2*k2+2*i,i)*alpha^(2*i);
                 od;  
                 p[d-2*k2+1]:=p[d-k2+1]^4+alpha*sum1+sum2;
              fi;
           fi; 
        fi;       
     od;
     return [p,elimvars,k,k1,k2];
   end;

#subfunction for computing isogenies between elliptic curves over fields of  characteristic 2
   _Solve:=function(p,m,elimvars,M,BR,PAMult)
     local i,elim,j,le,f,L,c,easy,ellist,k,g,summ1,sum2,Ls,cs,pos,sols,sum1,sum2;
     i:=m;
     elim:=[];
     while i >1  do 
        while i in elimvars do
           i:=i-1; 
           m:=m-1; 
        od;  
        j:=0;
        le:=m-i;
        if i >1 then
           while j <Length(M) and (Length(elim) = le) do
              j:=j+1;
              f:=M[j];
              L:=Variables(f);
              if i in Last(L) then
                 c:=Coefficients(f);
                 easy:=true;
                 ellist:=[];
                 for k in [1..Length(L)-1] do
                    if i in L[k] then
                       easy:=false;
                       Append_(ellist,[k]);
                    fi;
                 od; 
                 if easy then
                    if Length(Last(L)) = 1 then
                       g:= f/Last(c);
                       g:=MultPol(PAMult,Remove(Variables(g),Length(L)),Remove(Coefficients(g),Length(L)));
                    else return Error("case not yet implemented");
                    fi;
                    for k in [1..Length(M)] do
                       M[k]:=Elimination_(M[k],Last(L)[1],g);
                    od;
                    Append_(elim,[[Last(L)[1],g]]);
                    Remove_(M,j);
                    i:=i-1;
                 else
                    Append_(ellist,[k+1]);
                    sum1:=0;
                    sum2:=0;
                    for k in [1..Length(ellist)] do
                       sum1:= sum1+MultPol(PAMult,L[ellist[k]],c[ellist[k]]);
                    od;
                    sum2:=MultPol(PAMult,Variables(f),Coefficients(f))-sum1; 
                    Ls:=Variables(sum1);
                    cs:=Coefficients(sum1);
                    sum1:=0;
                    for k in [1..Length(Ls)] do
                       pos:=Position(Ls[k],i);
                       sum1:=sum1+MultPol(PAMult,Remove(Ls[k],pos),cs[k]);
                    od;
                    g:=MultPol(PAMult,Variables(sum2),Coefficients(sum2),sum1);
                    for k in [1..Length(M)] do
                       M[k]:=Elimination_(M[k],Last(L)[1],g);
                    od;
                    Append_(elim,[[Last(L)[1],g]]);
                    Remove_(M,j);
                    i:=i-1;
                 fi;
              else
                 return Error("not yet implemented");
              fi;
           od;  
        fi;
     od;
     Assert(i=1," all variables except for one eliminated ");
     sols:=[0,1];
     for i in [1..Length(M)] do
        if sols = [] then 
           return Error(" no solution found "); 
        fi;
        j:=0;
        while j <Length(sols) do    
           j:=j+1;
           if "denom" in RecFields(M[i]) then
              if Coerce(BR,Elimination_(M[i],1,sols[j]))/Coerce(BR,Elimination_(M[i].denom,1,sols[j])) <> 0 then
                 Remove_(sols,j);
                 j:=j-1;
              fi;
           else 
              if Coerce(BR,Elimination_(M[i],1,sols[j])) <> 0 then
                 Remove_(sols,j);
                 j:=j-1;
              fi;
           fi; 
        od;
     od;
     Assert(Length(sols) =1 , " only one solution found ");
     sols[1]:=Tuple([1,sols[1]]);
     i:=Length(elim);
     for j in [1..Length(elimvars)] do
        sols[elimvars[j]]:= Tuple([elimvars[j],"nothing"]);
     od;
     while i >=1 do
        k:=elim[i][1];
        f:=elim[i][2];
        if "denom" in RecFields(f) then
           j:=k-1;
           while j in elimvars do
              j:=j-1;
           od;
           while j >1 do
              while j in elimvars do
                 j:=j-1;
              od;
              g:=Elimination_(f,j,sols[j][2]);
              f:=MultPol(PAMult,Variables(g),Coefficients(g),Elimination_(f.denom,j,sols[j][2]));
              j:=j-1;
           od;
           f:=Coerce(BR,Elimination_(f,1,sols[1][2]));
           Assert(f = 1 or f = 0 , " right solution");
           sols[k]:=Tuple([k,f]);
        else 
           j:=k-1;
           while j in elimvars do
              j:=j-1;
           od;
           while j >1 do
              while j in elimvars do
                 j:=j-1;
              od;
              f:=Elimination_(f,j,sols[j][2]);
              j:=j-1;
           od;
           f:=Coerce(BR,Elimination_(f,1,sols[1][2]));
           Assert(f = 1 or f = 0 , " right solution");
           sols[k]:=Tuple([k,f]);
        fi;
        i:=i-1;
     od;
     for i in [1..Length(p)] do
        while not IsConstant(p[i]) do
           k:=Last(Variables(p[i]))[1];
           p[i]:=Elimination_(p[i],k,sols[k][2]);
        od;
     od;
     return p;
   end;


#computes isogeny of degree l, starting from EE, to curve whose j-invariant is the m-th root of the l-th modular polynomial, characteristic 2, Tr(L[3])=0
   IsoFromJChar2new:=function(EE,l,m) 
     local d,L,size,E,Iso1,E4,BR,s,PA,x,PA2,mo,j,ro,PAMult,E2,b,a,alpha,beta,j2,mytime, p,elimvars,k,k1,k2,M,sum1,sum2,i,q,Q,P,spliq,testP,G,K,R,H,xmap,ymap,_Iso,re,PP,fr,lambda;
     if not IsOdd(l) or not IsPrime(l) then 
        return Error("Algorithm implemented only for odd primes"); 
     fi;
     d:=Coerce(Integers(),(l-1)/2);
     L:=Coefficients(EE);
     size:=EE.size;
     E:=NormalFormChar2(EE,1);
     Iso1:=E.ext1;
     E4:=Base(E);
     L:= Coefficients(E4);
     BR:=BaseRing(EE);
     s:=Size(BR);
     E:=E4;
     E.size:=size;
     PA:=PolynomialAlgebra(BR);
     x:=PA.1;
     PA2:=PolynomialAlgebra(PA);
     mo:=Coerce(PA2,ModularPolynomial(l));
     j:=jInvariant(E);
     ro:=Roots(Evaluate(mo,j));
     if Length(ro) = 0 then 
        return Error("There is no ",l,"-isogenous curve");
     fi;
     if m =1 then
        j2:=ro[1][1];
     elif Length(ro) <m then 
        return Error("The modular polynomial does not have k roots");
     else j2:= ro[m][1];
     fi;
     PAMult:=MultPolAlg(BR,Floor(d/3));
     E2:=EllipticCurveFromJInvariant(j2);
     b:=j2^(-1);
     L:= Coefficients(E);
     Assert(L[1]=1 and L[2]=0 and L[3]=0 and L[4] = 0 and L[5] <>0, " curve in normal form");
     a:= L[5];
     alpha:=AllRoots(a,4)[1];
     beta:=AllRoots(b,4)[1];
     p:=[];
     p[d+1]:=1;
     p[d]:=alpha+beta;
     if d >1 then
        if IsOdd(d) then
           p[d-1]:= p[d]^4+alpha*p[d]+alpha^2;
        else 
           p[d-1]:= p[d]^4+alpha*p[d];
        fi;
     fi;
     p[1]:=AllRoots(alpha^(2*d)+alpha^(2*d-1)*p[d],4)[1];
     p:=Fillp(p,alpha,beta,PAMult,d,x);
     elimvars:=p[2];
     k:=p[3];
     k1:=p[4];
     k2:=p[5];
     p:=p[1];
     M:=[];
     m:=k-1;
     if m > 0 then
        for k in [k1+1..Floor((d-1)/2)] do
           sum1:=0;
           sum2:=0;
           for i in [1..k+1] do
              sum1:=sum1+p[d-2*k+2*(i-1)]*B(d-2*k-1+2*(i-1),i-1)*alpha^(2*(i-1));
              sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1));
           od;
           Append_(M,[p[k+1]^4+alpha^(2*d-4*(k)-1)*sum1+alpha^(2*d-4*(k))*sum2]);
        od;
        for k in [k2+1..Floor(d/2)] do
           sum1:=0;
           sum2:=0;
           for i in [1..k] do
              sum1:=sum1+p[d-2*k+2*i]*B(d+1-2*k+2*(i-1),i-1)*alpha^(2*(i-1));
              sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1));
           od;
           sum2:= sum2+p[d+1]*B(d,k)*alpha^(2*k);
           Append_(M,[p[d-k+1]^4+alpha*sum1+sum2]);
        od;
        Assert(IsDense(p)," all coefficients of p computed ");
        p:=_Solve(p,m,elimvars,M,BR,PAMult);
     fi;
     q:=[];
     q[d+1]:=1;
     for i in [1..d+1] do
        q[i]:=Coerce(BR,(AllRoots(alpha,4)[1]/AllRoots(beta,4)[1])*Sqrt(alpha)^(d-2*(i-1))*p[d-(i-1)+1]);
     od;
     Q:=q[1]^2;
     for i in [2..d+1] do
        Q:=Q+q[i]^2*x^(i-1);
     od;
     P:=Coerce(BR,p[1])^2;
     for i in [2..d+1] do
        P:=P+Coerce(BR,p[i])^2*x^(i-1);
     od;
     spliq:=SplittingField(Q);
     testP:=MakePoints(ChangeBaseRing(E,Extension(spliq,2)),HasRoot(Q,spliq).ext1)[1];
     G:=x*P^2;
     K:=(P^2)*Q;
     R:=x*Derivative(P*Q);
     H:=x*R*P+Sqrt(b)*Q^3+Sqrt(a)*P^2*Q;
     xmap:=function(P)
        local xx,q;
        xx:=PX(P);
        q:=Evaluate(Q^2,xx);
        if q = 0 then
           return "INFTY";
        else 
           return Evaluate(G,xx)/q;
        fi;
     end;
     ymap:= function(P)
        local xx,yy;
        xx:=PX(P);
        yy:=PY(P);
        return (Evaluate(H,xx)+yy*Evaluate(K,xx))/Evaluate(Q^3,xx);
     end;
     E2:=ChangeBaseRing(E2,BR);
     E2.size:=E.size;
     _Iso:=function(P)
        local x;
        if P = Zero(E) then
           return Zero(E2);
        else
           x:=xmap(P);
           if x = "INFTY" then
              return Zero(E2);
           else 
              return MakePoint(E2,[x,ymap(P)]);
           fi;
        fi;
     end;
     re:= Map(E,E2,_Iso);
     re.type:= iso;
     re.operations:= _iso_ops;
     PP:=GenericPoint(E);
     re.xpol:=xmap(PP);
     re.ypol:=ymap(PP);
     re.degree:= l;
     i:= Degree(BaseRing(Parent(testP)),BaseRing(E));
     fr:= MakePoint(Parent(testP),[PX(testP)^s,PY(testP)^s]);
     ro:= Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E)));
     lambda:= Coerce(Integers(),ro[1][1]);
     if not (lambda*testP=fr) then 
        lambda:=Coerce(Integers(),ro[2][1]);
        Assert(lambda*testP = fr, " eigenvalue found ");
     fi; 
     re.ei:= lambda;
     re:=Composition(re,Iso1);
     Codomain(re).size:=EE.size;
     return re;
   end;


#computes isogeny of degree l, starting from EE, to curve whose j-invariant is the m-th root of the l-th modular polynomial, characteristic 2
   IsoFromJChar2:=function(EE,l,m) 
     local d,L,E,Iso1,E3,BR,s,PA,x,PA2,mo,j,ro,j2,PAMult,E2,b,alpha,beta,p,a,k,k1,k2,elimvars,M,m,sum1,sum2,i,q,Q,K,P,Iso11,Iso2,Iso,re,Kxy,Alf,defi,Alf2,fr,lambda;
     if not IsOdd(l) or not IsPrime(l) then
        return Error("Algorithm implemented only for odd primes"); 
     fi;
     d:=Coerce(Integers(),(l-1)/2);
     L:=Coefficients(EE);
     if Trace(L[3],GF(2)) = 0 then
        return IsoFromJChar2new(EE,l,m);
     else
        E:=NormalFormChar2(EE,2);
     fi;
     Iso1:=E.ext1;
     E3:=Base(E);
     L:= Coefficients(E3);
     E2:=QuadraticTwist(E3);
     E:=E2;
     BR:=BaseRing(E);
     s:=Size(BR);
     PA:=PolynomialAlgebra(BR);
     x:=PA.1;
     PA2:=PolynomialAlgebra(PA);
     mo:=Coerce(PA2,ModularPolynomial(l));
     j:=jInvariant(E);
     ro:=Roots(Evaluate(mo,j));
     if Length(ro) = 0 then 
        return Error("There is no ",l,"-isogenous curve");
     fi;
     if m =1 then
        j2:=ro[1][1];
     elif Length(ro) <m then
        return Error("The modular polynomial does not have k roots");
     else 
        j2:= ro[m][1];
     fi;
     PAMult:=MultPolAlg(BR,Floor(d/3));
     E2:=EllipticCurveFromJInvariant(j2);
     b:=j2^(-1);
     L:= Coefficients(E);
     Assert(L[1]=1 and L[2]=0 and L[3]=0 and L[4] = 0 and L[5] <>0, " Form of the curve ");
     a:= L[5];
     alpha:=AllRoots(a,4)[1];
     beta:=AllRoots(b,4)[1];
     p:=[];
     p[d+1]:=1;
     p[d]:=alpha+beta; 
     if d >1 then
        if IsOdd(d) then
           p[d-1]:= p[d]^4+alpha*p[d]+alpha^2;
        else 
           p[d-1]:= p[d]^4+alpha*p[d];
        fi;
     fi;
     p[1]:=AllRoots(alpha^(2*d)+alpha^(2*d-1)*p[d],4)[1];
     p:=Fillp(p,alpha,beta,PAMult,d,x);
     elimvars:=p[2];
     k:=p[3];
     k1:=p[4];
     k2:=p[5];
     p:=p[1];
     M:=[];
     m:=k-1;
     if m > 0 then
        for k in [k1+1..Floor((d-1)/2)] do
           sum1:=0;
           sum2:=0;
           for i in [1..k+1] do
              sum1:=sum1+p[d-2*k+2*(i-1)]*B(d-2*k-1+2*(i-1),i-1)*alpha^(2*(i-1));
              sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1));
           od;
           Append_(M,[p[k+1]^4+alpha^(2*d-4*(k)-1)*sum1+alpha^(2*d-4*(k))*sum2]);
        od;
        for k in [k2+1..Floor(d/2)] do
           sum1:=0;
           sum2:=0;
           for i in [1..k] do
              sum1:=sum1+p[d-2*k+2*i]*B(d+1-2*k+2*(i-1),i-1)*alpha^(2*(i-1));
              sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1));
           od;
           sum2:= sum2+p[d+1]*B(d,k)*alpha^(2*k);
           Append_(M,[p[d-k+1]^4+alpha*sum1+sum2]);
        od;
        Assert(IsDense(p)," all coefficients of p computed ");
        p:=_Solve(p,m,elimvars,M,BR,PAMult);
     fi;
     q:=[];
     q[d+1]:=1;
     for i in [1..d+1] do
        q[i]:=Coerce(BR,(AllRoots(alpha,4)[1]/AllRoots(beta,4)[1])*Sqrt(alpha)^(d-2*(i-1))*p[d-(i-1)+1]);
     od;
     Q:=q[1]^2;
     for i in [2..d+1] do
        Q:=Q+q[i]^2*x^(i-1);
     od;
     K:=Extension(SplittingField(Q),2);
     ro:=HasRoot(Q,K);
     P:=MakePoints(ChangeBaseRing(E,K),ro.ext1)[1];
     Iso11:=NormalFormChar2(ChangeBaseRing(EE,K),2);
     Iso2:=Composition(Iso11.ext1,Isomorphism(ChangeBaseRing(E,K),Base(Iso11),1));
     P:=Iso2(P);
     E:= EE; 
     Assert(l*P=Zero(Parent(P))," successful computation of Q ");
     Iso:=IsogenyFromSubgroup(P,E,l);
     E3:=Codomain(Iso);
     re:= Map(ChangeBaseRing(E,BR),ChangeBaseRing(Codomain(Iso),BR),Base(Iso));
     re.type:= iso;
     re.operations:= Iso.operations;
     Kxy:= PolynomialAlgebra(PolynomialAlgebra(BR));
     Alf:= Parent(Iso.xpol);
     defi:= DefiningPolynomial(Alf);
     Alf2:= FunctionField(Coerce(Kxy,defi));
     re.xpol:= Coerce(Alf2,Eltseq(Iso.xpol));
     re.ypol:= Coerce(Alf2,Eltseq(Iso.ypol));
     re.degree:= Iso.degree;
     i:= Degree(BaseRing(Parent(P)),BaseRing(E));
     fr:= MakePoint(Parent(P),[PX(P)^s,PY(P)^s]);
     ro:= Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E)));
     lambda:= Coerce(Integers(),ro[1][1]);
     if not (lambda*P=fr) then 
        lambda:=Coerce(Integers(),ro[2][1]);
        Assert(lambda*P = fr, " eigenvalue found ");
     fi; 
     re.ei:= lambda;
     if not InNormalFormChar2(E3) then
        re:=Composition(NormalFormChar2(E3,1).ext1,re);
        Codomain(re).size:=EE.size;
        return re;
     fi;
     re:=Composition(re,Iso1);
     Codomain(re).size:=EE.size;
     return re;
   end;


#computes isogeny of degree 2, starting from E to curve whose j-invariant is the k-th root of the second modular polynomial
   IsoFromJ2:= function(E,k) #berechnet Isogenie vom Grad 2
     local to,Iso,re,Kxy,Alf,defi,Alf2,l,P,i,R,fr,fro,lambda,eifound,BR,f,E2,j,found,m,mo,fo,ro,fro;
     to:= CompleteTorsionSubgroup(E,2);
     mo:= ModularPolynomial(2);
     j:=Roots(Evaluate(Coerce(PolynomialAlgebra(PolynomialAlgebra(BaseRing(E))),mo),jInvariant(E)))[k][1];
     for i in [1..Length(to)] do
        P:=to[i];
        fro:=extFrob(E,Degree(BaseRing(Parent(P)),BaseRing(E)));
        if fro(P) = P then
           Iso:= IsogenyFromSubgroup(P,E,2);
           if jInvariant(Codomain(Iso)) = j then 
              Iso.ei:=1;
              return Iso; 
           fi;
        fi;
     od;
     Iso:= FrobeniusToOtherCurve(E,2);
     Iso.ei:=0;
     Codomain(Iso).size:=E.size;
     return Iso;
   end;


#subfunction for endomorphism ring computation, computes conductor of End(E)
   FindEndRingIndexFromEqOrder:= function(E)
     local jE,f,O,i,l,fact,a,m,ind,d,divpol,x,BR,divap,divam,divas,h,las,p,j,r,t,atmp,modpoly,PA,g,nr,ro,n,co,n1,n2,mi,n3,prev,p,n2set,n,to,fro,k,atmp,P,done;
     f:= CharacteristicFrobeniusPolynomial(E);
     t:= -Coefficient(f,1);
     O:= EquationOrder(f);
     m:= Conductor(O);
     if m = 1 then
        return 1;
     fi;
     d:= Discriminant(MaximalOrder(O));
     if d mod 4 = 0 then
        a:= t/2;
     elif d mod 4 = 1 then
        a:= (t+m)/2;
     else return Error("something went wrong");
     fi;
     a:= Coerce(Integers(),a);
     ind:=1;
     fact:= Factorization(m);
     BR:= BaseRing(E);
     x:= PolynomialAlgebra(BR).1;
     p:= Size(BR);
     PA:= PolynomialRing(PolynomialRing(BR));
     jE:= jInvariant(E);
     prev:=[];
     n2set:=false;
     for i in [1..Length(fact)] do
        l:= fact[i][1];
        r:= fact[i][2];
        n:=0;
        j:=jE;
        modpoly:= Coerce(PA,ModularPolynomial(l));
        las:=false;
        while n <=r and not las do
           prev:=Append(prev,[j]);
           n:=n+1;
           g:= Evaluate(modpoly,j); 
           ro:= Roots(g);
           nr:= 0;
           for co in [1..Length(ro)] do
              nr:= nr+ ro[co][2];
           od;
           if nr <> l+1 then
              las:= true;
              n:=n-1;
           else 
              j:= ro[1][1];
              p:=1;
              while (j in prev or j =0) and not las do
                 p:= p+1;
                 if p>Length(ro) then
                    las:=true;
                    n:=n-1;
                 else 
                    j:= ro[p][1];
                 fi;
              od;
           fi;
        od;
        n1:=n;
        n:=0;
        j:=jE;
        las:=false;
        prev:=[];
        while n <=r and not las do
           prev:=Append(prev,[j]) ;
           n:=n+1;
           g:= Evaluate(modpoly,j); 
           ro:= Roots(g);
           nr:= 0;
           for co in [1..Length(ro)] do
              nr:= nr+ ro[co][2];
           od;
           if nr <> l+1 then
              las:= true;
              n:=n-1;
           else
              j:= ro[2][1];
              p:= 1;
              while (j in prev or j = 0) and n2set=false  do
                 p:=p+1;
                 if p>Length(ro) then
                    n2:=n1;
                    n2set:=true;
                 else
                    j:= ro[p][1];
                 fi;
              od;
           fi;
        od;
        if not n2set then
           n2:= n;
        fi;
        if not (n1>r and n2>r) then 
           mi:= Minimum([n1,n2]);
           if not mi = 0 then las:= true;
             ind:=ind*l^mi;
           fi;
        else ind:= ind*(l^r);    
        fi;
     od;
     return Coerce(Integers(),m/ind);
   end;


#computes End(E)
   EndRing:= function(E)
     local c,Omax,f,g,O,OE,alpha,beta,fact,D,x;
     if IsSupersingular(E) then 
        return Error("Cannot compute Endomorphism Ring for supersingular curve");
     fi;
     f:= CharacteristicFrobeniusPolynomial(E);
     O:= EquationOrder(f);
     c:=Conductor(O);
     if not c = 1 then
        fact:=Factorization(c);
        if Last(fact)[1] > 59 then
           return FAILURE;
        fi;
     fi;
     c:= FindEndRingIndexFromEqOrder(E);
     Omax:= MaximalOrder(O);
     OE:=EquationOrder(MinimalPolynomial(c*Omax.2));
     Assert(Conductor(OE) = c, "Order has expected conductor");
     return OE;
   end;



#computes isogeny of degree l, starting from E, to curve whose j-invariant is k-th root of l-th modular polynomial
   IsogenyFromJInvariant:= function(E,l,k)
     local sum,jTilDer,p1,coeffs,k,PA2,m,C,A,K,CFun,ro,AFun,unf,gen,bo,f, xval,sub,BR,jDer,a,b,ModPoly,PA,j2,j,E2,a,b,a2,d,b2,Phix,Phiy,Phixx,Phiyy,Phixy,J,E4,E6,E4Til,E6Til,c,cTil,i,j,sum1,sum2,divpolyfact,insum,Kxy,Alf,Alf2,defi,Iso,re,R,fro,fr,P,lambda,eifound,E2,i,ro,r, mytime,esize,xx,ff,kk,sumA,jj,f2;
     esize:=Size(E);
     if l = 2 then 
        return IsoFromJ2(E,k);
     fi;
     if Characteristic(BaseRing(E)) = 2 then
        return IsoFromJChar2(E,l,k);
     fi;
     if IsSupersingular(E) then
        return Error("cannot compute isogeny for supersingular curve");
     fi;
     Assert(Length(Coefficients(E)) =2, "curve in short WeierstrassNormalForm");
     a:= E.L[1];
     b:= E.L[2];
     c:=[];
     cTil:=[];
     BR:= BaseRing(E);
     j:=jInvariant(E);
     if j = Coerce(BR,0) or j = Coerce(BR,1728) then
        return Error("cannot compute isogeny for j-Invariant ",j);
     fi;
     PA:=PolynomialAlgebra(BR);
     PA2:= PolynomialAlgebra(PA);
     ModPoly:= Coerce(PA2,ModularPolynomial(l));
     ro:=Roots(Evaluate(ModPoly,j));#Print(ro, "\n");
     if Length(ro) = 0 then 
        return Error("There is no ",l,"-isogenous curve");
     fi;
     E:= ChangeBaseRing(E,BR);
     E.size:=esize;
     PA:= PolynomialAlgebra(BR);
     if k =1 then
        j2:=ro[1][1];
     elif Length(ro) <k then 
        return Error("The modular polynomial does not have k roots");
     else
        j2:= ro[k][1];
     fi;
     if j2 = 0 or j2 = 1728 then
        return FAILURE;
     fi;
     Phix:= Derivative(ModPoly);
     if Evaluate(Evaluate(Phix,j),j2) = 0 then
        return FAILURE;
     fi;
     Phiy:=(Derivative(_SwapVars(ModPoly)));
     Phiyy:= _SwapVars(Derivative(Phiy));
     Phiy:= _SwapVars(Phiy);
     Phixx:= Derivative(Phix);
     Phixy:= _SwapVars(Derivative(_SwapVars(Phix)));
     E4:=-48*a;
     E6:= 864*b;
     jDer:= -j*E6/E4;
     jTilDer:= -((jDer*Evaluate(Evaluate(Phix,j),j2))/(l*Evaluate(Evaluate(Phiy,j),j2)));
     jTilDer:= Coerce(BR,jTilDer);
     a2:= -(jTilDer^2)/(48*j2*(j2-1728));
     b2:= -(jTilDer^3)/(864*j2^2*(j2-1728));
     E4Til:= -48*a2;
     E6Til:= 864*b2;
     J:=(jDer^2*Evaluate(Evaluate(Phixx,j),j2));
     J:= J+2*l*jDer*jTilDer*Evaluate(Evaluate(Phixy,j),j2);
     J:= J+l^2*jTilDer^2*Evaluate(Evaluate(Phiyy,j),j2);
     J:= -J/(jDer*Evaluate(Evaluate(Phix,j),j2));
     p1:= l/2*J+l/4*(E4^2/E6-l*E4Til^2/E6Til)+l/3*(E6/E4-l*E6Til/E4Til);
     E2:=EllipticCurve(BR,[l^4*a2,l^6*b2]);
     c[1]:=-a/5;
     cTil[1]:=-(l^4*a2)/5;
     c[2]:= -b/7;
     cTil[2]:= -(l^6*b2)/7;
     d:= Coerce(Integers(),(l-1)/2);
     for i in [3..d+1] do
        sum1 :=0;
        sum2:=0;
        for j in [1..i-2] do
           sum1:= sum1 + c[j]*c[i-1-j];
	   sum2:= sum2 + cTil[j]*cTil[i-1-j];
        od;
        c[i] :=  sum1*3/((i-2)*(2*i+3));
        cTil[i] :=  sum2*3/((i-2)*(2*i+3));
     od;
     f:= Coerce(PA,Append([0],c));
     C:= function(k,j)
        local f2,i,elseq;
        f2:=f^(k-j);
        elseq:= Eltseq(f2);
        if elseq=Sequence([]) then 
           return 0;
        elif Length(elseq)<=j then 
           return 0;
        else
           return elseq[j+1];
        fi;
     end;
     A:=[];
     xx:= PA.1;
     ff:=-1/2*p1*xx;
     for kk in [1..d] do
        ff:= ff-(cTil[kk]-l*c[kk])/((2*kk+1)*(2*kk+2))*xx^(kk+1);
     od;
     sumA:=0;
     for jj in [0..d] do
        sumA:= sumA+ff^jj/(Factorial(jj));
     od;
     A:= Eltseq(sumA);
     coeffs:=[];
     coeffs[d+1]:=1;
     for i in [1..d] do
        m:= d-i+1;
        sum:=0;
        for k in [1..i] do
           insum:=0;
           for j in [0..k] do
	      insum:= insum+Binomial(d-i+k,k-j)*C(k,j);
	   od;
           sum:=sum+insum*coeffs[d-i+k+1];
        od;
        coeffs[m]:=A[i+1]-sum;
     od;
     divpolyfact := Coerce(PA,coeffs);
     if Characteristic(BaseRing(E))>6*d+1 then
        return IsogenyFromKernel(E,E2,divpolyfact);
     fi;
     K := SplittingField( divpolyfact );
     PA := PolynomialAlgebra(K);
     xval:= HasRoot(Coerce(PA,divpolyfact)).ext1;
     E2:= ChangeBaseRing(E,Extension(K,2));
     E:= ChangeBaseRing(E,BR);
     E.size:=esize;
     P:=MakePoints(E2,xval)[1];
     Assert( l * P = Zero(E2), "l * P = O");
     Iso:=IsogenyFromSubgroup(P,E,l);
     re:= Map(E,ChangeBaseRing(Codomain(Iso),BR),Base(Iso));
     Codomain(re).size:=esize;
     re.type:= iso;
     re.operations:= Iso.operations;
     Kxy:= PolynomialAlgebra(PolynomialAlgebra(BR));
     Alf:= Parent(Iso.xpol);
     defi:= DefiningPolynomial(Alf);
     Alf2:= FunctionField(Coerce(Kxy,defi));
     re.xpol:= Coerce(Alf2,Eltseq(Iso.xpol));
     re.ypol:= Coerce(Alf2,Eltseq(Iso.ypol));
     re.degree:= Iso.degree;
     i:= Degree(BaseRing(Parent(P)),BaseRing(E));
     if i =1 then 
        re.ei:=1;
     else
        fro:= extFrob(E,i);
        fr:= fro(P);
        ro:= Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E)));
        lambda:= Coerce(Integers(),ro[1][1]);
        if not (lambda*P=fr) then 
           lambda:=Coerce(Integers(),ro[2][1]);
           Assert(lambda*P = fr, " eigenvalue found ");
        fi; 
        re.ei:= lambda;
     fi;
     Codomain(re).size:=E.size;
     return re;
   end;


#computes eigenvalue of frobenius on subgroup in kernel of isogeny of degree 2, starting from E to curve whose j-invariant is the k-th root of the second modular polynomial
   SubgroupFromJ2:=function(E,k)
     local to,j,j2,ModPoly,PA,PA2,ro,Iso,BR,i;
     to:= TorsionSubgroup(E,2);
     j:=jInvariant(E);
     BR:=BaseRing(E);
     ModPoly:= ModularPolynomial(2);
     PA:=PolynomialAlgebra(BR);
     PA2:= PolynomialAlgebra(PA);
     ModPoly:= Coerce(PA2,ModPoly);
     ro:=Roots(Evaluate(ModPoly,j));
     j2:= ro[k][1];
     i:=0;
     while i<(Length(to) - 1) do
        i:=i+1;
        Iso:=IsogenyFromSubgroup(to[i],E,2);
        if jInvariant(Codomain(Iso)) = j2 then
           return [2,1];
        fi;
     od;
     if GCD(2,Characteristic(BR))<>1 then
        Iso:=FrobeniusToOtherCurve(E,2);
        if jInvariant(Codomain(Iso)) = j2 then
           return [2,0];
        fi;
     fi;
     return Error("something went wrong");
   end;

#computes eigenvalue of frobenius on subgroup in kernel of isogeny of degree l, starting from E to curve whose j-invariant is the m-th root of the l-th modular polynomial, characteristic of base ring is 2
   SubFromJChar2:=function(EE,l,m) 
     local L,a,j,mo,PA,PA2,BR,j2,ro,K,k1,P,k2,d,E2,b,q,alpha,beta,x,f,ok,k,i,g,pi,b1,xlist,PAtmp,xtmp,q,Q,Iso,E3,p,sum,sum2,co,ck,bk,PAMult,coeffs,T,vars,fel,i2,elim,sum1,M,PXM,V,W,n,pi,L2,c,easy,ellist,sols,le,Ls,cs,pos,re,Kxy,Alf,defi,Alf2,fr,sub,fro,lambda,r,s,E,quad,Iso1,Iso2,E4,Iso11,EC,elimlist,elimvars,mytime,ro,lambda,lambdax,li,de,nu,roro,poi,xval,gen;
     if not IsOdd(l) or not IsPrime(l) then
        return Error("Algorithm implemented only for odd primes"); 
     fi;
     d:=Coerce(Integers(),(l-1)/2);
     L:=Coefficients(EE);
     if Trace(L[3],GF(2)) = 0 then
        E:=NormalFormChar2(EE,1);
     else
        E:=NormalFormChar2(EE,2);
     fi;
     Iso1:=E.ext1;
     E4:=Base(E);
     L:= Coefficients(E4);
     quad:=L[3]<>0;
     if quad then
        quad:=-1;
        E2:=QuadraticTwist(E4);
        E:=E2;
     else 
        E:=E4;
        quad:=1;
     fi;
     BR:=BaseRing(E);
     s:=Size(BR);
     PA:=PolynomialAlgebra(BR);
     x:=PA.1;
     PA2:=PolynomialAlgebra(PA);
     mo:=Coerce(PA2,ModularPolynomial(l));
     j:=jInvariant(E);
     ro:=Roots(Evaluate(mo,j));
     if Length(ro) = 0 then 
        return Error("There is no ",l,"-isogenous curve");
     fi;
     if m =1 then
        j2:=ro[1][1];
     elif Length(ro) <m then
        return Error("The modular polynomial does not have k roots");
     else j2:= ro[m][1];
     fi;
     EC:=ECFromJGivenNumber(j2,Size(EE));
     PAMult:=MultPolAlg(BR,Floor(d/3));
     E2:=EllipticCurveFromJInvariant(j2);
     b:=j2^(-1);
     L:= Coefficients(E);
     if not L[1]=1 and L[2]=0 and L[3]=0 and L[4] = 0 and L[5] <>0 then
        return Error("Algorithm implemented for curve of the form Y^2+XY=X^3+a");
     fi;
     a:= L[5];
     alpha:=AllRoots(a,4)[1];
     beta:=AllRoots(b,4)[1];
     p:=[];
     p[d+1]:=1;
     p[d]:=alpha+beta;
     if d>1 then
        if IsOdd(d) then
           p[d-1]:= p[d]^4+alpha*p[d]+alpha^2;
        else 
           p[d-1]:= p[d]^4+alpha*p[d];
        fi;
     fi;
     p[1]:=AllRoots(alpha^(2*d)+alpha^(2*d-1)*p[d],4)[1];
     p:=Fillp(p,alpha,beta,PAMult,d,x);
     elimvars:=p[2];
     k:=p[3];
     k1:=p[4];
     k2:=p[5];
     p:=p[1];
     M:=[];
     m:=k-1;
     if m > 0 then
        for k in [k1+1..Floor((d-1)/2)] do
           sum1:=0;
           sum2:=0;
           for i in [1..k+1] do
              sum1:=sum1+p[d-2*k+2*(i-1)]*B(d-2*k-1+2*(i-1),i-1)*alpha^(2*(i-1));
              sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1));
           od;
           Append_(M,[p[k+1]^4+alpha^(2*d-4*(k)-1)*sum1+alpha^(2*d-4*(k))*sum2]);
        od;
        for k in [k2+1..Floor(d/2)] do
           sum1:=0;
           sum2:=0;
           for i in [1..k] do
              sum1:=sum1+p[d-2*k+2*i]*B(d+1-2*k+2*(i-1),i-1)*alpha^(2*(i-1));
              sum2:=sum2+p[d-2*k+2*(i-1)+1]*B(d-2*k+2*(i-1),i-1)*alpha^(2*(i-1));
           od;
           sum2:= sum2+p[d+1]*B(d,k)*alpha^(2*k);
           Append_(M,[p[d-k+1]^4+alpha*sum1+sum2]);
        od;
        Assert(IsDense(p)," all coefficients of p computed ");
        p:=_Solve(p,m,elimvars,M,BR,PAMult);
     fi;
     q:=[];
     q[d+1]:=1;
     for i in [1..d+1] do
        q[i]:=Coerce(BR,(AllRoots(alpha,4)[1]/AllRoots(beta,4)[1])*Sqrt(alpha)^(d-2*(i-1))*p[d-(i-1)+1]);
     od;
     Q:=q[1]^2;
     for i in [2..d+1] do
        Q:=Q+q[i]^2*x^(i-1);
     od;
     ro:=Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E)));
     if Length(ro) = 2 and ro[1][1] = - ro[2][1] then
        K := SplittingField( Q );
        PA := PolynomialAlgebra(K);
        xval:= HasRoot(Coerce(PA,Q)).ext1;
        E2:= ChangeBaseRing(E,Extension(K,2));
        P:=MakePoints(E2,xval)[1];
        fro:=extFrob(E,Degree(BaseRing(E2),BR));
        lambda:=Coerce(Integers(),ro[1][1]);
        if lambda*P = fro(P) then
           return [l,quad*lambda];
        else 
           return [l,quad*(Coerce(Integers(),ro[2][1]))];
        fi;
     fi;
     lambda:=ro[1][1];
     if lambda = 1 then
        roro:= HasRoot(Q);
        if roro then
           poi:=MakePoints(E,roro.ext1);
           if poi <> [] then
             return [l,quad];
           else 
             return [l,quad*Coerce(Integers(),ro[2][1])];
           fi;
        else return [l,quad*Coerce(Integers(),ro[2][1])];
        fi;
     fi;
     lambdax:=MultMap(E,lambda).xpol;
     li:=List(lambdax);
     de:=Denominator(li[1]);
     nu:=de*li[1];
     de:=Coerce(Parent(Q),de);
     nu:=Coerce(Parent(Q),nu);
     if not nu mod Q = de*(Modexp(x,Size(BR),Q)) mod Q then
        lambda:=ro[2][1];
     fi;
     return [l,quad*Coerce(Integers(),lambda)];
   end;

#computes eigenvalue of frobenius on subgroup in kernel of isogeny of degree l, starting from E to curve whose j-invariant is the k-th root of the second modular polynomial
   SubgroupFromJInvariant:= function(E,l,k)
     local sum,jTilDer,p1,coeffs,k,PA2,m,C,A,K,CFun,ro,AFun,unf,gen,bo,f,ENeu,xx,ff,kk,sumA,jj, xval,sub,BR,jDer,a,b,ModPoly,PA,j2,j,E2,a,b,a2,d,b2,Phix,Phiy,Phixx,Phiyy,Phixy,J,E4,E6,E4Til,E6Til,c,cTil,i,j,sum1,sum2,divpolyfact,insum,Kxy,Alf,Alf2,defi,Iso,re,R,fro,fr,P,lambda,eifound,E2,i,gen,ro,lambda,lambdax,li,de,nu,x,roro,poi;
     Assert(IsPrime(l)," computing prime order subgroup");
     if IsSupersingular(E) then
        return Error("cannot compute isogeny for supersingular curve");
     fi;
     if l = 2 then 
        return SubgroupFromJ2(E,k); 
     fi;
     if Characteristic(BaseRing(E))=2 then
        return SubFromJChar2(E,l,k);
     fi;
     Assert(Length(E.L)=2,"curve in short WNF");
     a:= E.L[1];
     b:= E.L[2];
     c:=[];
     cTil:=[];
     BR:= BaseRing(E);
     j:=jInvariant(E);
     if j = Coerce(BR,0) or j = Coerce(BR,1728) then
        return Error("cannot compute isogeny for j-Invariant ",j);
     fi;
     ModPoly:= ModularPolynomial(l);
     PA:=PolynomialAlgebra(BR);
     PA2:= PolynomialAlgebra(PA);
     ModPoly:= Coerce(PA2,ModPoly);
     ro:=Roots(Evaluate(ModPoly,j));#Print(ro, "\n");
     if Length(ro) = 0 then
        return Error("There is no ",l,"-isogenous curve");
     fi;
     PA:= PolynomialAlgebra(BR);
     if k =1 then
        j2:=ro[1][1];
     elif Length(ro) =1 then
        return Error("The modular polynomial has only one root, try k=1");
     else j2:= ro[2][1];
     fi;
     Phix:= Derivative(ModPoly);
     if Evaluate(Evaluate(Phix,j),j2) = 0 then
        return FAILURE;
     fi;
     Phiy:=(Derivative(_SwapVars(ModPoly)));
     Phiyy:= _SwapVars(Derivative(Phiy));
     Phiy:= _SwapVars(Phiy);
     Phixx:= Derivative(Phix);
     Phixy:= _SwapVars(Derivative(_SwapVars(Phix)));
     E4:=-48*a;
     E6:= 864*b;
     jDer:= -j*E6/E4;
     jTilDer:= -((jDer*Evaluate(Evaluate(Phix,j),j2))/(l*Evaluate(Evaluate(Phiy,j),j2)));
     jTilDer:= Coerce(BR,jTilDer);
     a2:= -(jTilDer^2)/(48*j2*(j2-1728));
     b2:= -(jTilDer^3)/(864*j2^2*(j2-1728));
     E4Til:= -48*a2;
     E6Til:= 864*b2;
     J:=(jDer^2*Evaluate(Evaluate(Phixx,j),j2));
     J:= J+2*l*jDer*jTilDer*Evaluate(Evaluate(Phixy,j),j2);
     J:= J+l^2*jTilDer^2*Evaluate(Evaluate(Phiyy,j),j2);
     J:= -J/(jDer*Evaluate(Evaluate(Phix,j),j2));
     p1:= l/2*J+l/4*(E4^2/E6-l*E4Til^2/E6Til)+l/3*(E6/E4-l*E6Til/E4Til);
     ENeu:=EllipticCurve(BR,[l^4*a2,l^6*b2]);
     c[1]:=-a/5;
     cTil[1]:=-(l^4*a2)/5;
     c[2]:= -b/7;
     cTil[2]:= -(l^6*b2)/7;
     d:= Coerce(Integers(),(l-1)/2);
     for i in [3..d+1] do
        sum1 :=0;
        sum2:=0;
        for j in [1..i-2] do
           sum1:= sum1 + c[j]*c[i-1-j];
	   sum2:= sum2 + cTil[j]*cTil[i-1-j];
        od;
        c[i] :=  sum1*3/((i-2)*(2*i+3));
        cTil[i] :=  sum2*3/((i-2)*(2*i+3));
     od;
     C:= function(k,j)
        local f,i,elseq;
        f:= Coerce(PA,Append([0],Sublist(c,[1..j])));
        f:=f^(k-j);
        elseq:= Eltseq(f);
        if elseq=Sequence([]) then 
           return 0;
        elif Length(elseq)<=j then
           return 0;
        else
           return elseq[j+1];
        fi;
     end;
     A:=[];
     xx:= PA.1;
     ff:=-1/2*p1*xx;
     for kk in [1..d] do
        ff:= ff-(cTil[kk]-l*c[kk])/((2*kk+1)*(2*kk+2))*xx^(kk+1);
     od;
     sumA:=0;
     for jj in [0..d] do
        sumA:= sumA+ff^jj/(Factorial(jj));
     od;
     A:= Eltseq(sumA);
     coeffs:=[];
     coeffs[d+1]:=1;
     for i in [1..d] do
        m:= d-i+1;
        sum:=0;
        for k in [1..i] do
           insum:=0;
           for j in [0..k] do
	      insum:= insum+Binomial(d-i+k,k-j)*C(k,j);
	   od;
           sum:=sum+insum*coeffs[d-i+k+1];
        od;
        coeffs[m]:=A[i+1]-sum;
     od;
     divpolyfact := Coerce(PA,coeffs);
     ro:=Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E)));
     if Length(ro) = 2 and ro[1][1] = - ro[2][1] then
        K := SplittingField( divpolyfact );
        PA := PolynomialAlgebra(K);
        xval:= HasRoot(Coerce(PA,divpolyfact)).ext1;
        E2:= ChangeBaseRing(E,Extension(K,2));
        P:=MakePoints(E2,xval)[1];
        fro:=extFrob(E,Degree(BaseRing(E2),BR));
        lambda:=Coerce(Integers(),ro[1][1]);
        if lambda*P = fro(P) then
           return [l,lambda];
        else 
           return [l,Coerce(Integers(),ro[2][1])];
        fi;
     fi;
     lambda:=Coerce(Integers(),ro[1][1]);
     if lambda = 1 then
        roro:= HasRoot(divpolyfact);
        if roro then
           poi:=MakePoints(E,roro.ext1);
           if poi <> [] then 
              return [l,1];
           else 
              return [l,Coerce(Integers(),ro[2][1])];
           fi;
        else 
           return [l,Coerce(Integers(),ro[2][1])];
        fi;
     fi;
     lambdax:=MultMap(E,lambda).xpol;
     li:=List(lambdax);
     de:=Denominator(li[1]);
     nu:=de*li[1];
     de:=Coerce(Parent(divpolyfact),de);
     nu:=Coerce(Parent(divpolyfact),nu);
     x:=PolynomialAlgebra(BR).1;
     if not nu mod divpolyfact = de*(Modexp(x,Size(BR),divpolyfact)) mod divpolyfact then
        lambda:=ro[2][1];
     fi;
     return [l,Coerce(Integers(),lambda)];
   end;




        

#computes isogeny of degree 2 from an ideal
   IsoIdeals2:= function(fact,E1,E2,isoli)
      local E,n,mo,ma,I,i,lev,k,jtmp,j,j2,ro,path,Iso,c,O;
      if Length(isoli)=0 then
         E:=E1;
      else
         E:= Codomain(isoli[Length(isoli)]);
      fi;
      O:=EndRing(E);
      c:= Conductor(O);
      n:= fact[2];
      I:=fact[1];
      mo:= Coerce(PolynomialAlgebra(PolynomialAlgebra(BaseRing(E))),ModularPolynomial(2));
      j:=jInvariant(E);
      j2:= jInvariant(E2);
      Extend(j);
      j.prev:=false;
      ro:=Roots(Evaluate(mo,j));
      j.next:=[];
      for i in [1..Length(ro)] do
        if Conductor(EndRing(EllipticCurveFromJInvariant(ro[i][1])))=c then
           Append_(j.next,[ro[i][1]]);
           k:=Length(j.next);
           Extend(j.next[k]);
           j.next[k].level:=1;
           j.next[k].prev:=j;
           j.next[k].ro:=i;
           j.next[k].checked:=0;
           j.checked:=0;
        fi;
     od;
     j.level:=0;
     while not j.base=j2 do
        i:= Length(j.next);
        if j.checked=i then 
           j:=j.prev;
        else
           lev:=j.level;
           j.checked:=j.checked+1;
           j:= j.next[j.checked];
           j.level:=lev+1;
           if not j.level =n then
              ro:=Roots(Evaluate(mo,j));
              j.next:=[];
              for i in [1..Length(ro)] do
                 if Conductor(EndRing(EllipticCurveFromJInvariant(ro[i][1])))=c then
                    Append_(j.next,[ro[i][1]]);
                    k:= Length(j.next);
                    Extend(j.next[k]);
                    jtmp:=j.next[k];
                    jtmp.level:=lev+2;
                    jtmp.prev:=j;
                    jtmp.ro:=i;
                    jtmp.checked:=0;
                 fi;
              od;
           elif j.base = j2 then
              path:=[];
              path[1]:=Tuple([j.base,j.ro]);
              while not j.prev = false do
                 j:=j.prev;
                 if not j.prev = false then
                    Append_(path,[Tuple([j.base,j.ro])]);
                 fi;
              od;
              j.base:=j2;
           else 
              j:=j.prev; 
           fi;
        fi;
     od;
     if path=[] then
        path[1]:=Tuple([j.base,j.ro]);
        while not j.prev = false do
           j:=j.prev;
           if not j.prev=false then
              Append_(path,[Tuple([j.base,j.ro])]);
           fi;
        od;
     fi;
     path:= Reversed(path);
     for i in [1..Length(path)] do
        Append_(isoli,[IsoFromJ2(E,path[i][2])]);
        E:= Codomain(isoli[Length(isoli)]);
     od;
     if not E = E2 then 
        Append_(isoli,[Isomorphism(E,E2,1)]);
     fi;
     Iso:=isoli[1];
     for i in [2..Length(isoli)] do
        Iso:=Composition(isoli[i],Iso);
     od;
     return Iso;
   end;

#returns eigenvalue of frobenius on kernel subgroup of f, if it has been assigned
   FrobeniusEigenvalue:=function(f)
     if "ei" in RecFields(f) then
        return f.ei;
     else
        return Error("eigenvalue of frobenius has not been assigned to ", f);
     fi;
   end;


#computes isogeny between E1 and E2 from a factored product of ideals
   IsoFromIdeals:= function(fact,E1,E2)
     local bo,nr2,O,PA,ch,pi,I,r,no, elk,Iso,j,f,BR,l,i,P,isoli,n,lambda,Iso,fro,to,E,li,iso1,iso2,j,ro,k,m,tmp,c,PA2,O,c,pi1,ff,Omax1,OE,Omax,pi,rlist,to2,ll;
     nr2:=false;
     E:= E1;
     isoli:=[];
     j:= jInvariant(E2);
     ff:=CharacteristicFrobeniusPolynomial(E1);
     O:=EndRing(E);
     c:=Conductor(O);
     if c > 1 then
        pi1:=EquationOrder(ff).2;
        Omax1:=MaximalOrder(Parent(pi1));
        OE:=EquationOrder(MinimalPolynomial(c*Omax1.2));
        Omax:=MaximalOrder(OE);
        ff:=Coerce(PolynomialAlgebra(Omax),ff);
        pi:=Roots(ff);
        rlist:=List(pi[1][1]);
        if rlist[2]>0 then
           pi:=Coerce(OE,pi[1][1]);
        else 
           pi:=Coerce(OE,pi[2][1]);
        fi; 
     else
        OE:=EquationOrder(ff);
        Omax:=MaximalOrder(OE);
        pi:=Coerce(Omax,OE.2);
        OE:=Omax;
     fi;
     i:=1;
     while i < Length(fact) do
        if Norm(fact[i][1]) = Norm(fact[i+1][1]) then
           if fact[i][2] = fact[i+1][2] then
              Remove_(fact,i);
              Remove_(fact,i);
           elif fact[i][2]>fact[i+1][2] then
              fact[i][2]:=fact[i][2]-fact[i+1][2];
              Remove_(fact,i+1);
              i:=i+1;
           else 
              fact[i+1][2]:=fact[i+1][2]-fact[i][2];
              Remove_(fact,i);
              i:=i+1;
           fi;
        else i:=i+1;
        fi;
     od;
     i:=0;
     PA2:= PolynomialAlgebra(PolynomialAlgebra(BaseRing(E1)));
     while i <Length(fact) do
        i:=i+1;
        I:=(List(Basis(fact[i][1])[1])[1])*OE+(List(Basis(fact[i][1])[2])[1]+List(Basis(fact[i][1])[2])[2]*OE.2)*OE;
        l:= Norm(I);
        r:=fact[i][2];
        no:=0;
        while no<r do
           no:=no+1;
           if l > 2 then
              ll:=List(pi-Basis(I)[2]);
              if ll[2] = 0 then
                 lambda:=Coerce(Integers(),ll[1]);
              else
                 ro:=Roots(Coerce(PolynomialAlgebra(GF(l)),CharacteristicFrobeniusPolynomial(E1)));
                 lambda:=Coerce(Integers(),ro[1][1]);
                 if not pi -lambda in I then 
                    lambda:=Coerce(Integers(),ro[2][1]);
                 fi;
              fi;
              iso1:= SubgroupFromJInvariant(E,l,1);
              if iso1[2] mod l  = lambda mod l then
                 isoli:= Append(isoli,[IsogenyFromJInvariant(E,l,1)]);
              else 
                 iso2:= IsogenyFromJInvariant(E,l,2);
                 isoli:= Append(isoli,[iso2]);
              fi;
              E:= Codomain(isoli[Length(isoli)]);
              if jInvariant(E) = j then 
                 i := Length(fact);
                 no:=r+1;
              fi;
           else 
              to2:=TorsionSubgroup(E,2);
              if Length(to2)<>2 then
                 return IsoIdeals2(Tuple([fact[i][1],r-(no-1)]),E1,E2,isoli);
              fi;
              lambda:=Coerce(Integers(),List(pi-Basis(I)[2])[1]);
              if lambda mod 2 = 0  then 
                 iso1:=FrobeniusToOtherCurve(E,2);
              else 
                 iso1:=IsogenyFromSubgroup(to2[1],E,2);
              fi;
              isoli:= Append(isoli,[iso1]);
              E:= Codomain(isoli[Length(isoli)]);
              if jInvariant(E) = j then 
                 i := Length(fact);
              fi;
           fi;
        od;
     od;
     if not E = E2 then
        if not jInvariant(E) = j then
           nr2:=true;
        else
           isoli:= Append(isoli,[Isomorphism(E,E2,1)]);
        fi;
     fi;
     Iso:= isoli[1];
     for i in [2..Length(isoli)] do
        Iso:= Composition(isoli[i],Iso);
     od;
     return Iso;
   end; 



#The composition f1(f2())
   IsoComp:= function(f2,f1) 
     local i,xelseq2,xdenom2,xnom2,xpol,re,yelseq2,ynom2,ydenom2;
     re:=[];
     xelseq2:= Eltseq(f2.xpol);
     xdenom2:= Denominator(xelseq2[1]);
     xnom2:= xdenom2*xelseq2[1];
     re[1]:= Evaluate(xnom2,f1.xpol)/Evaluate(xdenom2,f1.xpol);
     yelseq2:= Eltseq(f2.ypol);
     re[2]:= 0;
     for i in [1..Length(yelseq2)] do
        ydenom2:= Denominator(yelseq2[i]);
        ynom2:= ydenom2*yelseq2[i];
        re[2]:= (re[2]+(Evaluate(ynom2,f1.xpol)/Evaluate(ydenom2,f1.xpol)))*(f1.ypol^(i-1));
     od;
     return re;
   end;


#subfunction for composition
   _Decrease:=function(Iso,K)
     local E1,E2,E11,E22,Iso2,yp,el,K2,el1,el2;
     yp:=Iso.ypol;
     el:=Eltseq(yp);
     K2:=BaseRing(Domain(Iso));
     el1:=Coerce(PolynomialAlgebra(K2),el[1]*Denominator(el[1]));
     el2:=Coerce(PolynomialAlgebra(K2),el[2]*Denominator(el[2]));
     if not Base(IsCoercible(PolynomialAlgebra(K),el2)) or not Base(IsCoercible(PolynomialAlgebra(K),el1)) then
        return false;
     fi;
     E1:=Domain(Iso);
     E2:=Codomain(Iso);
     E11:=ChangeBaseRing(E1,K);
     E22:=ChangeBaseRing(E2,K);
     Iso2:=Copy(Iso);
     Iso2.domain:=ChangeBaseRing(Domain(Iso),K);
     Iso2.codomain:=ChangeBaseRing(Codomain(Iso),K);
     return Iso2;
   end;

#subfunction for composition
   _Increase:=function(Iso,K)
     local E1,E2,E11,E22,Iso2;
     E1:=Domain(Iso);
     E2:=Codomain(Iso);
     E11:=ChangeBaseRing(E1,K);
     E22:=ChangeBaseRing(E2,K);
     Iso2:=Copy(Iso);
     Iso2.domain:=ChangeBaseRing(Domain(Iso),K);
     Iso2.codomain:=ChangeBaseRing(Codomain(Iso),K);
     return Iso2;
   end;

#composition Iso1(Iso2)
   _Composition:= function(Iso2,Iso1)
     local  Iso,_Iso,P,EE,comp,Iso22,Iso11;
     if not Codomain(Iso1) = Domain(Iso2) then
        if IsSubset(BaseRing(Codomain(Iso1)),BaseRing(Domain(Iso2))) then
           Iso11:=Iso1;
           Iso22:=_Decrease(Iso2,BaseRing(Codomain(Iso1))); 
           if Iso22 = false then
              Iso11:=_Increase(Iso1,BaseRing(Codomain(Iso2)));
              Iso22:=Iso2;
           fi; 
        elif IsSubset(BaseRing(Domain(Iso2)),BaseRing(Codomain(Iso1))) then
           Iso11:=_Decrease(Iso1,BaseRing(Domain(Iso2)));
           Iso22:=Iso2;
           if Iso11 = false then
              Iso22:=_Increase(Iso2,BaseRing(Codomain(Iso1)));
              Iso11:=Iso1;
           fi;
        else
           return Error("composition: domain of ", Iso2," is not codomain of ",Iso1);
        fi;
        if not Codomain(Iso11) = Domain(Iso22) then
           return Error("composition: domain of ", Iso2," is not codomain of ",Iso1);
        fi;
     else 
        Iso11:=Iso1;
        Iso22:=Iso2;
     fi;
     if IsZero(Iso11) or IsZero(Iso22) then
        return NuMap(Iso11.domain,Iso22.codomain); 
     fi;
     P:= GenericPoint(Iso11.domain);
     EE:=Parent(P);
     _Iso:=function(P1)
        return Iso22(Iso11(P1));
     end;
     Iso:= Map(Iso11.domain,Iso22.codomain,_Iso);
     Iso.type:=iso;
     Iso.degree:=Iso11.degree*Iso22.degree;
     comp:= IsoComp(Iso22,Iso11);
     Iso.xpol:= comp[1];
     Iso.ypol:= comp[2];
     Iso.operations:= _iso_ops;
     if Degree(Iso1) = 1 then
        if "ei" in RecFields(Iso2) then
           Iso.ei:=Iso2.ei;
        fi;
     elif Degree(Iso2) = 1 then
        if "ei" in RecFields(Iso1) then
           Iso.ei:=Iso1.ei;
        fi;
     fi;
     return Iso;
   end;


#returns true if I is L-smooth, i.e. factors have norms given in L
   IsSmooth:= function(I,O,L,E)
     local i,m,smooth,fact,Idoben,ba,Omax;
     fact:=IdealFactorization(I);
     if fact = FAILURE then 
        return false;
     fi;
     smooth:= true;
     for i in [1..Length(fact)] do
        if not Norm (fact[i][1]) in L then
           smooth:= false;
           i:= Length(fact)+1;
        fi;
     od;
     return smooth;
   end;


   MakeSmooth:=function(I,O,L,E)
     local Omax,Cl,F,G,b,L1,bool,c,L2,M,fact,r,j,Ma,H,L3,e,II,ten,primeIds,i,l,p,ba,L3,k,pp,fact,gens,hnf,a,l,MM,m;
     if IsSmooth(I,O,L,E) then
        return I;
     fi;
     if Conductor(O)<>1 then return Error("smoothing ideals only implemented for ideals of maximal order");
     fi;
     Omax:=MaximalOrder(O);
     Cl:=ClassGroup(Omax);
     G:=Cl.ext1;
     F:=Inverse(G);
     ba:=Basis(I);
     a:=Coerce(Omax,ba[1])*Omax+(Coerce(Omax,ba[2]))*Omax;
     L1:=List(F(a));
     bool:=true;
     r:=Length(L1);
     gens:=[];
     for i in [1..r] do
        fact:=Factorization(G(Cl.(i)));
        gens[i]:=1*Omax;
        for j in [1..Length(fact)] do
           gens[i]:=gens[i]*fact[j][1];
        od;
     od;
     M:=[];
     primeIds:=[];
     c:= Conductor(EquationOrder(CharacteristicFrobeniusPolynomial(E)));
     for i in [1..Length(L)] do
        l:=L[i];
        fact:= IdealFactorization(l*Omax);
        if  GCD(Norm(fact[1][1]),c)=1 then
           Append_(primeIds,[fact[1][1]]);
           if Length(fact) = 2 then
              Append_(primeIds,[fact[2][1]]);
           fi;
        fi;
     od;
     MM:=[];
     for i in [1..r] do
        L2:=List(F(primeIds[i]));
        for j in [1..r] do
           if j=i then
              Append_(L2,[1]);
           else 
              Append_(L2,[0]);
           fi;
        od;
        MM[i]:=L2;
     od;
     M:=[];
     for i in [1..r] do
        Append_(M,MM[i]);
     od;
     Ma:=Matrix(Z,r,2*r,M);
     H:=Base(HermiteForm(Ma));
     bool:=true;
     k:=r;
     hnf:=false;
     while k<Length(primeIds) and not hnf do
        for i in [1..r] do
           if H[i][i]<>1 then
              bool:=false;
           fi;
        od;
        if not bool then
           k:=k+1;
           L2:=List(F(primeIds[k]));
           for j in [1..k-1] do
              Append_(L2,[0]);
              Append_(MM[j],[0]);
           od;
           Append_(L2,[1]);
           MM[k]:=L2;
           M:=[];
           for i in [1..k] do
              Append_(M,MM[i]);
           od;
           Ma:=Matrix(Z,k,k+r,M);
           H:=Base(HermiteForm(Ma));
           bool:=true;
        else hnf:=true;
        fi;
     od;
     L3:=[];
     for i in [1..r] do
        m:=r-i+1;
        e:=H[m][m];
        pp:=1*Omax;
        for j in [r+1..r+k] do
           pp:=pp*primeIds[j-r]^H[m][j];
        od;
        for l in [1..i-1] do
           pp:=pp*L3[l]^(-H[m][i-l]);
        od;
        pp:=pp^(-e);
        L3[i]:=Denominator(pp)*pp;
     od;        
     p:=1*Omax;
     Reverse_(L3);
     for i in [1..r] do
        e:=L1[i];
        if e <>0 then
           p:=p*L3[i]^e;
        fi;
     od;
     return p;
   end;

#old MakeSmooth, tries to make I L-smooth
   MakeSmooth2:= function(I,O,L,E)
     local Omax, fact,a,b,l,ro,f,phi,primeIds,h,ra,found,factored,pos,i,c,bli,posl,j,tr,fact,c,L2,Idsoben,ba,Iinv,m,p,n,r,k,p1,t,e,p2,e1,e2,Omax,Cl,G,F,primes;
     if IsSmooth(I,O,L,E) then
        return I;
     fi;
     c:= Conductor(EquationOrder(CharacteristicFrobeniusPolynomial(E)));
     tr:= 0;
     primeIds:=[];
     Idsoben:=[];
     Iinv:=IdInvert(I);
     Omax:= MaximalOrder(O);
     h:= PicardNumber(O);
     phi:= Basis(Omax)[2];
     found:=false;
     primes:=PrimesUpTo(1000);
     for i in [1..Length(primes)] do
        if IsElkiesPrime(E,primes[i]) then
           fact:= IdealFactorization(primes[i]*O);
           if fact<> FAILURE then
              Append_(primeIds,[fact[1][1]]);
              if Length(fact) = 2 then
                 Append_(primeIds,[fact[2][1]]);
              fi;
           fi;
        fi;
     od;
     m:=Length(primeIds);
     t:=Ceiling(Root(Root(Size(BaseRing(E)),4),m));
     a:=[];
     i:=1;
     while i <=Length(primeIds) do
        n:=Norm(primeIds[i]);
        Append_(a,[n]);
        Append_(a,[[primeIds[i],0]]);
        i:=i+1;
        if Norm(primeIds[i]) = n then
           Append_(a[Length(a)],[primeIds[i],0]);
           i:=i+1;
        fi; 
     od;
     b:= I;
     i:=0;
     while not found do
        p:=Random(primeIds);
        n:=Norm(p);
        r:=Random(t);
        b:= IdReduction(b*(p^r));
        pos:=Position(a,n);
        if a[pos+1][1] = p then
           a[pos+1][2]:=a[pos+1][2]+ r;
        else a[pos+1][4]:=a[pos+1][4]+ r;
        fi;
        fact:=IdealFactorization(b);
        i:=i+1;
        if fact <> FAILURE then
           if IsSmooth(b,O,L,E) then
              found:=true;
           fi;
        fi; 
     od;
     p:=1*O;
     k:=1;
     i:=1;
     while i <=Length(a) do
        p1:=a[i+1][1];
        e1:=a[i+1][2];
        if IsPrincipal(p1^2) then
           p2:=p1;
           e2:=e1;
        else 
           p2:=a[i+1][3];
           e2:=a[i+1][4];
        fi;
        if k<=Length(fact) then
           if p1 = fact[k][1] then
              e:=fact[k][2]-e1;
              k:=k+1;
              if e>=0 then
                 p:=p*(p1^e);
              else p:=p*p2^e;
              fi;  
           else p:=p*p2^e1;
           fi;
        else p:=p*p2^e1;
        fi;
        if k<=Length(fact) then
           if p2 = fact[k][1] then
              e:=fact[k][2]-e2;
              k:=k+1;
              if e>=0 then
                 p:=p*(p2^e);
              else p:=p*p1^e;
              fi;  
           else p:=p*p1^e2;
           fi;
        else p:=p*p1^e2;
        fi;
        i:=i+2;
     od;
     return p;
   end;  


#returns the domain of a map
   _Domain:= function(f)
     return f.domain;
   end;

#returns the codomain of a map
   _Codomain:= function(f)
     return f.codomain;
   end;




#true if there is more than one root of Psi_l(j,x)
   ExistsSecondRoot:= function(l,j)
     local ro,f,F,PA;
     F:= Parent(j);
     PA:= PolynomialAlgebra(PolynomialAlgebra(F));
     f:= Coerce(PA,ModularPolynomial(l));
     f:= Evaluate(f,j);
     ro:= Roots(f);
     if Length(ro) =1 then
        return FALSE;
     elif Length(ro) >= 2 then
        return TRUE;
     else 
        return Error("ExistsSecondRoot: l must be an Elkies Prime");
     fi;
   end;

#returns the n-th root of Psi_l(j,x)
   ModPolyRoot:= function(l,j,n)
     local ro,f,F,PA;
     F:= Parent(j);
     PA:= PolynomialAlgebra(PolynomialAlgebra(F));
     f:= Coerce(PA,ModularPolynomial(l));
     f:= Evaluate(f,j);
     ro:= Roots(f);
     if n = 1 then
        return [ro[n][1],1];
     elif n > Length(ro) then
        return [ro[1][1],1];
     else 
        return [ro[2][1],n];
     fi;
   end;


#tests if an isogeny of E is computable
   TestElkiesPrimes:= function(E)
     local el,i,j,c;
     el:= ElkiesPrimes59(E);
     j:= jInvariant(E);
     for i in [1..Length(el)] do
        if  Type(SubgroupFromJInvariant(E,el[i],1)[1])=elt-ord^rat then
           return [el[i],1];
        fi;
        if ExistsSecondRoot(el[i],j) then
           if Type(SubgroupFromJInvariant(E,el[i],2)[1]) = elt-ord^rat then
              return [el[i],2];
           fi;
        fi;
     od;
     return FAILURE;
   end;

#computes the picard-number of O
   PicardNumber:=function(O)
     local h,p,i,fact,c,Omax,UE,UO,l,pic,Id,factid;
     if IsMaximal(O) then
        return ClassNumber(O);
     fi;
     Omax:=MaximalOrder(O);
     h:=ClassNumber(Omax);
     c:=Conductor(O);
     UE:=UnitGroup(O);
     UO:=UnitGroup(Omax);
     if Size(Base(UO)) = INFTY and Size(Base(UE)) = INFTY then
        l:=1;
     else
        l:=Size(Base(UO))/Size(Base(UE));
     fi;
     pic:=1;
     fact:=Factorization(c);
     for i in [1..Length(fact)] do
        p:=fact[i][1];
        Id:=p*Omax;
        factid:=Factorization(Id);
        if Length(factid) = 2 then
           pic:=fact[i][2]*pic*(1-1/p);
        elif fact[i][2] = 1 then
           pic:=fact[i][2]*pic*(1+1/p);
        fi;
     od;
     return Coerce(Integers(),pic*h*(c/l));
   end;


#tests if an isogeny of degree l between E and E2 (with j(E) = j and j(E2) = j2) is computable 
   SubComputable:= function(j,j2,l)
     local ModPoly,BR,PA,PA2,Phix;
     ModPoly:= ModularPolynomial(l);
     BR:= Parent(j);
     PA:=PolynomialAlgebra(BR);
     PA2:= PolynomialAlgebra(PA);
     ModPoly:= Coerce(PA2,ModPoly);
     Phix:= Derivative(ModPoly);
     if Evaluate(Evaluate(Phix,j),j2) = 0 then
        return false;
     else return true;
     fi;
   end;

#tries to find a connecting pseudo-random function between E1 and E2
   PseudoRandom2:= function(E1,E2,r,L)
     local j, i, jr, k, n, j2,f,found,l,f,m,o,tes,ma,O,c,L2;
     O:=EndRing(E1);
     c:=Conductor(O);
     k:=-1;
     m:=Length(L);
     found:= false;
     while not found do
        k:= k+1;
        if k =m then 
           L2:=L;
           for i in [Last(L)..59] do
              if IsPrime(i) then
                 if IsElkiesPrime(E1,i) then
                    L2:=Append(L,[i]);
                 fi;
              fi;
           od;
           if L2 = L then
              return Error(" no connectiong pseudo random function found");
           else 
              return PseudoRandom(E1,E2,r,L2);
           fi;
        fi;
        j:= jInvariant(E1);
        for i in [1..r] do
           l:= L[(Coerce(Integers(),(j+k+1)) mod m) +1];
           n:= (Coerce(Integers(),(j+k+1)) mod 2) +1;
           j2:= ModPolyRoot(l,j,n)[1];
           f:=0;
           ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c;
           while ((not ma) or j2 = 0 or j2 = 1728 or (not SubComputable(j,j2,l))) do
              f:=f+1;
              l:= L[(Coerce(Integers(),(j+k+f+1)) mod m) +1];
              n:= (Coerce(Integers(),(j+k+1)) mod 2) +1;
              j2:= ModPolyRoot(l,j,n)[1];
              ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c;
	      if f = 5 then 
                 tes:= TestElkiesPrimes(E);
                 if tes = FAILURE then
                    return Error("Isogeny cannot be computed");
                 else j2 := ModPolyRoot(tes[1],j,tes[2])[1];
                 fi;
              fi;
           od;
           j:= j2;
        od;
        jr:=j;
        j:= jInvariant(E2);
        o:=0;
        while o < r do
           o:= o+1;
           l:= L[(Coerce(Integers(),(j+k+1)) mod m) +1];
           n:= (Coerce(Integers(),(j+k+1)) mod 2) +1;
           j2:= ModPolyRoot(l,j,n)[1];
           f:=0;
           ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c;
           while ((not ma) or j2 = 0 or j2 = 1728 or (not SubComputable(j,j2,l))) do
              f:=f+1;
              l:= L[(Coerce(Integers(),(j+k+f+1)) mod m) +1];
              n:= (Coerce(Integers(),(j+k+1)) mod 2) +1;
              j2:= ModPolyRoot(l,j,n)[1];
              ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c;
	      if f = 5 then tes:= TestElkiesPrimes(E);
                 if tes = FAILURE then
                    return Error("Isogeny cannot be computed");
                 else j2 := ModPolyRoot(tes[1],j,tes[2])[1];
                 fi;
              fi;
           od;
           j:= j2;
           if j = jr then 
              return [k,2,true,jr]; 
           fi;
        od;
     od;
   end;

#tries to find a connecting pseudo-random function between E1 and E2
   PseudoRandom59:= function(E1,E2,r,L)
     local j, i, jr, k, n, j2,f,found,l,f,m,o,tes,ma,O,c;
     k:=-1;
     O:=EndRing(E1);
     c:= Conductor(O);
     L:= ElkiesPrimes59(E1);
     m:=Length(L);
     found:= false;
     while not found do
        k:= k+1;
        if k =m then 
           return PseudoRandom2(E1,E2,r,L);
        fi;
        j:= jInvariant(E1);
        for i in [1..r] do
           l:= L[(Coerce(Integers(),(j+k+1)) mod m) +1];
           n:= (Coerce(Integers(),(j+k)) mod 2) +1;
           j2:= ModPolyRoot(l,j,n)[1];
           f:=0;
           ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c;
           while ((not ma) or j2 = 0 or j2 = 1728 or ( not SubComputable(j,j2,l))) do
              f:=f+1;
              l:= L[(Coerce(Integers(),(j+k+f+1)) mod m) +1];
              n:= (Coerce(Integers(),(j+k)) mod 2) +1;
              j2:= ModPolyRoot(l,j,n)[1];
              ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c;
	      if f = 5 then tes:= TestElkiesPrimes(E);
                 if tes = FAILURE then
                    return Error("Isogeny cannot be computed");
                 else j2 := ModPolyRoot(tes[1],j,tes[2])[1];
                 fi;
              fi;
           od;
           j:= j2;
        od;
        jr:=j;
        j:= jInvariant(E2);
        o:=0;
        while o < r do
           o:= o+1;
           l:= L[(Coerce(Integers(),(j+k+1)) mod m) +1];
           n:= (Coerce(Integers(),(j+k)) mod 2) +1;
           j2:= ModPolyRoot(l,j,n)[1];
           f:=0; 
           ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c;
           while ((not ma) or j2 = 0 or j2 = 1728 or (not SubComputable(j,j2,l))) do
              f:=f+1;
              l:= L[(Coerce(Integers(),(j+k+f+1)) mod m) +1];
              n:= (Coerce(Integers(),(j+k)) mod 2) +1;
              j2:= ModPolyRoot(l,j,n)[1];
              ma:=Conductor(EndRing(EllipticCurveFromJInvariant(j2)))=c;
	      if f = 5 then 
                 tes:= TestElkiesPrimes(E);
                 if tes = FAILURE then
                    return Error("Isogeny cannot be computed");
                 else j2 := ModPolyRoot(tes[1],j,tes[2])[1];
                 fi;
              fi;
           od;
           j:= j2;
           if j = jr then 
              return [k,1,true,jr]; 
           fi;
        od;
     od;
   end;

#tries to find a connecting pseudo-random function between E1 and E2
   PseudoRandom:= function(E1,E2,r,L,polys)
     local j, i, jr, k, n,loop,c,xalt, j2,f,found,l,f,m,o,tes,ma,O,c,eval,BR,q,p,L2,poss,mi,ind,PA2,r2,degs,els;
     m:=Length(L);
     BR:= BaseRing(E1);
     found:= false;
     k:=-1;
     poss:=[];
     loop:=Last(L) > 37;
     eval:=function(l,j,n)
        local ro;
        ro:= Roots(Evaluate(polys[Position(L,l)],j));
        if n = 1 then
           return [ro[n][1],1,Length(ro)];
        elif (n =2 and Length(ro) =1) then
           return [ro[1][1],1,Length(ro)];
        else return [ro[2][1],2,Length(ro)];
        fi;
     end;
     q:=1;
     Print("von E1\n");
     while not found do
        k:= k+1;
        if k =m then 
           if poss=[] then
              L2:=L;
              PA2:=Parent(polys[1]);
              for i in [Last(L)+1..59] do
                 if IsPrime(i) then
                    if IsElkiesPrime(E1,i) then
                       Append_(L2,[i]);
                       Append_(polys,[Coerce(PA2,ModularPolynomial(i))]);
                    fi;
                 fi;
              od;
              if Length(L2) = m then 
                 return Error(" no connectiong pseudo random function found");
              else return PseudoRandom(E1,E2,r,L2,polys);
              fi;
           elif Length(poss)= 1 then
              return poss[1][1];
           else mi:=poss[1][3];
              ind:=1;
              for i in [2..Length(poss)] do
                 if poss[i][3] < mi then
                    mi:=poss[i][3];
                    ind:=i;
                 fi;
              od;
              return poss[ind][1];
           fi;
        fi;        
        j:= jInvariant(E1);
        i :=0;
        while i <r do
           i:=i+1;
           els:=Eltseq(j);
           q:=Sum(function(x) return Coerce(Integers(),x); end,els);
           l:= L[(q+k+1) mod m +1];
           n:= ((q+k) mod 2) +1;
           j2:= eval(l,j,n);
           f:=0;
           while (j2[1] = 0 or j2[1] = 1728 or (not SubComputable(j,j2[1],l)) or j2[1] = j) do
              f:=f+1;
              l:= L[(q+k+f+1) mod m +1];
              n:= ((q+k) mod 2) +1;
              j2:= eval(l,j,n);
	      if f = 5 then tes:= TestElkiesPrimes(E1);
                 if tes = FAILURE then
                    return Error("Isogeny cannot be computed");
                 else j2 := eval(tes[1],j,tes[2]);
                 fi;
              fi;
           od;
           j:= j2[1];
           Print(j,"\n");
        od;
        jr:=j;
        j:= jInvariant(E2);
        Print("von E2\n");
        o:=0;
        degs:=0;
        while o < r do
           o:= o+1;
           els:=Eltseq(j);
           q:=Sum(function(x) return Coerce(Integers(),x); end,els);
           l:= L[(q+k+1) mod m +1];
           n:= ((q+k) mod 2) +1;
           j2:= eval(l,j,n);
           f:=0;
           while (j2[1] = 0 or j2[1] = 1728 or not (SubComputable(j,j2[1],l))or j2[1]=j) do
              f:=f+1;
              l:= L[(q+f+k+1) mod m +1];
              n:= ((q+k) mod 2) +1;
              j2:= eval(l,j,n);
	      if f = 5 then tes:= TestElkiesPrimes(E1);
                 if tes = FAILURE then
                    return Error("Isogeny cannot be computed");
                 else j2 := eval(tes[1],j,tes[2]);
                 fi;
              fi;
           od;
           degs:=degs+l;
           j:= j2[1];
           Print(j,"\n");
           if j = jr then 
              return [k,1,loop,jr];
           fi;
        od;
     od;
   end;



#returns ideal that represents an isogeny between E1 and E2
   GetIsoIdeal:= function(E1,E2,c)
     local pos,p,r,k,I1,I2,n,E,I,f,h,BR,EC,fro,P,mu,ro,frob,R,lambda,sub,nu,su,tes,nr,eq,D,a,fact,pos,ma,PA2,eval,mopos,index,co,Etmp,m,cha2,Iso,l,el,Isos,Ids,dom,Id,Cl,O,Ge,Gen,ma,Id,pi1,pi,su,j,i,Id,ord,o,elli,pos,j1,j2,jInvs,found,jr,ch,id,PA,phi,c,Omax1,Omax,polys,pri,q,jr,L,mytime,loop,OE,w,r1,r2,ba,ff,rlist,xalt,els,mytime;
     BR:= BaseRing(E1);
     pri:=PrimitiveElement(BR);
     PA2:=PolynomialAlgebra(PolynomialAlgebra(BR));
     ff:=CharacteristicFrobeniusPolynomial(E1);
     if c > 1 then
        pi1:=EquationOrder(ff).2;
        Omax1:=MaximalOrder(Parent(pi1));
        OE:=EquationOrder(MinimalPolynomial(c*Omax1.2));
        Omax:=MaximalOrder(OE);
        ff:=Coerce(PolynomialAlgebra(Omax),ff);
        pi:=Roots(ff);
        rlist:=List(pi[1][1]);
        if rlist[2]>0 then
           pi:=Coerce(OE,pi[1][1]);
        else 
           pi:=Coerce(OE,pi[2][1]);
        fi; 
     else
        OE:=EquationOrder(ff);
        Omax:=MaximalOrder(OE);
        pi:=Coerce(Omax,OE.2);
        OE:=Omax;
     fi;
     phi:=Omax.2;
     PA:=PolynomialAlgebra(Omax);
     ch:= CharacteristicFrobeniusPolynomial(E1);
     el:= ElkiesPrimes(E1);
     j:= jInvariant(E1);
     m:=Length(el);
     h:= PicardNumber(OE);
     r:= Ceiling(Sqrt(h));
     k:= -1;
     found:=false;
     I1:=1*OE;
     I2:=1*OE;
     E:=E1;
     m:=Length(el);
     found:= false;
     polys:=[];
     for i in [1..Length(el)] do
        polys[i]:= Coerce(PA2,ModularPolynomial(el[i]));
     od;
     k:=PseudoRandom(E1,E2,r,el,polys);
     nr:= k[2]-1;
     jr:=k[4];
     if k[3] then el:=ElkiesPrimes59(E1);
        polys:=[];
        for i in [1..Length(el)] do
           polys[i]:= Coerce(PA2,ModularPolynomial(el[i]));
        od; 
     fi;
     k:=k[1];
     m:=Length(el);
     eval:=function(l,j,n)
        local ro;
        ro:= Roots(Evaluate(polys[Position(el,l)],j));
        if n = 1 then
           return [ro[n][1],1,Length(ro)];
        elif n > Length(ro) then
           return [ro[1][1],1,Length(ro)];
        else return [ro[n][1],n,Length(ro)];
        fi;
     end;
     Print("computing subgroups\n");
     while not found do
        E:=E1;
        j:=jInvariant(E1);
        i:=0;
        while i <r do
           i:=i+1;
           els:=Eltseq(j);
           q:=Sum(function(x) return Coerce(Integers(),x); end,els);
           l:= el[(q+k+1) mod m +1];
           n:= (Coerce(Integers(),(q+k)) mod 2) +1;
           j2:= eval(l,j,n);
           su:= (SubgroupFromJInvariant(E,l,j2[2]));
           f:=0;
           while (j2[1] = 0 or j2[1] = 1728 or not Type(su) = list or j2[1]=j) do
              f:=f+1;
              l:= el[(q+f+k+1) mod m +1];
              n:= (Coerce(Integers(),(q+k)) mod 2) +1;
              j2:= eval(l,j,n);
              su:= SubgroupFromJInvariant(E,l,j2[2]);
              if f = 5 then tes:= TestElkiesPrimes(E);
                 if tes = FAILURE then
                    return Error("Isogeny cannot be computed");
                 else su:=SubgroupFromJInvariant(E,tes[1],tes[2]);
                 fi;
              fi;
           od;
           lambda:=su[2];
           Id:=l*OE+(pi-lambda)*OE;
           I1:= IdReduction(Id*I1);
           j:=j2[1];
           Print("j: ",j, "l: ",l,"\n");
           E:=ECFromJGivenNumber(j,Size(E));
           if j = jr then
              i:=r;
           fi;
        od;
        Print("von E2\n");
        jr:=j;
        j:= jInvariant(E2);
        E:=E2;
        nu:=0;
        while nu <r do
           nu:=nu+1;
           els:=Eltseq(j);
           q:=Sum(function(x) return Coerce(Integers(),x); end,els);
           l:= el[(q+k+1) mod m +1];
           n:= ((q+k) mod 2) +1;
           j2:= eval(l,j,n); 
           su:= (SubgroupFromJInvariant(E,l,j2[2]));
           f:=0;
           while j2[1] =0 or j2[1] = 1728 or Type(su)<>list or j2[1]=j do 
              f:=f+1;
              l:= el[(q+f+k+1) mod m +1];
              n:= ((q+k) mod 2) +1;
              j2:= eval(l,j,n); 
              su:= SubgroupFromJInvariant(E,l,j2[2]);
              if f = 5 then 
                 tes:= TestElkiesPrimes(E);
                 if tes = FAILURE then 
                    return Error("Isogeny cannot be computed");
                 else su:=SubgroupFromJInvariant(E,tes[1],tes[2]);
                 fi;
              fi;
           od;
           Id:=l*OE+(pi-su[2])*OE;
           I1:=IdReduction(IdInvert(Id)*I1);
           j:=j2[1];
           Print("j: ",j, "l: ",l,"\n");
           E:=ECFromJGivenNumber(j,Size(E));
           if j = jr then 
              found:=true; 
              nu := r+1; 
           fi;
        od;
     od;
     Print(I1);
     return MakeSmooth(I1,OE,el,E1);
   end;




#converts a hex-number to an integer
   HexToInt:=function(L)
     local i,sum,k,tmp;
     sum:=0;
     k:=Length(L);
     for i in [1..k] do
        tmp:=L[k-i+1];
        if not (ASCII(tmp)-48) in [0..9] then
           tmp:=ASCII(tmp)-87;
        else tmp:=ASCII(tmp)-48;
        fi;
        sum:=sum+tmp*16^(i-1);
     od;
     return sum;
   end;







#factorization of an ideal prime to the conductor of its order
   IdealFactorization:=function(Id)
     local O,Omax,fact,Ioben,ba,i,fact2,II,ten,d,Id2;
     O:=Order(Id);
     d:=Denominator(Id);
     Assert(d = 1," integral ideal ");
     if IsMaximal(O) then
        return Factorization(Id);
     else 
        Omax:=MaximalOrder(O);
        ba:=Basis(Id);
        Ioben:=Coerce(Omax,ba[1])*Omax+(Coerce(Omax,ba[2]))*Omax;
        fact:=Factorization(Ioben);
        fact2:=[];
        for i in [1..Length(fact)] do
           II:=Intersection(O,fact[i][1]);
           if GCD(Norm(II),Conductor(O))<>1 then
              return FAILURE;
           fi;
           ten:=TwoElementNormal(II);
           II:=Ideal(O,[Base(ten),ten.ext1]);
           II:=Ideal(O,Basis(II));
           Append_(fact2,[Tuple([II,fact[i][2]])]);
        od;
     fi;
     return fact2;
   end;
                
  
#subfunction of FindIso
   _FindIso:= function(E1,E2,c)
     local O,Id,fact;
     if not BaseRing(E1) = BaseRing(E2) then
        return Error("FindIso: given curves have different base rings");
     elif not Size(E1) = Size(E2) then
        return Error("FindIso:= given curves are not isogenous");
     elif IsSupersingular(E1) or IsSupersingular(E2) then
        return Error("Algorithm does not work for supersingular curves");
     elif jInvariant(E1) = jInvariant(E2) then
        return Isomorphism(E1,E2,1);
     fi;
     Id:= GetIsoIdeal(E1,E2,c);
     if GCD(Norm(Id),c)<>1 then
        Id:=FindIso(E1,E2,true);
     fi;
     fact:= IdealFactorization(Id);
     if fact = Sequence([]) then
        return Error("something went wrong: Ideal representation of isogeny is ",Id);
     fi;
     return IsoFromIdeals(fact,E1,E2);
   end;



#tries to find an isogeny between E1 and E2
   FindIso:= function(E1,E2) 
     local O,EE1,EE2, PA,BR,E3,E4,c1,c2,OE1,OE2,fact1,fact2,len,i,k,j,Psi,ro,l,E,Iso,Iso1,Iso2,Iso11,Iso21,n,Otmp,r,c,s,kgV,mytime;
     mytime:=Realtime();
     if IsSupersingular(E1) or IsSupersingular(E2) then
        return Error("Algorithm does not work for supersingular curves");
     fi;
     if BaseRing(E1) <> BaseRing(E2) then
        return Error("Curves are defined over different Base Rings");
     fi;
     if Characteristic(BaseRing(E1))>3 then
        EE1:=ShortWNF(E1);
        EE2:=ShortWNF(E2);
     elif Characteristic(BaseRing(E1)) =2 then
        EE1:=NormalFormChar2(E1,1);
        EE2:=NormalFormChar2(E2,1);
     else return Error("Algorithm not implemented for Characteristic 3");
     fi;
     E1:=Base(EE1);
     E2:=Base(EE2);
     if EndRingIsMaximal(E1) and EndRingIsMaximal(E2) then        
        Iso:= _FindIso(E1,E2,1);
        return Composition(EE2.ext1,Composition(Iso,EE1.ext1));
     fi;
     OE1:= EndRing(E1);
     if OE1 = FAILURE then 
        return Error("cannot compute endomorphism ring: conductor too large"); 
     fi;
     OE2:= EndRing(E2);
     if OE2 = FAILURE then 
        return Error("cannot compute endomorphism ring: conductor too large");
     fi;
     c1:= Conductor(OE1);
     c2:= Conductor(OE2);
     kgV:=Lcm(c1,c2);
     fact1:=Factorization(Coerce(Integers(),kgV/c1));  
     fact2:=Factorization(Coerce(Integers(),kgV/c2));  
     Iso1:=[];
     Iso2:=[];
     BR:= BaseRing(E1);
     PA:=PolynomialAlgebra(PolynomialAlgebra(BR));
     j:= jInvariant(E1);
     E:= E1;
     O:=OE1;
     s:=Size(E1);
     if Length(fact1)>0 then
        for i in [1..Length(fact1)] do
           l:= fact1[i][1];
           Psi:= Coerce(PA,ModularPolynomial(l));
           for k in [1..fact1[i][2]] do
              ro:= Roots(Evaluate(Psi,j));
              if Length(ro) = 1 then
                 Iso:=IsogenyFromJInvariant(E,l,1);
              elif Length(ro) = l+1 then 
                 n:=0;
                 while n in [0..l] do
                    n:=n+1;
                    if ro[n][1] <>0 then
                       Otmp:= EndRing(ECFromJGivenNumber(ro[n][1],s));
                       if Conductor(Otmp)< Conductor(O) then
                          Iso:= IsogenyFromJInvariant(E,l,n);
                          if (not Iso = FAILURE) then
                             n:= l+3;
                          fi;
                       fi;
                    fi; 
                 od;
                 if n = l+1 then 
                    return Error("no suitable curve found");
                 fi;
              fi;
              if Iso= FAILURE then 
                 return Error("cannot compute isogenous maximal endRing curve for E1");
              fi;
              Iso1:= Append(Iso1,[Iso]);
              E:= Codomain(Iso);
              j:= jInvariant(E);
              O:= EndRing(E);
           od;
        od;
        Iso11:= Iso1[1];
        for i in [2..Length(Iso1)] do
           Iso11:= Composition(Iso1[i],Iso11);
        od;
     else Iso11:= Isomorphism(E1,E1,1);
     fi;
     if Length(fact2)>0 then
        E:= E2;
        j:= jInvariant(E2);
        O:=OE2;
        for i in [1..Length(fact2)] do
           l:= fact2[i][1];
           Psi:= Coerce(PA,ModularPolynomial(l));
           for k in [1..fact2[i][2]] do
              ro:= Roots(Evaluate(Psi,j));
              if Length(ro) = 1 then
                 Iso:=IsogenyFromJInvariant(E,l,1);
              elif Length(ro) = l+1 then 
                 n:=0;
                 while n in [0..l] do
                    n:=n+1;
                    if ro[n][1] <>0 then
                       Otmp:= EndRing(ECFromJGivenNumber(ro[n][1],s));
                       if Conductor(Otmp)< Conductor(O) then
                          Iso:= IsogenyFromJInvariant(E,l,n);
                          if not Iso= FAILURE then
                             n:= l+3;
                          fi;
                       fi;
                    fi;
                 od;
                 if n = l+1 then 
                    return Error("something went wrong 2");
                 fi;
              fi;
              if Iso= FAILURE then 
                 return Error("cannot compute isogenous maximal endRing curve for E2");
              fi;
              Iso2:= Append(Iso2,[Iso]);
              E:= Codomain(Iso);
              j:= jInvariant(E);
              O:=EndRing(E);
           od;
        od;
        len:= Length(Iso2);
        Iso21:= DualIsogeny(Iso2[len]);
        for i in [1..len-1] do
           Iso21:= Composition(DualIsogeny(Iso2[len-i]),Iso21);
        od;
     else Iso21:= Isomorphism(E2,E2,1);
     fi;
     E3:= Codomain(Iso11);
     E4:= Domain(Iso21);
     c:=Conductor(EndRing(E3));
     Print("init:",Realtime()-mytime,"\n");
     Iso:= _FindIso(E3,E4,c);
     Iso:= Composition(Iso,Iso11);
     Iso:= Composition(Iso21,Iso);
     return Composition(EE2.ext1,(Composition(Iso,EE1.ext1)));
   end;


#direkte magma-Uebersetzung

   WeierstrassFunction:= function(E,prec)
     local Q,z,coeffs,a4,a6,k,c,ctmp,h,ztmp;
     if (2*prec+1) >= Characteristic(BaseRing(E)) then return Error(
        "Request for division by zero in kernel polynomial computation.");
     fi;
     Q := LaurentSeriesRing(BaseRing(E),2*prec+2);
     z := Q.1;
     coeffs:=Coefficients(E);
     a4:=coeffs[1];
     a6:=coeffs[2];
     c := [ -a4/5, -a6/7 ];
     for k in [3..prec-1] do
        ctmp:=0;
        for h in [1..k-2] do
           ctmp:=ctmp+c[h]*c[k-1-h];
        od;
        Append_(c, [3*ctmp/ ((k-2) * (2*k + 3))]);
     od;
     ztmp:=0;
     for k in [1..prec-1] do
        ztmp:= ztmp+c[k]*z^(2*k);
     od;
     return z^-2 +ztmp;
   end;


   KernelPolynomial:=function(E,l,Eb,p1)
     local FF,xtr,deg,coeffs,coeffs2,e4,e6,e4_tilde,e6_tilde,c,c_tilde,P,x,psi2,dpsi2,mulist,k,pp,np,ii,i,M,B,v_raw,V,p0,v,ps,tlist,g,tmp;
     FF := BaseRing(E);
     xtr := 2;
     deg := Div(l-1, 2);
     if (2*(deg+xtr)+3) >= Characteristic(FF) then
        return Error(
        "Request for division by zero in kernel polynomial computation.");
     fi;
     coeffs:=Coefficients(E);
     coeffs2:=Coefficients(Eb);
     if Length(coeffs) <> 2 then 
        return Error;
     fi;
     e4:=coeffs[1];
     e6:=coeffs[2];
     e4_tilde:=coeffs2[1];
     e6_tilde := coeffs2[2];
     c := WeierstrassFunction(E,deg+xtr+1);
     c_tilde := WeierstrassFunction(Eb,deg+xtr+1);
     P := PolynomialRing(FF); x := P.1;
     psi2 := 4*x^3 + 4*e4*x + 4*e6;
     dpsi2 := 6*x^2 + 2*e4;
     mulist := [ 2*e4 + 6*x^2 ];
     for k in [2..deg+xtr] do
        pp := Eltseq(mulist[k-1]);
        np := pp[2] * dpsi2;
        for ii in [3..Length(pp)] do
           i := ii - 1;
           np := np+ i*pp[i+1]*(dpsi2*x^(i-1) + (i-1)*psi2*x^(i-2));
        od;
        Append_(mulist,[ np]);
     od;
     M:=[];
     for k in [1..deg+xtr] do
        pp := Eltseq(mulist[k]);
        for i in [1..Length(pp)] do
           M[(k-1)*(deg+2+xtr)+i] := pp[i] * 2 / Factorial(2*k);
        od;
        for i in [Length(pp)+1..deg+2+xtr] do
           M[(k-1)*(deg+2+xtr)+i]:=0;
        od;
     od;
     M :=Coerce( KMatrixSpace(FF, deg+xtr, deg+2+xtr),M);
     B:=[];
     for i in [1..deg+xtr] do
        B[i]:=Coefficient(c_tilde,2*i) -
        Coefficient(c,2*i);
     od;
     B:=Vector(FF,B);
     v_raw:= Solution(Transpose(M), B);
     V:=v_raw.ext1;
     if Dimension(V) <> 2 then
        return Error( "Trace term not uniquely determined in kernel polynomial computation.");
     fi;
     p0 := deg;
     v :=  v_raw + V.1 * (p0 - Eltseq(v_raw)[1]);
     v :=  v + V.2 * (p1 - Eltseq(v)[2]);
     ps := Eltseq(v);
     tlist:=[];
     for i in [1..deg+1+xtr] do
        tlist[i]:=Coerce(FF,0);
     od;
     tlist[deg+1+xtr] := Coerce(FF,1);
     for k in [1..deg+xtr] do
        tmp := Coerce(FF,0);
        for i in [1..k-1] do
           tmp := tmp + tlist[deg+xtr-i+1] * ps[k-i+1];
        od;
        tmp := -tmp - ps[k+1];
        tlist[deg+xtr-k+1] := tmp / k;
     od;
     g := Coerce(P,Sublist(tlist,[xtr+1..deg+xtr+1]));
     return g;
end;


 

   PolynomialExpression:=function(S2,S1)
     local SS,t,F,PP,x,vS1,vS2,degP,S1_powers_list,S1_i,i,powers_of_S1,P,S2tmp,S1coeff,S2coeff,Pcoeff;
     SS := Parent(S1);
     t := SS.1;
     F := BaseRing(Parent(S1));
     PP := PolynomialRing(F);
     x := PP.1; 
     vS1 := Valuation(S1);
     vS2 := Valuation(S2);
     degP := Div(vS2,vS1);
     if vS1 = 0 or vS2 = 0 or (vS2 mod vS1) <> 0 or degP < 0 then
        return Error( "Invalid data in polynomial reconstruction.");
     fi;
     S1_powers_list := [Coerce(Parent(S1),1)];
     S1_i := Coerce(Parent(S1),1);
     for i in [1..Div(vS2,vS1) ] do
        S1_i := S1_i * S1;
	Append_(S1_powers_list, [S1_i]);
     od;
     powers_of_S1 := function(i) 
        return S1_powers_list[i+1];
     end;
     P := Coerce(PP,0);
     S2tmp := S2;
     S1coeff := Coefficient(S1, vS1);
     for i in [0..degP] do   
        S2coeff := Coefficient(S2tmp, vS2 - i*vS1);
	Pcoeff := S2coeff / S1coeff;
	P := P + Pcoeff * x^(degP - i);
        S2tmp := S2tmp - Pcoeff * powers_of_S1(degP - i);
	if (i <> degP) and S2tmp = Coerce(Parent(S2tmp),0) then
	   Print("not enough terms to determine poly");
	   return Coerce(PP,0);
	fi;
     od;
     return P;
   end;

   IsogenyPhi:=function(E,F,psi,prec)
     local wE,wF,den,nu;
     wE := WeierstrassFunction(E,prec);
     wF := WeierstrassFunction(F,prec);
     den := Evaluate(psi,wE);
     nu := den^2*wF;
     return PolynomialExpression(nu,wE);
   end;


   IsogenyFromKernel:=function(E,F,psi)
     local wE,wF,den,nu,phi,wF2,wE2,omega,xmap,ymap,_Iso,re,prec,deg;
     deg:=Degree(psi);
     prec:=3*deg;
     wE := WeierstrassFunction(E,prec);
     wF := WeierstrassFunction(F,prec);
     den := Evaluate(psi,wE);
     nu := den^2*wF;
     phi:= PolynomialExpression(nu,wE);
     wF2:=Derivative(wF);
     wE2:=Derivative(wE);
     den:=Evaluate(psi,wE);
     nu:=den^3*wF2;
     omega:=PolynomialExpression(nu/wE2,wE);
     xmap:=function(P)
        local x;
        x:=PX(P);
        if x = INFTY then
           return INFTY;
        elif Evaluate(psi,x) = 0 then
           return INFTY;
        else return Evaluate(phi/psi^2,x);
        fi;
     end;
     ymap:=function(P)
        local x,y;
        x:=PX(P);
        y:=PY(P);
        return Evaluate(omega/psi^3,x)*y;
     end;
     _Iso:=function(P)
        if xmap(P) = INFTY then
           return Zero(F);
        else return MakePoint(F,[xmap(P),ymap(P)]);
        fi;
     end;
     re:=Map(E,F,_Iso);
     re.operations:= _iso_ops;
     re.type:= iso;
     re.ypol:=omega/psi^3*EllipticFunctionField(E).1;
     re.xpol:=Coerce(Parent(re.ypol),phi/psi^2);
     re.degree:=2*deg+1;
     F.size:=E.size;
     return re;
   end;






InstallMethod(
rec(
    kind:="FUNCTION",
    name:="Degree",
    sin:=[[iso,"f"]],
    sou:=[[elt-ord^rat,"d"]],
    short:=
      "The degree of f ",
    ex:=[],
    see:=[]), _IsoDegree);

InstallMethod(
rec(
    kind:="FUNCTION",
    name:="Composition",
    sin:=[[iso,"phi"],[iso,"psi"]],
    sou:=[[iso]],
    short:=
      "The composition phi(psi()). ",
    ex:=[],
    see:=[]), _Composition);
InstallMethod(
rec(
    kind:="FUNCTION",
    name:="Domain",
    sin:=[[iso,"phi"]],
    sou:=[[elt-grp^ell]],
    short:=
      "The domain of phi. ",
    ex:=[],
    see:=[]), _Domain);


InstallMethod(
rec(
    kind:="FUNCTION",
    name:="Codomain",
    sin:=[[iso,"phi"]],
    sou:=[[elt-grp^ell]],
    short:=
      "The codomain of phi. ",
    ex:=[],
    see:=[]), _Codomain);

InstallMethod(
rec(
    kind:="FUNCTION",
    name:="Conductor",
    sin:=[[ord^num,"O"]],
    sou:=[[elt-ord^rat]],
    short:=
      "The index of O in its maximal order. ",
    ex:=[],
    see:=[]), _IntCond);

