/************************************************************
 * Implementation of X. Wang's algorithms from              *
 *    "2-dimensional simple factors of J_0(N)."             *
 *                                                          *
 * William A. Stein, December 1999                          *
 * with the assistance of Adam Logan.                       *
 ************************************************************/

/*
function I2(a0,a1,a2,a3,a4,a5,a6)
  return -240*a0*a6+40*a1*a5-16*a2*a4+6*a3^2;
end function;

function I4(a0,a1,a2,a3,a4,a5,a6)
  return 48*a0*a4^3+48*a2^3*a6+4*a2^2*a4^2+1620*a0^2*a6^2+36*a1*a3^2*a5-12*a1*a3*a4^2-12*a2^2*a3*a5+300*a1^2*a4*a6+300*a0*a5^2*a2+324*a0*a6*a3^2-504*a0*a2*a4*a6-180*a0*a4*a3*a5-180*a1*a3*a2*a6+4*a1*a4*a2*a5-540*a0*a5*a1*a6-80*a1^2*a5^2;
end function;

function I6(a0,a1,a2,a3,a4,a5,a6) 
  return 176*a1^2*a5^2*a3^2+64*a1^2*a5^2*a4*a2+1600*a1^3*a5*a4*a6+1600*a1*a5^3*a0*a2-2240*a1^2*a5^2*a0*a6-160*a0*a4^4*a2-96*a0^2*a4^3*a6+60*a0*a4^3*a3^2+72*a1*a3^4*a5-24*a1*a3^3*a4^2+2250*a1^3*a3*a6^2-160*a2^4*a4*a6-96*a2^3*a0*a6^2+60*a2^3*a3^2*a6-24*a2^2*a3^3*a5+8*a2^2*a3^2*a4^2-900*a2^2*a1^2*a6^2-24*a2^3*a4^3-36*a2^4*a5^2-36*a1^2*a4^4+424*a0*a4^2*a2^2*a6+492*a0*a4^2*a2*a3*a5+20664*a0^2*a4*a6^2*a2+3060*a0^2*a4*a6*a3*a5-468*a0*a4*a3^2*a2*a6-198*a0*a4*a3^3*a5-640*a0*a4*a2^2*a5^2+3472*a0*a4*a2*a5*a1*a6-18600*a0*a4*a1^2*a6^2-876*a0*a4^2*a1*a6*a3+492*a1*a3*a2^2*a4*a6-238*a1*a3^2*a2*a4*a5+76*a1*a3*a2*a4^3+3060*a1*a3*a0*a6^2*a2+1818*a1*a3^2*a0*a6*a5-198*a1*a3^3*a2*a6+26*a1*a3*a2^2*a5^2-1860*a1^2*a3*a2*a5*a6+330*a1^2*a3^2*a6*a4+76*a2^3*a4*a3*a5-876*a2^2*a0*a6*a3*a5+616*a3^2*a5*a1*a6+2250*a0^2*a5^3*a3-900*a0^2*a5^2*a4^2-10044*a0^2*a6^2*a3^2+162*a0*a6*a3^4+28*a1*a4^2*a2^2*a5-640*a1^2*a4^2*a2*a6+26*a1^2*a4^2*a3*a5-1860*a1*a4*a0*a5^2*a3+616*a1*a4^3*a0*a5-18600*a0^2*a5^2*a6*a2+59940*a0^2*a5*a6^2*a1+330*a0*a5^2*a3^2*a2-119880*a0^3*a6^3-320*a1^3*a5^3;
end function;


function I10(a0,a1,a2,a3,a4,a5,a6) 
   R<x> := 
   disc. of a6*x^6+a5*x^5+a4*x^4+a3*x^3+a2*x^2+a1*x+a0
*/




/*
i1 = I2^5/I10
i2 = I2^3*I4/I10
i3 = I2^2*I6/I10


i1 = h12^5/h10^6
i2 = h12^3*h4/h10^4
i3 = h12^2*h16/h10^4
*/


/*********************************************************/
/*                                                       */
/* Compute the Theta null-values (see page 190 of Wang). */
/*                                                       */
/*********************************************************/

intrinsic ThetaNull(a::ModMatFldElt, b::ModMatFldElt,
                    Omega::AlgMatElt,
                    prec::RngIntElt) -> FldPrElt 
{ Compute the Theta null-values (see page 190 of Wang) summing
  in a box of radius prec about the origin. 
  Here a and b are in RMatrixSpace(ComplexField(),2,1), and
  Omega is in MatrixAlgebra(ComplexField(),2). }

   // Later re-program so prec is the number of decimal digits of precision?

   return ThetaNull(a, b, Omega, 
                    0.0, -1, prec);
end intrinsic;

