/* y^2 = x*(x-1)*(x-a) = x^3 - (a+1)x^2 + ax */

pi := Pi(RealField());

function Ell(a)
   return EllipticCurve([0, -a-1, 0, a,0]);
end function;

function anlist(E,prec)
   return Eltseq(qEigenform(E,prec));
end function;

function LValue(E, anlist, A, eps)
   N := Conductor(E);
   rootN := Sqrt(N);
   return &+[anlist[n]/n*(Exp(-2*pi*n*A/rootN) + eps*Exp(-2*pi*n/(A*rootN))) :
                n in [1..#anlist] | GCD(N,n) eq 1];
end function;


function L(E, anlist, stop)
   N := Conductor(E);
   rootN := Sqrt(N);
   //   return &+[anlist[n]/n*Exp(-2*pi*n/rootN) : n in [1..stop] | GCD(N,n) eq 1];
   return &+[anlist[n]/n*Exp(-2*pi*n/rootN) : n in [1..stop]];
end function;


function CongruentLValue(n,prec)
   E := EllipticCurve([-n^2,0]);
   an := anlist(E,prec);
   if IsOdd(n) then
      if KroneckerSymbol(-2,n) eq -1 then
	 return 0;
      end if;
   else
      if KroneckerSymbol(-2,n div 2) eq -1 then
	 return 0;
      end if;
   end if;
   return LValue(E,an,1,+1);
end function;

function PSeq(n1,n2)
   return [p : p in [n1..n2] |IsPrime(p)];
end function;

procedure Convergents(E,an,stop)
   for s in [q : q in [2..Min(#an,stop)]|IsPrime(q)] do
      s, ":  ", an[s], "  \t", L(E,an,s);
   end for;
end procedure;


/* I copied this from the mathscinet review of Ono's "Partitions into distinct
parts and elliptic curves */


function Waldspurger(prec)
   R<q> := PowerSeriesRing(Rationals());
   return (1 + &+[2*q^(2*n^2) : n in [1..Ceiling(Sqrt(prec+1))]]) *
      q * &*[(1-q^(8*n))*(1-q^(16*n)) : n in [1..prec div 8+10]]
      +O(q^(prec+1));
end function;

function CongruentNumberCurve(n)
   return EllipticCurve([-n^2,0]);
end function;


function IsSquareFree(n)
   return n le 2 select true else &*[F[2] : F in Factorization(n)] eq 1;
end function;

procedure DisplayCongruentLValues(start,prec)
   n := start;
   while true do
      if IsSquareFree(n) then
	 L := CongruentLValue(n,prec);
	 E := CongruentNumberCurve(n);
	 print n, ": \t", L, " \t", Abs(Periods(E)[1]);
      end if;
      n +:= 1;
   end while;
end procedure;


function TunnelsCriterion(n)
   n := SquareFree(n);
   bnd := Floor(Sqrt(n))+1;
   box := [-bnd..bnd]; 
   if IsOdd(n) then
      set8 :=  [<x,y,z> : x in box, y in box, z in box | 2*x^2 + y^2 + 8*z^2 eq n];
      set32 := [<x,y,z> : x in box, y in box, z in box | 2*x^2 + y^2 + 32*z^2 eq n];
   else
      set8 :=  [<x,y,z> : x in box, y in box, z in box | 4*x^2 + y^2 + 8*z^2 eq n div 2];
      set32 := [<x,y,z> : x in box, y in box, z in box | 4*x^2 + y^2 + 32*z^2 eq n div 2];
   end if;
//   "First set := ", set8;
//   "Second set := ", set32;

   return #set8 eq 2*#set32;
end function;


function abRepresentation(X,Y,Z)
   /*************************************************************
   Assume that (X, Y, Z) form a primitive Pythagorean triple
   with Y even.  This function returns relatively prime integers
   a, b with X = a^2 - b^2, Y = 2ab, and Z = a^2 + b^2.
   *************************************************************/
   assert IsEven(Y);
   assert X^2 + Y^2 eq Z^2;
   assert GCD([X, Y, Z]) eq 1;
   u := X/Z;
   v := Y/Z;
   alpha := v/(u+1);
   a := Denominator(alpha);
   b := Numerator(alpha);
   assert X eq a^2 - b^2;
   assert Y eq 2*a*b;
   assert Z eq a^2 + b^2;
   return a,b;
end function;


function RationalSqrt(x)
   assert Type(x) in {RngIntElt, FldRatElt};
   assert x ge 0;
   a := Numerator(x);
   assert IsSquare(a);
   b := Denominator(x);
   assert IsSquare(b);
   return Round(Sqrt(a))/Round(Sqrt(b));
end function;



function CongruentTriangle(n)
   // assuming n is a congruent number, find a triangle that
   // realizes it as a congruent number.

   E := CongruentNumberCurve(n);
   G,f := MordellWeilGroup(E);
   if Ngens(G) le 2 then
      l,r := RankBounds(E);
      if r gt 0 then
	 print n, "is probably a congruent number, but I can't find the representation.";
      else
	 print n, "is not a congruent number.";
      end if;
      return 0;
   end if;
   print "A solution is ", f(G.3);
   P := f(2*G.3);
   x := P[1]/P[3];
   y := P[2]/P[3];   
   u := RationalSqrt(x);
   v := y/u;
   t := Denominator(u);

   // print t^2*v, t^2*n, t^2*x;
   a,b := abRepresentation(t^2*v, t^2*n, t^2*x);
   return Abs(2*a/t), Abs(2*b/t), Abs(2*u);
end function;