import "../../../../modsym/core.m": LiftToCosetRep;
import "../../../../modsym/analytic.m": PeriodIntegral;
import "../../../../modsym/period.m" : RationalPeriodMapping;

forward LatticeCoordinates, LinearFractionalTransformation, Integral;

intrinsic TwistsOfBigRank(E::CrvEll, start::RngIntElt, stop::RngIntElt) -> SeqEnum
{}
   S := ModularSymbols(E);
   M := AmbientSpace(S);
   phi := RationalPeriodMapping(S);

   ans := [];

   for p in [q : q in [start..stop] |IsPrime(q) and Conductor(E) mod q ne 0] do
      if p mod 4 eq 1 then
	 G<eps> := DirichletGroup(p);
	 d := p;
      else
	 G<a,b> := DirichletGroup(4*p);
	 eps := a*b;
	 d := 4*p;
      end if;
      if phi(VectorSpace(M)!Eltseq(TwistedWindingElementModular(M,1,eps))) eq 0 then
	 Append(~ans,<Evaluate(eps,-11),-d>);
	 ans;
      end if;
   end for;
   
   return ans;
end intrinsic;
   
intrinsic ModularPoint(E::CrvEll, f::RngSerElt, D::RngIntElt) -> FldPrElt, FldPrElt
   {Computes the image in C/Lambda = E(C) of the trace of the CM point corresponding
   to the order of discriminant D.  The discriminant must be of the form D_K*p^2,
   where p is the conductor of E, and p is prime.  f must be the newform corresponding to E.}
   
   require D lt 0 : "Argument 2 must be less than 0.";

   require E eq MinimalModel(E) : "Argument 1 must be minimal.";
   p := Conductor(E);
   require IsPrime(p) : "Argument 1 must have prime conductor.";
   
//   f := qEigenform(E,prec);

   beta := Sqrt(D);
   Pic := ClassGroupRepresentatives(D,beta);
   Pic;

   z :=  &+ [ Integral(f, alpha) : alpha in Pic];

   return LatticeCoordinates(E,z);
end intrinsic;

intrinsic CMPoint(E::CrvEll, alpha::FldPrElt, prec::RngIntElt) -> FldPrElt, FldPrElt
   {}
   require E eq MinimalModel(E) : "Argument 1 must be minimal.";
   
   f := qEigenform(E,prec);
   N := Conductor(E);
   P1:= P1Classes(N);
   R := [LiftToCosetRep(x, N) : x in P1];
   r_of_alpha := [LinearFractionalTransformation(r,alpha) : r in R];
   //   return r_of_alpha;
   z := &+ [ Integral(f, ralp) : ralp in r_of_alpha];
   return LatticeCoordinates(E,z);
end intrinsic;

intrinsic CMModSym(M::ModSym, alpha::FldRatElt) -> ModSymElt
{}
   N := Level(M);
   P1:= P1Classes(N);
   R := [LiftToCosetRep(x, N) : x in P1];
   r_of_alpha := [Cusps()|LinearFractionalTransformation(r,alpha) : r in R];
   ans := M!0;
   A := AmbientSpace(M);
   for ralp in r_of_alpha do
      ans := ans + A!<PolynomialRing(Rationals(),2)!1, [Cusps()|Infinity(), ralp]>;
   end for;
   return ans;
end intrinsic;

intrinsic ClassGroupRepresentatives(D::RngIntElt) -> SeqEnum
   {}
   K<alp> := QuadraticField(D);
   r := Round(Sqrt(Integers()!(D/alp^2)));
   assert r^2 eq D/alp^2;
   beta := r * alp;   // so beta^2 = D.
   assert beta^2 eq D;
   /*   beta^2 = r^2 * alp^2 = D/alp^2 * alp^2 = D.    */
   return ClassGroupRepresentatives(D, beta);
end intrinsic;

intrinsic ClassGroupRepresentatives(D::RngIntElt, beta::FldElt) -> SeqEnum
   {Reduced representatives of the class group of
   discriminant D, as on page 21 of Number Theory II.
   Beta should be the square root of D, in some field.}

   require D lt 0 : "Argument 1 must be less than 0.";
   
   bound := Round(Sqrt(-D/3));
   ans := [];

   for a in [1..bound] do
      for b in [-a..a-1] do
	 c := (b^2 - D)/(4*a);
	 if Denominator(c) eq 1 then
	    c := Integers()!c;
	    if GCD([a,b,c]) eq 1 then
	       if b le 0 then
		  if c ge a then
		     Append(~ans,[a,b]);
		  end if;
	       else
		  if c gt a then
		     Append(~ans,[a,b]);
		  end if;
	       end if;
	    end if;
	 end if;
      end for;
   end for;

   return [(-z[2] + beta)/(2*z[1]) : z in ans];
end intrinsic;
      
function Integral(f, alpha)
   return PeriodIntegral(f, PolynomialRing(Rationals())!1, alpha);
end function;

function LinearFractionalTransformation(r, alpha)
   if r[3]*alpha + r[4] eq 0 then
      return Infinity();
   end if;
   return (r[1]*alpha + r[2])/(r[3]*alpha + r[4]);
end function;

function LatticeCoordinates(E,z)
   assert Type(E) eq CrvEll;
   om1, om2 := Explode(Periods(E));
   z := z/om1;
   t := om2/om1;
   x := Real(z);
   y := Imaginary(z);
   return y/Imaginary(t), x - Real(t)*(y/Imaginary(t));
end function;

intrinsic Go(E::CrvEll, prec::RngIntElt)
{}
   P := [p : p in [3..97] | IsPrime(p)];
   for p in P do
      if p mod 4 eq 1 then 
         alpha := Sqrt(-p);
      else
         alpha := (1+Sqrt(-p))/2;
      end if;
      CMPoint(E, alpha, prec);
   end for;
end intrinsic;

procedure TestCGR()
   D := -Random(1,2000)*4  + Random(0,1);
   print D, ClassNumber(D);
   time h := #ClassGroupRepresentatives(D);
   assert h eq ClassNumber(D);
end procedure;


