////////////////////////////////////////////////////////////////////////////////////////////////
// Auxiliary functions:
////////////////////////////////////////////////////////////////////////////////////////////////

function CM_EllipticCurve(Df2,d)
// Input: Df2=discriminant of the CM quadratic order and d=twist an integer
// return the d-twist of the elliptic curve E_cm
  case Df2: 
    when -3:   s:=[ 0, 1 ]; E:=EllipticCurve([0,d]); 
    when -12:  s:=[ -15, 22 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -27:  s:=[ -480, 4048 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -4:   s:=[ 1, 0 ]; E:=EllipticCurve([d,0]); 
    when -16:  s:=[ -11, 14 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -7:   s:=[ -2835, -71442 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -28:  s:=[ -595, 5586 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -8:   s:=[ -4320, 96768 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -11:  s:=[ -9504, 365904 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -19:  s:=[ -608, 5776 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -43:  s:=[ -13760, 621264 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -67:  s:=[ -117920, 15585808 ]; E:=QuadraticTwist(EllipticCurve(s),d);
    when -163: s:=[ -34790720, 78984748304 ]; E:=QuadraticTwist(EllipticCurve(s),d);
  end case;
  as:=aInvariants(E);
  A:=as[4];
  B:=as[5];
  return E,[A,B];
end function;


function PrimitiveDivisionPolynomial(E,m)
// m-th division polynomial of the elliptic curve E divided by the d-th division polynomial for any proper divisor d of m
// By Drew
  local f;
  f:=DivisionPolynomial(E,m);
  for d in Divisors(m) do 
  	if d gt 1 and d lt m then 
  		f := ExactQuotient(f,$$(E,d)); 
  	end if; 
  end for;
  return f;
end function;


////////////////////////////////////////////////////////////////////////////////////////////////

function CMInvariants(E)
   boo, D:=HasComplexMultiplication(E); 
   assert boo;                                                                                                                                                                         
   W:=WeierstrassModel(E);
   A:=aInvariants(W)[4];
   B:=aInvariants(W)[5];
   case D:
    when -3:   s:=[ 0, 1 ]; d:=B; d:=PowerFreePart(d,6);
    when -12:  s:=[ -15, 22 ]; d:=(B/s[2])/(A/s[1]); d:=PowerFreePart(d,2);
    when -27:  s:=[ -480, 4048 ]; d:=(B/s[2])/(A/s[1]);d:=PowerFreePart(d,2);
    when -4:   s:=[ 1, 0 ]; d:=A;d:=PowerFreePart(d,4);
    when -16:  s:=[ -11, 14 ]; d:=(B/s[2])/(A/s[1]);d:=PowerFreePart(d,2);
    when -7:   s:=[ -2835, -71442 ]; d:=(B/s[2])/(A/s[1]);d:=PowerFreePart(d,2);
    when -28:  s:=[ -595, 5586 ]; d:=(B/s[2])/(A/s[1]);d:=PowerFreePart(d,2);
    when -8:   s:=[ -4320, 96768 ]; d:=(B/s[2])/(A/s[1]);d:=PowerFreePart(d,2);
    when -11:  s:=[ -9504, 365904 ]; d:=(B/s[2])/(A/s[1]);d:=PowerFreePart(d,2);
    when -19:  s:=[ -608, 5776 ]; d:=(B/s[2])/(A/s[1]);d:=PowerFreePart(d,2);
    when -43:  s:=[ -13760, 621264 ]; d:=(B/s[2])/(A/s[1]);d:=PowerFreePart(d,2);
    when -67:  s:=[ -117920, 15585808 ]; d:=(B/s[2])/(A/s[1]);d:=PowerFreePart(d,2);
    when -163: s:=[ -34790720, 78984748304 ]; d:=(B/s[2])/(A/s[1]);d:=PowerFreePart(d,2);
   end case;
   EE:=CM_EllipticCurve(D,d);
   assert IsIsomorphic(E,EE);
   printf "CM Invariants:\n cm= %o   -   k=%o)\n", Abs(D),d;
   return Abs(D),d;
end function;