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

function TunnelsCriterion(n)
   n := SquareFree(n);
   bnd := Floor(Sqrt(n))+1;
   box := [-bnd..bnd]; 
   set8 := [];
   set32 := [];
   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];
      for x in box do
         for z in box do
            t, y := IsSquare(n - 2*x^2 - 8*z^2);
            if t and y in box then
               Append(~set8, <x,y,z>);
            end if;
            t, y := IsSquare(n - 2*x^2 - 32*z^2);
            if t and y in box then
               Append(~set32, <x,y,z>);
            end if;
         end for;
      end for;
   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];
      for x in box do
         for z in box do
            t, y := IsSquare((n div 2) - 4*x^2 - 8*z^2);
            if t and y in box then
               Append(~set8, <x,y,z>);
            end if;
            t, y := IsSquare((n div 2) - 4*x^2 - 32*z^2);
            if t and y in box then
               Append(~set32, <x,y,z>);
            end if;
         end for;
      end for;
   end if;
   return #set8 eq 2*#set32, set8, 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;
   R := RealField(2000);
   a := R!Numerator(x);
   assert IsSquare(a);
   b := R!Denominator(x);
   assert IsSquare(b);
   sqrt := Round(Sqrt(a))/Round(Sqrt(b));
   assert sqrt^2 eq x;
   return sqrt;
end function;


function CongruentTriangle(P)   
   assert Type(P) eq PtEll;
   // P should be a point on y^2 = x^3 - n^2*x. 
   n := Integers()!RationalSqrt(-aInvariants(Curve(Parent(P)))[4]);
   P := 2*P;
   x := P[1]/P[3];
   y := P[2]/P[3];   
   u := RationalSqrt(x);
   v := y/u;
   t := Denominator(u);

   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;


function CongruentRep(n)
   // Assume that n is a congruent number.  Try to 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 didn't find a representation.";
      else
	 print n, "is not a congruent number.";
      end if;
      return 0;
   end if;
   return CongruentTriangle(f(G.3)), f(G.3);
end function;