intrinsic ThetaNull(a::ModMatFldElt, b::ModMatFldElt, 
                    Omega::AlgMatElt,
                    x::FldPrElt,  // sum so far
                    r1::RngIntElt, r2::RngIntElt) -> FldPrElt 
{  r1 is the radius summed to so far. }

   V  := RMatrixSpace(ComplexField(),2,1);
   Z2 := [ V![c,d] : c in [-r2..r2], d in [-r2..r2] |
                     Abs(c) gt r1 or Abs(d) gt r1 ] ; 

   C<i> := ComplexField();
   pi   := Pi(C);

   return x + 
      &+[Exp(pi*i*(Transpose(n+a)*Omega*(n+a))[1,1]+
             2*pi*i*(Transpose(n+a)*b)[1,1]) : 
         n in Z2];
end intrinsic;

intrinsic ThetaNullEvens(Omega::AlgMatElt, prec::RngIntElt) -> SeqEnum
{Return a sequence [v0, v1, v3, v5, v21, v23, v25, v41, v43, v45]
 consisting of the 10 even theta-null values.}
   V  := RMatrixSpace(ComplexField(),2,1);
   return [ ThetaNull(V![0,0],V![0,0], Omega, prec),
            ThetaNull(V![0.5,0],V![0,0], Omega, prec),
            ThetaNull(V![0,0.5],V![0.5,0], Omega, prec),
            ThetaNull(V![0,0],V![0.5,0.5], Omega, prec),
            ThetaNull(V![0,0],V![0.5,0], Omega, prec),
            ThetaNull(V![0.5,0.5],V![0,0], Omega, prec),
            ThetaNull(V![0.5,0],V![0,0.5], Omega, prec),
            ThetaNull(V![0.5,0.5],V![0.5,0.5], Omega, prec),
            ThetaNull(V![0,0],V![0,0.5], Omega, prec),
            ThetaNull(V![0,0.5],V![0,0], Omega, prec) 
          ];
end intrinsic;

intrinsic SiegelForms(v::SeqEnum) -> SeqEnum
{Given the sequence of v's computed by ThetaNullEvans, returns
the values [h4, h10, h12, h16] of certain Siegel modular forms.}
   v0, v1, v3, v5, v21, v23, v25, v41, v43, v45 
       := Explode(v);

   h4 := v0^8+v1^8+v21^8+v23^8+v25^8+v3^8+v41^8+v43^8+v45^8+v5^8;
   h10 := (v0*v1*v21*v23*v25*v3*v41*v43*v45*v5)^2;
   h12 := (v0*v1*v21*v23*v25*v41)^4 + (v0*v21*v23*v25*v3*v43)^4 + (v1*v23*v25*v3*v41*v45)^4 + (v1*v21*v25*v3*v43*v45)^4 +
          (v0*v1*v21*v41*v43*v45)^4 + (v0*v23*v3*v41*v43*v45)^4 + (v0*v1*v21*v3*v41*v5)^4 + (v0*v1*v23*v3*v43*v5)^4 + 
          (v1*v23*v25*v41*v43*v5)^4 + (v21*v25*v3*v41*v43*v5)^4 + (v0*v21*v23*v25*v45*v5)^4 + (v0*v1*v25*v3*v45*v5)^4 + 
          (v21*v23*v3*v41*v45*v5)^4 + (v1*v21*v23*v43*v45*v5)^4 + (v0*v25*v41*v43*v45*v5)^4;
   h16 := 
v0^4*v1^4*v21^4*v23^4*v25^4*v3^8*v41^4 + v0^4*v1^8*v21^4*v23^4*v25^4*v3^4*v43^4 + v0^4*v21^4*v23^4*v25^4*v3^4*v41^8*v43^4 + v0^4*v1^4*v21^4*v23^4*v25^4*v41^4*v43^8 + 
v0^8*v1^4*v23^4*v25^4*v3^4*v41^4*v45^4 + v1^4*v21^8*v23^4*v25^4*v3^4*v41^4*v45^4 + v0^8*v1^4*v21^4*v25^4*v3^4*v43^4*v45^4 + v1^4*v21^4*v23^8*v25^4*v3^4*v43^4*v45^4 + 
v0^4*v1^4*v21^4*v23^8*v41^4*v43^4*v45^4 + v0^4*v1^4*v21^4*v25^8*v41^4*v43^4*v45^4 + v0^4*v1^8*v23^4*v3^4*v41^4*v43^4*v45^4 + v0^4*v21^8*v23^4*v3^4*v41^4*v43^4*v45^4 + 
v0^4*v23^4*v25^8*v3^4*v41^4*v43^4*v45^4 + v0^4*v1^4*v21^4*v3^8*v41^4*v43^4*v45^4 + v1^4*v21^4*v25^4*v3^4*v41^8*v43^4*v45^4 + v1^4*v23^4*v25^4*v3^4*v41^4*v43^8*v45^4 + 


