create_primes:=function(p,i)							/* This function creates the list of primes such that the product of these primes exeeds 4*Sqrt(p^i) */
        if p eq 3 then			
                list:=[5];
        else
                list:=[3];
        end if;
        product:=list[1];
        while product le 4*Sqrt(p^i) do
                s:=NextPrime(list[#list]);
                if not s eq p then
                        Append(~list,s);
                else
                        Append(~list,NextPrime(s));
                end if;
                product := product*list[#list];
        end while;
        return product,list;
end function;

s:=function(n)									/* This function is made such that f[s(n)] is actually the n'th reduced division polynomial */
        return n+2;
end function;

create_polynomials:=function(number,x,A,B)					/* This function creates a list of the first `number' reduced division polynomials */
        f:=[-1,0,1,2,3*x^4+6*A*x^2+12*B*x-A^2,4*(x^6+5*A*x^4+20*B*x^3-5*A^2*x^2-4*A*B*x-8*B^2-A^3)];
        k:=5;
        while k le number do
                n:= k div 2;
                if IsEven(k) then
                        Append(~f,1/2*f[s(n)]*(f[s(n+2)]*f[s(n-1)]^2-f[s(n-2)]*f[s(n+1)]^2));
                else
                        if IsEven(n) then
                                Append(~f,f[s(n+2)]*f[s(n)]^3*(x^3+A*x+B)^2-f[s(n+1)]^3*f[s(n-1)]);
                        else
                                Append(~f,f[s(n+2)]*f[s(n)]^3-f[s(n+1)]^3*f[s(n-1)]*(x^3+A*x+B)^2);
                        end if;
                end if;
                k:=k+1;
        end while;
        return f;
end function;

case_1:=function(f,l,q,x,A,B)  							/* Case 1 is the case where phi^2(P)=+/-[q]P for some point P in E[l] */
        leg:=LegendreSymbol(q,l);
        k:=q mod l;
        if leg eq -1 then
                return 0;
        end if;
        w:=1;
        while not (w^2 mod l) eq k do
                w:=w+1;
        end while;
        P<y>:=PolynomialRing(GF(q));
        Q<y>:=quo<P|P!f[s(l)]>;							/* This is the polynomial ring over GF(q) modulo f_l */
        gw:=[Q!(P!f[s(w-2)]),Q!(P!f[s(w-1)]),Q!(P!f[s(w)]),Q!(P!f[s(w+1)]),Q!(P!f[s(w+2)])]; 			/* Here the necessary reduced division polynomials are computed modulo f_l */
        if IsEven(w) then
                Pol:=(y^q-y)*gw[3]^2*(y^3+A*y+B)+gw[2]*gw[4];
        else
                Pol:=(y^q-y)*gw[3]^2+gw[2]*gw[4]*(y^3+A*y+B);
        end if;
        gcd:=Gcd(Parent(x)!Pol,f[s(l)]);					/* This is the first gcd that has to be computed */
        if gcd eq 1 then
                return 0;
        end if;
        if IsEven(w) then
                Pol:=4*(y^3+A*y+B)^(Integers()!((q+3)/2))*gw[3]^3-gw[5]*gw[2]^2+gw[1]*gw[4]^2;
        else
                Pol:=4*(y^3+A*y+B)^(Integers()!((q-1)/2))*gw[3]^3-gw[5]*gw[2]^2+gw[1]*gw[4]^2;
        end if;
        gcd:=Gcd(Parent(x)!Pol,f[s(l)]);					/* This is the second gcd that has to be computed */	
        if gcd eq 1 then
                return -2*w mod l;
        else
                return 2*w mod l;
        end if;
end function;

case_2:=function(f_old,l,q,x,A,B)						/* Case 2 is the case where phi^2(P) is not equal to +/-[q]P for all points P in E[l] */
        k:= q mod l;
        x:=quo<Parent(x)|f_old[s(l)]>!x;
        f:=[];
        for i in [1..l+3] do
                Append(~f,Parent(x)!(f_old[i]));
        end for;
        S:=x^3+A*x+B;
        Sq_half:=S^(Integers()!((q-1)/2));
        Sq:=Sq_half^2*S;                                			/*This is S^q*/
        Sq_square_half:=Parent(S)!(Sq_half^(q+1));      			/*This is S^((q^2-1)/2)*/
        if IsEven(k) then
                alpha:=f[s(k+2)]*f[s(k-1)]^2-f[s(k-2)]*f[s(k+1)]^2-4*Sq_square_half*S^2*f[s(k)]^3;
                beta:=4*S*f[s(k)]*((x-x^q^2)*S*f[s(k)]^2-f[s(k-1)]*f[s(k+1)]);
                delta1:=beta^2*(f[s(k-1)]*f[s(k+1)]-(x^q^2+x^q+x)*S*f[s(k)]^2)+alpha^2*S^2*f[s(k)]^2;
                delta3:=4*Sq_half*S*(alpha*beta^2*(S*f[s(k)]^2*(2*x^q^2+x)-f[s(k-1)]*f[s(k+1)])-S^2*f[s(k)]^2*(alpha^3+Parent(S)!(Sq_square_half/S)*beta^3));
        else
                alpha:=S*(f[s(k+2)]*f[s(k-1)]^2-f[s(k-2)]*f[s(k+1)]^2)-4*Sq_square_half*S*f[s(k)]^3;
                beta:=4*f[s(k)]*((x-x^q^2)*f[s(k)]^2-S*f[s(k-1)]*f[s(k+1)]);
                delta1:=beta^2*S*(S*f[s(k-1)]*f[s(k+1)]-(x^q^2+x^q+x)*f[s(k)]^2)+alpha^2*f[s(k)]^2;
                delta3:=4*Sq_half*(S*alpha*beta^2*(f[s(k)]^2*(2*x^q^2+x)-S*f[s(k-1)]*f[s(k+1)])-f[s(k)]^2*(alpha^3+Sq_square_half*S^2*beta^3));
        end if;
        delta2:=beta^2*S*f[s(k)]^2;
	delta4:=beta*delta2;
        fq:=[];
        for i in [1..4] do
                Append(~fq,f[i]^q);
        end for;
        for t in [1..l-1] do
		Append(~fq,f[t+4]^q);
                if IsEven(t) then						/* Here the first polynomial is computed */
                        Pol1:=delta1*fq[s(t)]^2*Sq+delta2*fq[s(t-1)]*fq[s(t+1)];
                else
                        Pol1:=delta1*fq[s(t)]^2+delta2*fq[s(t-1)]*fq[s(t+1)]*Sq;
                end if;
		if Pol1 eq 0 then
			if IsEven(t) then					/* Here the second polynomial is computed */
				Pol2:=delta3*fq[s(t)]^3*S^(Integers()!((3*q+1)/2))-delta4*(fq[s(t+2)]*fq[s(t-1)]^2-fq[s(t-2)]*fq[s(t+1)]^2)*S^(Integers()!((q+1)/2));
			else
				Pol2:=delta3*fq[s(t)]^3-delta4*(fq[s(t+2)]*fq[s(t-1)]^2-fq[s(t-2)]*fq[s(t+1)]^2)*Sq;
			end if;
                	if Pol2 eq 0 then
                        	return t;
                	end if;
		end if;
        end for;
        return "No trace has been found";
end function;

calculate_trace_modl:=function(f,l,q,x,A,B)					/* This function computes the trace mod l for a given prime number l */
        k:=q mod l;
        P<y>:=PolynomialRing(GF(q));
        Q<y>:=quo<P|P!f[s(l)]>;
        gk:=[Q!(P!f[s(k-1)]),Q!(P!f[s(k)]),Q!(P!f[s(k+1)])];
        if IsEven(k) then
                Pol:=(y^q^2-y)*gk[2]^2*(y^3+A*y+B)+gk[1]*gk[3];
        else
                Pol:=(y^q^2-y)*gk[2]^2+gk[1]*gk[3]*(y^3+A*y+B);
        end if;
        Pol:=Parent(x)!Pol;
        gcd:=Gcd(Pol,f[s(l)]);
        if gcd eq 1 then
                return case_2(f,l,q,x,A,B);
        else
                return case_1(f,l,q,x,A,B);
        end if;
end function;

main:=function(p,i,A,B)	
        product,primes:=create_primes(p,i);
        F:=GF(p^i);
        if F!(4*A^3+27*B^2) eq 0 then
                return "The curve is singular.";
        end if;
        P<x>:=PolynomialRing(AlgebraicClosure(F));
        f:=create_polynomials(primes[#primes]+1,x,A,B);
        traces:=[];
        for l in primes do
                Append(~traces,calculate_trace_modl(f,l,p^i,x,A,B));
        end for;
        trace_mod:=ChineseRemainderTheorem(traces,primes);
        if trace_mod lt 2*Sqrt(p^i) then
                trace:=trace_mod;
        else
                trace:=trace_mod-product;
        end if;
        return p^i+1-trace;
end function;