v0^4*v1^4*v21^4*v23^4*v25^4*v41^4*v45^8 + v0^4*v21^4*v23^4*v25^4*v3^4*v43^4*v45^8 + v0^4*v1^4*v21^4*v23^8*v3^4*v41^4*v5^4 + v0^4*v1^4*v21^4*v25^8*v3^4*v41^4*v5^4 +
 v0^4*v1^4*v21^8*v23^4*v3^4*v43^4*v5^4 + v0^4*v1^4*v23^4*v25^8*v3^4*v43^4*v5^4 + v0^8*v1^4*v23^4*v25^4*v41^4*v43^4*v5^4 + v1^4*v21^8*v23^4*v25^4*v41^4*v43^4*v5^4 +
v0^8*v21^4*v25^4*v3^4*v41^4*v43^4*v5^4 + v1^8*v21^4*v25^4*v3^4*v41^4*v43^4*v5^4 + v21^4*v23^8*v25^4*v3^4*v41^4*v43^4*v5^4 + v1^4*v23^4*v25^4*v3^8*v41^4*v43^4*v5^4 +
v0^4*v1^4*v23^4*v3^4*v41^8*v43^4*v5^4 + v0^4*v1^4*v21^4*v3^4*v41^4*v43^8*v5^4 + v0^4*v1^8*v21^4*v23^4*v25^4*v45^4*v5^4 + v0^4*v1^4*v21^8*v25^4*v3^4*v45^4*v5^4 +


v0^4*v1^4*v23^8*v25^4*v3^4*v45^4*v5^4 + v0^4*v21^4*v23^4*v25^4*v3^8*v45^4*v5^4 + v0^8*v21^4*v23^4*v3^4*v41^4*v45^4*v5^4 + v1^8*v21^4*v23^4*v3^4*v41^4*v45^4*v5^4 +
v21^4*v23^4*v25^8*v3^4*v41^4*v45^4*v5^4 + v0^4*v21^4*v23^4*v25^4*v41^8*v45^4*v5^4 + v0^4*v1^4*v25^4*v3^4*v41^8*v45^4*v5^4 + v0^8*v1^4*v21^4*v23^4*v43^4*v45^4*v5^4 +
v1^4*v21^4*v23^4*v25^8*v43^4*v45^4*v5^4 + v1^4*v21^4*v23^4*v3^8*v43^4*v45^4*v5^4 + v0^4*v1^8*v25^4*v41^4*v43^4*v45^4*v5^4 + v0^4*v21^8*v25^4*v41^4*v43^4*v45^4*v5^4 +
v0^4*v23^8*v25^4*v41^4*v43^4*v45^4*v5^4 + v0^4*v25^4*v3^8*v41^4*v43^4*v45^4*v5^4 + v1^4*v21^4*v23^4*v41^8*v43^4*v45^4*v5^4 + v0^4*v21^4*v23^4*v25^4*v43^8*v45^4*v5^4 +
v0^4*v1^4*v25^4*v3^4*v43^8*v45^4*v5^4 + v21^4*v23^4*v3^4*v41^4*v43^8*v45^4*v5^4 + v0^4*v1^4*v21^4*v3^4*v41^4*v45^8*v5^4 + v0^4*v1^4*v23^4*v3^4*v43^4*v45^8*v5^4 +
 v1^4*v23^4*v25^4*v41^4*v43^4*v45^8*v5^4 + v21^4*v25^4*v3^4*v41^4*v43^4*v45^8*v5^4 + v0^4*v1^4*v21^4*v23^4*v25^4*v41^4*v5^8 + v0^4*v21^4*v23^4*v25^4*v3^4*v43^4*v5^8 +
 v1^4*v23^4*v25^4*v3^4*v41^4*v45^4*v5^8 + v1^4*v21^4*v25^4*v3^4*v43^4*v45^4*v5^8 + v0^4*v1^4*v21^4*v41^4*v43^4*v45^4*v5^8 + v0^4*v23^4*v3^4*v41^4*v43^4*v45^4*v5^8;

   return [h4, h10, h12, h16];
end intrinsic;

intrinsic AbsoluteInvariants(h::SeqEnum) -> SeqEnum
{Given the sequence of values [h4, h10, h12, h16], this function
returns the three absolute invariants i1, i2, and i3 of the principally
polarized abelian variety.}
   h4, h10, h12, h16 := Explode(h);
   i1 := h12^5/h10^6;
   i2 := h12^3*h4/h10^4;
   i3 := h12^2*h16/h10^4;
   return [i1, i2, i3];
end intrinsic;

intrinsic AbsoluteInvariants(Omega::AlgMatElt,
                             prec::RngIntElt) -> SeqEnum
{This function returns the three absolute invariants i1, i2, and i3 
 of the principally polarized abelian variety defined by the
 period matrix Omega.}
   v := ThetaNullEvens(Omega, prec);
   h := SiegelForms(v);
   return AbsoluteInvariants(h);
end intrinsic;



/*

a:=-0.726444449+i*0.897558797;  
b:=-0.363222225+i*0.448779399;
c:=-0.363222225+i*0.448779399;
d:=0.2735555509+i*0.897558797;  
*/
