/****-*-magma-* EXPORT DATE: 2004-03-08 ************************************
                                                                            
                     MODABVAR: Modular Abelian Varieties in MAGMA
 
                              William A. Stein        
                         
   FILE: modabvar.m
   DESC: Basic ModAbVar functionality.
      
 ***************************************************************************/


/*
HANDBOOK_TITLE: Creation and Basic Functions

BEGIN_HANDBOOK_INTRO

The functions described below are for creating modular abelian abelian
varieties, combining them together in various ways, and obtaining
simple information about them. 

Modular abelian varieties are much less restricted than spaces of
modular symbols, in that one can take arbitrary finite direct sum.


END_HANDBOOK_INTRO



BEGIN_HANDBOOK_EXAMPLES
\beginex{modabvar1}%\>--------------
The simplest abelian variety is an abelian variety of dimension $0$,
i.e., a point.
\begincode
> A := ZeroModularAbelianVariety(); A;
Modular abelian variety ZERO of dimension 0 and level 1 over Q
\endcode
We create the Jacobian $J_0(22)$ of the modular curve
$X_0(22)$ as follows:
\begincode
> J := JZero(22); J;
Modular abelian variety JZero(22) of dimension 2 and level 2*11 over Q
\endcode
Notice that $A$ is a subset of $J_0(22)$.
\begincode
> A subset J;
true
\endcode
We can also create the higher weight analogues $J_0(N,k)$ of $J_0(N)$,
which are motives defined over~$\Q$.  Many computations that make
sense for $J_0(N)$ also make sense for these higher weight analogues.
\begincode
> J4 := JZero(22,4); J4;
Modular motive JZero(22,4) of dimension 7 and level 2*11 over Q
> IsOnlyMotivic(J4);
true
\endcode

One can also create $J_1(N)$:
\begincode
> JOne(22); 
Modular abelian variety JOne(22) of dimension 6 and level 2*11 over Q
\endcode

For efficiency purposes, it is often much quicker to do computations
working only with the $+1$ quotient of $H_1(J_0(N),\Z)$ by $*=1$, 
or using only the $-1$ quotient.  These computations will be off
by powers of $2$, or subgroups will be halved and off by a power of~$2$.
Nonetheless, if you know what you are doing, such computations can
be very useful.  Create $J_0(N)$, but working only with the $+1$
quotient of homology as follows:
\begincode
> Jplus := JZero(22,2,+1); Jplus;
Modular abelian variety JZero(22) of dimension 2 and level 2*11 over Q 
with sign 1
> Sign(Jplus);
1
> Sign(JZero(22));
0
\endcode
Notice that the sign is printed out when it is $1$ or $-1$.  Also,
sign $0$ is shorthand for ``no sign'', i.e., working with the full
homology.  It is not possible to take the direct sum of abelian
varieties with different signs.

Let $\eps:(\Z/N\Z)^*\to \Q(\zeta_n)^*$ be Dirichlet character. The command 
{\tt ModularAbelianVariety(eps)} creates the modular abelian
variety {\em over~$\Q$} corresponding to the cusp forms with
Dirichlet character any Galois conjugate of~$\eps$.
\begincode
> G<eps> := DirichletGroup(22,CyclotomicField(EulerPhi(22)));
> Order(eps);
10
> Conductor(eps);
11
> A := ModularAbelianVariety(eps); A;
Modular abelian variety of dimension 0 and level 2*11 over Q
> A := ModularAbelianVariety(eps^2); A;
Modular abelian variety of dimension 4 and level 2*11 over Q
\endcode

Let~$H$ be a sugroup of $G=(\Z/N\Z)^*$.
The group $(\Z/N\Z)^*$ acts by diamond bracket operators as
a group of automorphisms on $X_1(N)$, and 
the modular curve $X_H(N)$ is the quotient of $X_1(N)$
by the action of~$H$.  Thus if $H=1$, then $X_H(N)=X_1(N)$,
and if $H=G$, then $X_H(N)=X_0(N)$.
The command {\tt JH(N,d)} creates an abelian
variety {\em isogenous to} the Jacobian of $X_H(N)$, 
where~$d$ is the index of $H$ in $G$ (thus $d=1$ 
corresponds to $J_0(N)$). A weight~$k$ and sign (either 
$0$ or $\pm 1$) can be provided as optional parameters.
More precisely, {\tt JH(N,d)} is the product of 
$J(\eps)$, where~$\eps$ varies over Dirichlet characters
such that $\eps(H)=\{1\}$.
\begincode
> JH(22,1);
Modular abelian variety JZero(22) of dimension 2 and level 2*11 over Q
> JH(22,10);
Modular abelian variety Js(22) of dimension 6 and level 2*11 over Q
> JH(22,2);
Modular abelian variety J_H(22) of dimension 2 and level 2*11 over Q
\endcode
As a shortcut, the command {\tt Js(N)} creates a 
modular abelian variety that is isogenous to $J_1(N)$.  
Thus {\tt Js(N)} is the same as {JH(N,d)}, where $d$ is the
order of $(\Z/N\Z)^*$.
\begincode
> Js(22);
Modular abelian variety Js(22) of dimension 6 and level 2*11 over Q
> JH(22,10) eq Js(22);
true
\endcode

We can also create modular abelian varieties attached to 
spaces of modular forms and spaces of modular symbols:
\begincode
> ModularAbelianVariety(ModularForms(22));
Modular abelian variety of dimension 2 and level 2*11 over Q
> ModularAbelianVariety(ModularSymbols(22));
Modular abelian variety of dimension 2 and level 2*11 over Q
\endcode
Here is another example, in which the space of modular forms
is for $\Gamma_1(N)$.
\begincode
> M := ModularForms(Gamma1(25)); M;
Space of modular forms on Gamma_1(25) of weight 2 and dimension 
39 over Integer Ring.
> S := CuspidalSubspace(M); S;
Space of modular forms on Gamma_1(25) of weight 2 and dimension 
12 over Integer Ring.
> A := ModularAbelianVariety(S); A;
Modular abelian variety of dimension 12 and level 5^2 over Q
\endcode

We can also construct abelian varieties attached to newforms.
\begincode
> S := CuspForms(43);
> N := Newforms(S); N;
[* [*
q - 2*q^2 - 2*q^3 + 2*q^4 - 4*q^5 + 4*q^6 + O(q^8)
*], [*
q + a*q^2 - a*q^3 + (-a + 2)*q^5 - 2*q^6 + (a - 2)*q^7 + O(q^8),
q + b*q^2 - b*q^3 + (-b + 2)*q^5 - 2*q^6 + (b - 2)*q^7 + O(q^8)
*] *]
> f := N[2][1]; f;
q + a*q^2 - a*q^3 + (-a + 2)*q^5 - 2*q^6 + (a - 2)*q^7 + O(q^8)
> A := ModularAbelianVariety(f); A;
Modular abelian variety Af of dimension 2 and level 43 over Q
> Newform(A);
q + a*q^2 - a*q^3 + (-a + 2)*q^5 - 2*q^6 + (a - 2)*q^7 + O(q^8)
\endcode
When possible, we can also obtain a newform that gives rise to
an abelian variety.
\begincode
> J := JZero(43);
> D := Decomposition(J); D;
[
    Modular abelian variety 43A of dimension 1, level 43 and 
    conductor 43 over Q,
    Modular abelian variety 43B of dimension 2, level 43 and 
    conductor 43^2 over Q
]
> Newform(D[2]);
q + a*q^2 - a*q^3 + (-a + 2)*q^5 - 2*q^6 + (a - 2)*q^7 + O(q^8)
\endcode

Modular abelian varieties may be described by a label, which is a
string like "43B".  Continuing the previous code, we have:
\begincode
> A := ModularAbelianVariety("43B");
> A eq D[2];
true
\endcode
Notice that we created $A$ above using the command
{\tt ModularAbelianVariety} with a string as an argument.
This is only guaranteed to work on labels of the form
a level followed by an isogeny code. 

The integral, rational and real homology of a modular abelian variety~$A$
is $H_1(A,R)$, where $R=\Z,\R,\Q$, respectively.  
Thus this homology is just a free module over~$R$ 
of dimension twice $\dim(A)$ (unless the sign of~$A$ is $\pm 1$, 
in which case the homology is of dimension $\dim(A)$). 
\begincode
> J := JZero(22); J;
Modular abelian variety JZero(22) of dimension 2 and level 2*11 over Q
> IntegralHomology(J);
Standard Lattice of rank 4 and degree 4
> RationalHomology(J);
Full Vector space of degree 4 over Rational Field
> RealHomology(J);
Full Vector space of degree 4 over Real Field
> J := JZero(22, 2, +1); J;
Modular abelian variety JZero(22) of dimension 2 and level 2*11 over Q 
with sign 1
> RationalHomology(J);
Full Vector space of degree 2 over Rational Field
\endcode

Using the product operator we can take the direct product
of any two modular abelian varieties, even ones of different
weights or levels.  We first illustrate taking the product
of two abelian subvarieties of $J_0(65)$, then taking the
product of one of the subvarieties of $J_0(65)$ with
the weight~$4$ motive $J_1(11,4)$.
\begincode
> J := JZero(65); 
> D := Decomposition(J); D;
[
    Modular abelian variety 65A of dimension 1, level 5*13 and 
    conductor 5*13 over Q,
    Modular abelian variety 65B of dimension 2, level 5*13 and 
    conductor 5^2*13^2 over Q,
    Modular abelian variety 65C of dimension 2, level 5*13 and 
    conductor 5^2*13^2 over Q
]
> A := D[1];
> B := D[2];
> A*B; 
Modular abelian variety 65A x 65B of dimension 3 and level 5*13 
over Q
Homomorphism from 65A to 65A x 65B given on integral homology by:
[1 0 0 0 0 0]
[0 1 0 0 0 0]
Homomorphism from 65B to 65A x 65B given on integral homology by:
[0 0 1 0 0 0]
[0 0 0 1 0 0]
[0 0 0 0 1 0]
[0 0 0 0 0 1]
Homomorphism from 65A x 65B to 65A (not printing 6x2 matrix)
Homomorphism from 65A x 65B to 65B (not printing 6x4 matrix)
> M := JZero(11,4);M;
Modular motive JZero(11,4) of dimension 2 and level 11 over Q
> P := A*M; P;
Modular motive 65A x JZero(11,4) of dimension 3 and level 5*11*13 over Q
\endcode
The product also returns inclusions of each factor into the
product and projection from the product onto each factor.
\begincode
> C,f,g := A*B; 
> f;
Homomorphism from 65A to 65A x 65B given on integral homology by:
[1 0 0 0 0 0]
[0 1 0 0 0 0]
\endcode

The command {\tt DirectSum}, which takes either two arguments,
or a sequence of modular abelian varieties, can be used to compute
arbitrary finite direct sums.
Note that the {\tt +} operator applied to two modular abelian
varieties is not the direct sum.  It is the sum of the two abelian
varieties in a common ambient abelian variety.  Thus if $A$ is
as above, then
\begincode
> Dimension(A);
1
> Dimension(A*A);
2
> Dimension(A+A);
1
\endcode

If you take a direct sum of abelian varieties that are defined over
different base rings, then the program will first attempt to
coerce them to a common over-ring.
\begincode
> A := JZero(11);
> B := BaseExtend(JZero(14),CyclotomicField(3));
> C := A*B; C;
Modular abelian variety JZero(11) x JZero(14) of dimension 2 and level 
2*7*11 over Q(zeta_3)
\endcode
The above would not work if {\tt CyclotomicField(3)}
were replaced by {\tt GF(3)}, since the base ring $\Q$ 
of~$A$ is not contained in {\tt GF(3)}.




\endex


END_HANDBOOK_EXAMPLES
*/
import "complements.m":
   FindHomologyDecompositionOfSubspaceUsingDecomposition;

import "endo_alg.m":
   Decide_If_Quaternionic;

import "elt.m": 
   Create_ModAbVarElt,
   Create_ModAbVarElt_Zero;

import "homology.m": 
   BasisChange_Rational_to_Lattice,
   Create_ModAbVarHomol,
   Create_ModAbVarHomol_Quotient,
   Create_ModAbVarHomol_Subspace,
   HomologySignZeroDimension;

import "inner_twists.m":
   Compute_Inner_Twists_Over_K;

import "linalg.m":
   HorizontallyStackMatrices,
   Intersect_Vector_Space_With_Lattice,
   LatticeDirectSum,
   MakeLattice,
   RestrictionOfScalars,
   VerticallyStackMatrices;

import "misc.m":
   FactorPrint,
   NoSpaces,
   OddPart;

import "morphisms.m": 
   Copy_MapModAbVar,
   Create_MapModAbVar,
   Create_MapModAbVar_MultiplyToMorphism;

import "rings.m": 
   CC, QQ, Qbar, RR, 
   CommonFieldOfDefinition,
   CommonBaseRing,
   RingName,
   IsDefnField,
   IsNumberField,
   OverCompositum,
   OverRing;

import "torsion.m":
   ComputeSize;

forward 
   AllCharactersTrivial,
   Can_Compute_Inner_Twists,
   Copy_ModAbVar,
   Create_ModAbVar,
   Create_ModAbVar_AssociatedToModularSymbolsSubspace,
   Create_ModAbVar_ModSym,
   InCreation,
   MStoHom_Coercion,
   Name,
   SetName,
   ShortDesc,
   Verbose;



declare verbose ModAbVar, 4;

declare attributes ModAbVar:
   base_ring,              // ring over which A is defined
   associated_abvar_over_Q, // if A is a base extension of something over Q, this is it.
   component_group_order,  // order of component group of A, when base char p
   conductor,		   // conductor of A as an abvar
   dimension,		   // dimension of A as an abvar
   endomorphism_ring,	   // ring over Z of endomorphisms of A
   endomorphism_algebra,   // algebra over Q of endomorphisms of A
   hecke_algebra,          // hecke algebra
   field_of_definition,	   // a minimal (as far as we know) field of defn of A
   frobpolys,              // cached charpolys of frobenius.
   homology,		   // the real, rational, and integral homology of A
   inner_twists,           
   weights,                // weights of modular forms spaces corr. to A
   num_rational_points,    // number of rational points.
   cm_twists,
   isogeny_decomp,         // full decomposition of A up to isogeny; a refinement
                           // of isogeny_newform_decomp. 
   tamagawa_numbers,       // Tamagawa numbers 
   decomposition,          // similar to isogeny_decomp, but subvarieties and different format.
   is_simple,	   // true <==> A is simple up to isogeny
   is_selfdual,	           // true if A is known to be self-dual, e.g., A = JZero(N), JOne(N), or JH(N).
   is_Af,                  // true if and only if A is an abvar of the form A_f.
   is_quaternionic,        
   is_only_motivic,        // attached to modular forms of weight bigger than 2.
   newform_number,
   newform, 
   modular_param,	   // if a surj. map from modular jacobian
			   // to A is known, it is stored here.
   modular_embedding,	   // if a finite-kernel map to a modular jacobian
			   // is known, it is stored here.
   modular_degree,
   embeddings,             // * A list of morphisms from A into abelian varieties,
                           //   which are used in making sense of intersections, sums, etc.
   only_maximal_newform_blocks,  // if in newform decomp of A the newform blocks are as big as possible.
   level,                  
   point_precision,        // used when deciding equality of points in real homology mod integral homology. 
   name;		   // a helpful name to remember what A is, e.g.,"JZero(37)"
   

/***************************************************************************

  << Creating the Modular Jacobian $J_0(N)$ >>

Use the {\tt JZero} command to create the Jacobian $J_0(N)$ of the
modular curve $X_0(N)$ (which parameterizes pairs consisting of an
elliptic curve and a cyclic subgroup of order $N$).  You can also
create higher weight motivic analogues of this Jacobian, and you can
compute in the $+1$ or $-1$ quotient of homology for efficiency,
though certain results will be off by factors of $2$.

EXAMPLES
> JZero(23);
Modular abelian variety JZero(23) of dimension 2 and level 23 over Q
> JZero(23 : Sign := +1);
Modular abelian variety JZero(23) of dimension 2 and level 23 over Q
with sign 1
> JZero(23,4);
Modular motive JZero(23,4) of dimension 5 and level 23 over Q
> JZero(23,4 : Sign := -1);
Modular motive JZero(23,4) of dimension 5 and level 23 over Q with 
sign -1
> JZero(389,2,+1);
Modular abelian variety JZero(389) of dimension 32 and level 389 
over Q with sign 1

 ***************************************************************************/

intrinsic JZero(N::RngIntElt : Sign := 0) -> ModAbVar
{Create the modular abelian variety J_0(N), i.e., the Jacobian
 of the modular curve X_0(N).}
   return JZero(N,2 : Sign := Sign);
end intrinsic;

intrinsic JZero(N::RngIntElt, k::RngIntElt, Sign::RngIntElt) -> ModAbVar
{Create the modular abelian variety J_0(N) of weight k>=2.}
   requirege k, 2;
   require Sign in {-1,0,1} : 
       "Optional sign must be -1, 0, or 1."; 
   return JZero(N,k:Sign := Sign);
end intrinsic;

intrinsic JZero(N::RngIntElt, k::RngIntElt: Sign := 0) -> ModAbVar
{Create the modular abelian variety J_0(N) of weight k>=2.}
   Verbose("JZero",Sprintf("Creating JZero(%o,%o,%o).",N,k,Sign),"");
   requirege k, 2;
   if IsOdd(k) then
      return ZeroModularAbelianVariety(k);
   end if;
   require Sign in {-1,0,1} : 
       "Optional sign must be -1, 0, or 1."; 
   name := Sprintf("JZero(%o%o)",N,k eq 2 select "" else Sprintf(",%o",k));
   J := Create_ModAbVar_ModSym(name, ModularSymbols(N,k,Sign));
   J`is_selfdual := true;
   J`weights := {k};
   return J;
end intrinsic;

/***************************************************************************

  << Creating the Modular Jacobians $J_1(N)$ and $J_H(N)$ >>

The {\tt JOne} command create the Jacobian of the modular curve
$X_1(N)$ (which parameterizes pairs consisting of an elliptic curve
and a point of order $N$).  The command {\tt Js} creates an abelian
variety isogenous to $J_1(N)$; more precisely it is the product of
abelian variteis $J_{\eps}(N)$, where $J_{\eps}(N)$ is the abelian
variety attached to all modular forms that have character a Galois
conjugate of $\eps$.  Creating $J_s(N)$ much faster than creating
$J_1(N)$, since less time is spent finding the integral structure on
homomology.

The  {\tt  JH} command  create  the  Jacobian  $J_H(N)$ of  the  curve
$X_H(N)$, which is the quotient of $X_1(N)$ by the subgroup $H$ of the
integers modulo $N$.

EXAMPLES

> JOne(13);
Modular abelian variety JOne(13) of dimension 2 and level 13 over Q
> JOne(13,4);
Modular motive JOne(13,4) of dimension 15 and level 13 over Q
> JOne(13,4 : Sign := 1);
Modular motive JOne(13,4) of dimension 15 and level 13 over Q
> JH(13,6);
Modular abelian variety J_H(13) of dimension 2 and level 13 over Q
> JH(13,3);
Modular abelian variety J_H(13) of dimension 0 and level 13 over Q
> JH(13,[-1]);
Modular abelian variety J_H(13) of dimension 2 and level 13 over Q
> JOne(17);
Modular abelian variety JOne(17) of dimension 5 and level 17 over Q
> Js(17);
Modular abelian variety Js(17) of dimension 5 and level 17 over Q
> IsIsogenous(JOne(17),Js(17));
true
> Degree(NaturalMap(JOne(17),Js(17)));
16
> JH(17,2);
Modular abelian variety J_H(17) of dimension 1 and level 17 over Q
> JH(17,4);
Modular abelian variety J_H(17) of dimension 1 and level 17 over Q
> JH(17,8);
Modular abelian variety J_H(17) of dimension 5 and level 17 over Q


 ***************************************************************************/

intrinsic JOne(N::RngIntElt : Sign := 0) -> ModAbVar
{}
   return JOne(N,2 : Sign := Sign);
end intrinsic;

intrinsic JOne(N::RngIntElt, k::RngIntElt, Sign::RngIntElt) -> ModAbVar
{}
   requirege k,2;
   return JOne(N,k : Sign := Sign);
end intrinsic;

intrinsic JOne(N::RngIntElt, k::RngIntElt: Sign := 0) -> ModAbVar
{Create the modular abelian variety J_1(N), i.e., the Jacobian
of the modular curve X_1(N).  Note that creating and finding
the integral structure on J_s(N), which is isogenous to J_1(N), 
is much faster.   Take great care in computing with J_1(N), since
it could be very slow.}
/* SOLUTION:
    Create JOne as below, then make *another* modular abelian variety 
    which has the correct lattice, but is not associated to modular 
    symbols.   That is, make prod J(eps) and construct JOne(N)
    as an abvar isogenous to J(eps) with this modular param.
*/

   Verbose("JOne",Sprintf("Creating JOne(%o,%o,%o).",N,k,Sign),"");

   requirege k,2;
   J := Js(N, k: Sign := Sign);
   if Dimension(J) eq 0 then
      name := Sprintf("JOne(%o)",N); 
      SetName(J,name);
      J`weights := {k};
      return J;
   end if;

   // Find sub-lattice L of H_1(J,Z) generated by images of integral
   // Modular symbols for Gamma1(N).  These are symbols of the
   // form X^i*Y^(k-2-i){0,gamma(0)}, where i=0,...,k-2 and
   // gamma runs through generators for Gamma1(N).  Fortunately,
   // Helena Verrill implemented an excellent algorithm for listing
   // a quite small generating set for Gamma1(N), so it is reasonable
   // to use this algorithm in practice!

   t := Cputime();
   // For some reason the following results 
   // in an error in the Version 2.11-BETA-1 of MAGMA:
   //R<X,Y> := PolynomialRing(QQ,2); 
   // so I instead use the following three lines, which is fine
   // for my application. (nothing is printed anyways).
   R := PolynomialRing(QQ,2);
   X := R.1;
   Y := R.2;
   Gamma1Gens := Generators(Gamma1(N));
   symbols := [ <X^i*Y^(k-2-i), [Cusps()|0,g*0]> : 
                               i in [0..k-2], g in Gamma1Gens];

   L := MakeLattice([ModularSymbolToRationalHomology(J,s) : s in symbols]);

   embed_morphism := IdentityMatrix(QQ,Dimension(Homology(J)));
   name := Sprintf("JOne(%o%o)",N,k eq 2 select "" else Sprintf(",%o",k));
   K := QQ;
   R := QQ;
   H := Create_ModAbVarHomol(L);
   J1N := Create_ModAbVar(name, H, K, R, <J, embed_morphism, false>);
   J1N`weights := {k};
   J1N`is_selfdual := true;
   return J1N;
end intrinsic;


intrinsic JH(N::RngIntElt, gens::[RngIntElt]: k:=2, Sign:=0) -> ModAbVar
{Let H be the subgroup of (Z/NZ)^* generated by gens. 
Create the modular abelian variety J_H(N), where eps is a Dirichlet character 
mod N with kernel H and J_H(N) is isogenous to the Jacobian of the modular 
curve X_H(N) associated to the subgroup of SL_2(Z) of matrices [a,b;c,d] with c
divisible by N and a in H modulo N.  It is the product of modular
symbols variety J(eps) for all Dirichlet characters eps that are trivial on H.}
   requirege N,1;
   require #[d : d in gens | GCD(N,d) ne 1] eq 0 : 
                   "Arguments 1 and 2 must be coprime.";
   require Type(k) eq RngIntElt and k ge 2 : 
       "Optional argument k must be an integer >= 2.";
   require Sign in {-1,0,1} : 
       "Optional sign must be -1, 0, or 1."; 

   F := CyclotomicField(EulerPhi(N));
   D := DirichletGroup(N,F);
   X := [ModularSymbols(MinimalBaseRingCharacter(eps),k,Sign) : 
                  eps in GaloisConjugacyRepresentatives(D) | 
                    #[d : d in gens | Evaluate(eps,d) ne 1] eq 0];
   name := Sprintf("J_H(%o)",N);
   J := Create_ModAbVar_ModSym(name, X);
   J`weights := {k};
   J`is_selfdual := true;
   return J;
end intrinsic;

intrinsic JH(N::RngIntElt, d::RngIntElt : k := 2,  Sign := 0) -> ModAbVar
{Let H be some subgroup of G=(Z/NZ)^* such that G/H has order d.
Create the modular abelian variety J_H(N), where eps is a Dirichlet
character mod N with kernel H and J_H(N) is *isogenous* to the
Jacobian of the modular curve X_H(N) associated to the subgroup of
SL_2(Z) of matrices [a,b;c,d] with c divisible by N and a in H modulo
N.  It is the product of modular symbols variety J(eps) for all
Dirichlet characters eps that are trivial on H.}
   requirege d, 1;
   require EulerPhi(N) mod d eq 0 : 
           "Argument 2(",d,") must divide",EulerPhi(N);
   ord := EulerPhi(N) div d;
   if ord eq 1 then
      return Js(N, k  :  Sign := Sign);
   elif ord eq EulerPhi(N) then
      return JZero(N, k : Sign := Sign);      
   end if;
   U,phi := UnitGroup(IntegerRing(N));
   for d in U do
      if Order(d) eq ord then
         return JH(N, [Integers()|phi(d)] : k := k, Sign := Sign);
      end if;
   end for;
   require false : "There is no cyclic subgroup of (Z/N Z)^* of order ord.";
end intrinsic;

intrinsic Js(N::RngIntElt, k::RngIntElt :  Sign := 0) -> ModAbVar
{A modular abelian variety that is Q-isogenous to the weight k version
of J_1(N).  More precisely, Js is  the direct sum of the modular abelian 
varieties attached to modular symbols spaces with Nebentypus.}
   requirege k, 2;
   name := Sprintf("Js(%o%o)",N,k eq 2 select "" else Sprintf(",%o",k));

   F := CyclotomicField(EulerPhi(N));
   D := DirichletGroup(N,F);
   X := [ModularSymbols(MinimalBaseRingCharacter(eps),k) : 
                        eps in GaloisConjugacyRepresentatives(D) | 
                        DimensionCuspForms(eps,k) ne 0 and 
                        (IsEven(eps) and IsEven(k)) or (IsOdd(eps) and IsOdd(k))];
   if #X eq 0 then
      return JZero(N, k : Sign := Sign);
   end if;
   J := Create_ModAbVar_ModSym(name, X);
   J`weights := {k};
   return J;
end intrinsic;


intrinsic Js(N::RngIntElt : Sign := 0) -> ModAbVar
{}
   return Js(N,2 : Sign := Sign);
end intrinsic;

/***************************************************************************

  << Abelian Varieties Attached to Modular Forms >>

The following commands create abelian varieties attached to spaces of modular
forms, sequences of spaces of forms, newforms, and characters.     If an input space
of modular forms is not cuspidal, MAGMA automatically replaces it with its
cuspidal subspace.

EXAMPLES
We first create the modular abelian variety attached to the spaces
$S_2(\Gamma_0(11))$ and $S_2(\Gamma_1(13))$.  This is the direct
sum of $J_0(11)$ with $J_1(13)$.
\begincode
> X := [ModularForms(11,2), ModularForms(Gamma1(13),2)];
> A := ModularAbelianVariety(X); A;
Modular abelian variety of dimension 3 and level 11*13 over Q
> IsIsomorphic(A, JZero(11)*JOne(13));
true Homomorphism N(1) from modular abelian variety of dimension 
3 to JZero(11) x JOne(13) (not printing 6x6 matrix)
\endcode
Next we create the modular abelian variety $A$ attached to
$S_2(\Gamma_1(17))$ along with $J_1(17)$ and $J_s(17)$.
We then note that $A$ is isomorphic to $J_1(17)$, but there
is no reason that $A$ should be isomorphic to $J_s(17)$ (they
are probably only isogenous).   This example illustrates
the fact that the abelian variety computed by MAGMA attached
to $S_2(\Gamma_1(17))$ is $J_1(17)$ rather than $J_s(17)$.
(Recall that $J_s(N)$ is a product of copies of abelian
varieties corresponding to conjugacy classes of characters.)
\begincode
> A := ModularAbelianVariety(ModularForms(Gamma1(17),2)); A;
Modular abelian variety of dimension 5 and level 17 over Q
> B := JOne(17); B;
Modular abelian variety JOne(17) of dimension 5 and level 17 over Q
> C := Js(17); C;
Modular abelian variety Js(17) of dimension 5 and level 17 over Q
> IsIsomorphic(A,B);
true Homomorphism from modular abelian variety of dimension 5 to 
JOne(17) (not printing 10x10 matrix)
> Degree(NaturalMap(A,C));
16
\endcode

If $\eps$ is a Dirichlet character and $k\geq 2$ is an integer,
let  $S$ be the space of modular forms with weight $k$ and
character a Galois conjugate of $\eps$.  The command
{\tt ModularAbelianVariety(eps,k)} computes the modular abelian
variety attached to~$S$.
\begincode
> G<eps> := DirichletGroup(13,CyclotomicField(12));
> Order(eps^2);
6
> ModularAbelianVariety(eps^2);
Modular abelian variety of dimension 2 and level 13 over Q
> ModularAbelianVariety(eps,3);
Modular motive of dimension 4 and level 13 over Q
\endcode
Next we compute the modular abelian variety attached
to a newform in $S_2(\Gamma_1(25))$.
\begincode
> S := CuspForms(Gamma1(25),2);
> N := Newforms(S);
> #N;
2
> f := N[1][1];
> PowerSeries(f,4);
q + a*q^2 + 1/1355*(941*a^7 + 4820*a^6 + 11150*a^5 + 11522*a^4 + 
    3582*a^3 + 10041*a^2 + 24432*a - 5718)*q^3 + O(q^4)
> A_f := ModularAbelianVariety(f);
> A_f;
Modular abelian variety Af of dimension 8 and level 5^2 over Q
\endcode
The abelian variety $A_f$ also determines the newform:
\begincode
> PowerSeries(Newform(A_f),4);
q + a*q^2 + 1/1355*(941*a^7 + 4820*a^6 + 11150*a^5 + 11522*a^4 + 
    3582*a^3 + 10041*a^2 + 24432*a - 5718)*q^3 + O(q^4)
\endcode
The {\tt Newform} command works even if $A$ wasn't explicitly created using a newform.
\begincode
> A := Decomposition(JZero(37))[1];
> Newform(A);
q - 2*q^2 - 3*q^3 + 2*q^4 - 2*q^5 + 6*q^6 - q^7 + O(q^8)
\endcode

 ***************************************************************************/

intrinsic Newform(A::ModAbVar) -> ModFrmElt
{A newform f so that A is isogenous to the newform abelian variety A_f.  It is
an error if A is not attached to a newform.}
   if not assigned A`newform then
      Verbose("Newform",Sprintf("Finding newform attached to %o.",A),"");
      t, Af := IsAttachedToNewform(A);
      require t : "Argument 1 must be attached to a newform.";
      e := DirichletCharacter(ModularSymbols(Af)[1]);
      S := CuspForms(e,Weight(ModularSymbols(Af)[1]));
      f := Newform(S,Af`newform_number);
      A`newform := f;
   end if;
   return A`newform;
end intrinsic;


intrinsic ModularAbelianVariety(f::ModFrmElt) -> ModAbVar
{The abelian variety attached to the newform f.}
   require IsNewform(f) : "Argument 1 must be a newform.";
   Verbose("ModularAbelianVariety",Sprintf("Finding modular abelian variety attached to %o.",f),"");
   A := Create_ModAbVar_AssociatedToModularSymbolsSubspace("Af", ModularSymbols(f));
   A`weights := {Weight(f)};
   A`newform := f;
   return A;
end intrinsic;

intrinsic ModularAbelianVariety(eps::GrpDrchElt : Sign := 0) -> ModAbVar
{}
   Verbose("ModularAbelianVariety",
     Sprintf("Modular abelian variety with character %o and sign %o.",eps, Sign),"");
   eps := MinimalBaseRingCharacter(eps);
   return ModularAbelianVariety(CuspidalSubspace(ModularSymbols(eps,2,Sign)));
end intrinsic;

intrinsic ModularAbelianVariety(eps::GrpDrchElt, k::RngIntElt : 
      Sign := 0) ->  ModAbVar
{The abelian variety associated to eps.  This corresponds to the space
of modular forms of weight k and character any Galois conjugate of eps.
We include all Galois conjugates in order to obtain an abelian variety
that is defined over Q.}
   eps := MinimalBaseRingCharacter(eps);
   Verbose("ModularAbelianVariety",
     Sprintf("Modular abelian variety of weight %o with character %o and sign %o.",k,eps, Sign),"");
   return ModularAbelianVariety(CuspidalSubspace(ModularSymbols(eps,k,Sign)));
end intrinsic;

intrinsic ModularAbelianVariety(X::[ModFrm] : Sign := 0) -> ModAbVar
{The abelian variety attached to the sequence X of modular forms spaces. 
This is the direct sum of the spaces attached to each element of the sequence.}
   Verbose("ModularAbelianVariety",
     Sprintf("Modular abelian variety associated to modular forms spaces %o.",X), "");
   if #X eq 0 then
      return ZeroModularAbelianVariety();
   end if;
   return ModularAbelianVariety([ModularSymbols(M, Sign): M in X]);
end intrinsic;

intrinsic ModularAbelianVariety(M::ModFrm : Sign := 0) -> ModAbVar
{The abelian variety attached to the modular forms space M.}
   Verbose("ModularAbelianVariety",
     Sprintf("Modular abelian variety associated to %o.",M),"");
   return ModularAbelianVariety(ModularSymbols(M, Sign));
end intrinsic;

/***************************************************************************

  << Abelian Varieties Attached to Modular Symbols >>

The commands below associated modular abelian varieties to
spaces of modular symbols and to sequences of spaces of modular symbols.
Conversely, the associate spaces of modular symbols to modular
abelian varieties.   If an input space of modular symbols is not
cuspidal, it is replaced by its cuspidal subspace.

EXAMPLES
We create modular abelian varieties attached to several spaces
of modular symbols.
\begincode
> M := ModularSymbols(37,2);
> ModularAbelianVariety(M);
Modular abelian variety of dimension 2 and level 37 over Q
> M := ModularSymbols(Gamma1(17));
> ModularAbelianVariety(M);
Modular abelian variety of dimension 5 and level 17 over Q
\endcode
Note that the sign of the space of modular symbols determines
the sign of the corresponding abelian variety.
\begincode
> M := ModularSymbols(Gamma1(17),2,+1);
> A := ModularAbelianVariety(M); A;
Modular abelian variety of dimension 5 and level 17 over Q with sign 1
\endcode
We can also create an abelian variety attached to any sequence 
of modular symbols spaces.
\begincode
> ModularAbelianVariety([ModularSymbols(11), ModularSymbols(Gamma1(13))]);
Modular abelian variety of dimension 3 and level 11*13 over Q
\endcode
Conversely, there is a sequence of modular symbols spaces
associated to any abelian variety defined over $\Q$.  These
need not be the same as the spaces used to define the 
modular abelian variety; instead they are what is used internally
in computatations on that abelian variety.
\begincode
> ModularSymbols(JOne(13));
[
    Modular symbols space of level 13, weight 2, character $.1, 
    and dimension 2 over Cyclotomic Field of order 6 and degree 2
]
> ModularSymbols(JZero(37));
[
    Modular symbols space for Gamma_0(37) of weight 2 and 
    dimension 4 over Rational Field
]
> A := JOne(17);
> ModularSymbols(A);
[
    Modular symbols space for Gamma_0(17) of weight 2 and 
    dimension 2 over Rational Field,
    Modular symbols space of level 17, weight 2, character $.1, 
    and dimension 2 over Cyclotomic Field of order 8 and degree 4
]
\endcode

 ***************************************************************************/

intrinsic ModularSymbols(A::ModAbVar) -> SeqEnum
{A sequence of spaces of modular symbols associated to A.}
   Verbose("ModularSymbols", 
     Sprintf("Modular symbols spaces attached to %o.",A),"");
   if IsAttachedToModularSymbols(A) then
      return (A`homology)`modsym;
   end if;
   if Type(BaseRing(A)) eq FldRat then
      return &cat[ModularSymbols(Af[1]) : Af in Factorization(A)];
   end if;
   require false :  "Argument 1 must be associated to modular symbols or defined over Q.";
end intrinsic;


intrinsic ModularAbelianVariety(M::ModSym) -> ModAbVar
{The abelian variety attached to the modular symbols space M.} 
   Verbose("ModularAbelianVariety", 
     Sprintf("Modular abelian variety attached to %o.",M),"");

   if not IsMultiChar(M) then
      eps := DirichletCharacter(M);   
      require BaseField(M) cmpeq BaseRing(Parent(MinimalBaseRingCharacter(eps))) :
           "Must have base field equal the field generated the character values.";
   end if;
   return Create_ModAbVar_AssociatedToModularSymbolsSubspace("",M);
end intrinsic;

intrinsic ModularAbelianVariety(X::[ModSym]) -> ModAbVar
{The abelian variety attached to the sequence X of modular 
 symbols spaces.}
   Verbose("ModularAbelianVariety", 
     Sprintf("Modular abelian variety attached to %o.",X),"");
   if #X eq 0 then
      return ZeroModularAbelianVariety();
   end if;
   sign := Sign(X[1]);
   require #Set([Sign(M) : M in X]) eq 1 : 
               "The signs must all be the same.";
   Y := [];
   has_multichar := &or[IsMultiChar(M) : M in X];
   if has_multichar and #X gt 1 then
      X1 := [M : M in X | not IsMultiChar(M)];
      X2 := [M : M in X | IsMultiChar(M)];
      A1 := ModularAbelianVariety(X1);
      A2 := [ModularAbelianVariety(M) : M in X2];
      A := DirectSum([A1] cat A2);
      return A;
   end if;

   for M in X do 
      require Sign(M) eq sign : 
            "The modular symbols spaces must all have the same sign.";
      require Type(BaseField(M)) in {FldRat, FldCyc} : 
            "Each modular symbols space must be defined over"*
             " the rationals or a cyclotomic field.";
      if IsMultiChar(M) then
         for MM in MultiSpaces(M) do
            Append(~Y,CuspidalSubspace(MM));
         end for;
      else
         eps := DirichletCharacter(M);   
         require BaseField(M) cmpeq 
              BaseRing(Parent(MinimalBaseRingCharacter(eps))) :
             "Each modular symbols space must have base field the " *
             "field generated the character values.";
         Append(~Y,CuspidalSubspace(M));
      end if;
   end for;
   name := "";
   if #Y eq 1 and Y[1] eq CuspidalSubspace(AmbientSpace(Y[1])) then
      eps := DirichletCharacter(Y[1]);
      if IsTrivial(eps) then
         name := Sprintf("JZero(%o)",Level(Y[1]));
      else
         name := Sprintf("J[%o,%o]", Level(Y[1]), NoSpaces(Sprintf("%o",Eltseq(eps))));
      end if;
   end if;
   A := Create_ModAbVar_ModSym(name, Y);
   if has_multichar then
      /* Because one of the modular symbols spaces is a multicharacter space,
         the lattice that defines A is not correct.  What we have to do
         is construct a B which has the correct lattice, which is the 
         sublattice of the lattice that defines A generated by integral
         basis for each homology space, and make the modular embedding the 
         natural inclusion map to A.

         Also, not that by what we did above, we may assume at this point
         that the sequence Y is just the factors of a single multicharacter
         space.
       */
       assert #X eq 1;
       t := Cputime();
       I := IntegralBasis(CuspidalSubspace(X[1]));
       G := [MStoHom_Coercion(A,MultiTuple(x)) : x in I];
       L := MakeLattice(G);
       embed_morphism := IdentityMatrix(QQ,Dimension(Homology(A)));
       K := QQ;
       R := QQ;
       H := Create_ModAbVarHomol(L);
       H`sign := sign;
       B := Create_ModAbVar(name, H, K, R, <A, embed_morphism, false>);
       B`weights := {Weight(X[1])};
       return B;

   end if;
   return A;
end intrinsic;


/***************************************************************************

  << Creation of Abelian Subvarieties >>

Suppose $A$ is an abelian variety and $V$ is a vector subspace of
the rational homology $H_1(A,\Q)$.    Then the {\tt DefinesAbelianSubvariety}
command determines whether or not $V$ is the rational homology of an abelian 
subvariety of $A$, and if so computes that abelian subvariety. 
The {\tt DefinesAbelianSubvariety} command relies on knowning
a complete decomposition of $A$ as a product of simple abelian
varieties (so it is currently restricted to abelian varieties
for which such a decomposition can be computed in MAGMA, e.g., all
modular abelian varieties over $\Q$).  

The other commands below are used to create the zero-dimensional 
abelian variety.

EXAMPLES
We define two subspaces of the rational homology of $J_0(33)$;
one defines an abelian subvariety and the other does not.
\begincode
> A := JZero(33);
> w3 := AtkinLehnerOperator(A,3);
> W := Kernel(Matrix(w3)+1);
> DefinesAbelianSubvariety(A,W);
true Modular abelian variety of dimension 1 and level 3*11 over Q
> V := RationalHomology(A);
> DefinesAbelianSubvariety(A,W + sub<V|[V.1]>);
false
\endcode

We create several zero-dimensional abelian varieties.
\begincode
> ZeroModularAbelianVariety();
Modular abelian variety ZERO of dimension 0 and level 1 over Q
> ZeroModularAbelianVariety(2);
Modular abelian variety ZERO of dimension 0 and level 1 over Q
> ZeroSubvariety(JZero(11));
Modular abelian variety ZERO of dimension 0 and level 11 over Q
\endcode

 ***************************************************************************/

intrinsic ZeroModularAbelianVariety(k::RngIntElt) -> ModAbVar
{The zero-dimensional abelian variety of weight k.}
   requirege k,2;
   A := ModularAbelianVariety(ZeroSubspace(ModularSymbols(1,k)));
   SetName(A,"ZERO");
   return A;
end intrinsic;

intrinsic ZeroModularAbelianVariety() -> ModAbVar
{The zero-dimensional abelian variety.}
   A := ZeroModularAbelianVariety(2);
   SetName(A,"ZERO");
   return A;
end intrinsic;

intrinsic ZeroSubvariety(A::ModAbVar) -> ModAbVar
{}
   Verbose("ZeroSubvariety",
     Sprintf("Creating zero subvariety of %o.",A),"");

   e := ModularEmbedding(A);
   J := Codomain(e);
   V := RationalHomology(A);
   W := sub<V|[V!0]>;
   H,embed := Create_ModAbVarHomol_Subspace(Homology(A),W);
   emb := embed*Matrix(e);
   Z := Create_ModAbVar("ZERO", H, FieldOfDefinition(A), BaseRing(A), 
                 <J, emb, Transpose(emb)>);
   return Z;
end intrinsic;

intrinsic DefinesAbelianSubvariety(A::ModAbVar, V::ModTupFld) -> 
                            BoolElt, ModAbVar
{True if and only if the subspace V of rational homology defines an abelian
subvariety of A.  If true, also returns the abelian subvariety.}
   require IsField(BaseRing(A)) : "Argument 1 must be over a field.";
   require Characteristic(BaseRing(A)) eq 0 : "Argument 1 must be over a field of characteristic 0.";
   require Type(BaseRing(V)) eq FldRat : "Argument 2 must be over the rational numbers.";
   require V subset RationalHomology(A) : 
              "Argument 1 must be contained in the rational homology of argument 1.";
   if Dimension(V) eq 0 then
      return ZeroSubvariety(A);
   end if;
      
   D, SUM, full_only := FindHomologyDecompositionOfSubspaceUsingDecomposition(A, V);
   if #D eq 0 then
      return false, _;
   end if;
   return true, SUM;
end intrinsic;


/***************************************************************************

  << Creation Using a Label >>

As a useful shorthand, it is sometimes possible to create modular
abelian varieties by giving a short string.  If the string contains a
single integer $N$, e.g., 37, then the corresponding abelian variety
is $J_0(N)$.  If it is of the form {\tt "<level>k<weight>"}, then it is
the possibly motivic $J_0(N)$ of weight~$k$.  If it is of the form
{\tt "<level>k<weight><isogeny code>"}, where {\tt <isogeny code>} is one
of {\tt "A"}, {\tt "B"}, ..., {\tt "Z"}, {\tt "AA"}, {\tt "BB"}, ..., {\tt "ZZ"},
{\tt "AAA"}, {\tt "BBB"}, ..., then the corresponding abelian variety is
{\tt JZero(N,k)(iso)}, where iso is a positive integer, and
{\tt "A"} corresponds to iso=1, {\tt "Z"} to iso=26, 
{\tt "AA"} to iso=27, {\tt "ZZ"} to iso=52, {\tt "AAA"} to 
iso=53, etc.

This labeling convention is the same as the one used for modular
symbols, and extends the one used for Cremona's database of elliptic
curves, except that Cremona's database contains some random scrambling
for levels between $56$ and $450$.  If the weight part of the label is
omitted, the weight is assumed to be~$2$.  To get the optimal quotient
of $J_0(N)$ with Cremona label $s$, set the optional parameter {\tt
Cremona} equal to true.

EXAMPLES

> ModularAbelianVariety("37");
Modular abelian variety 37 of dimension 2 and level 37 over Q
> ModularAbelianVariety("37A");
Modular abelian variety 37A of dimension 1 and level 37 over Q
> ModularAbelianVariety("11k4A");
Modular motive 11k4A of dimension 2 and level 11 over Q
> ModularAbelianVariety("65C");
Modular abelian variety 65C of dimension 2 and level 5*13 over Q
> ModularDegree(ModularAbelianVariety("56A"));
4
> ModularDegree(ModularAbelianVariety("56A" : Cremona := true));
2

 ***************************************************************************/

function CremonaCurve(s)
   return true, EllipticCurve(CremonaDatabase(),s);
/*   for i in [1..#s] do
      if s[i] in {"k", "K"} and i lt #s and not (s[i+1] in {"k","K"}) then
         if s[i+1] ne "2" then
            return false;
         end if;
     end if;
   end for;
*/
end function;

function CreateFromLabel(s, Cremona, sign)
   Verbose("ModularAbelianVariety", Sprintf(
       "Creating modular abelian variety from label '%o'.",s), 
       Sprintf("sign = %o, Cremona = %o", sign, Cremona));
   if Cremona then
      t, E := CremonaCurve(s);
      if t then
         return true, ModularAbelianVariety(E, sign);
      end if;
   end if;
   M := ModularSymbols(s, sign);
   A := Create_ModAbVar_AssociatedToModularSymbolsSubspace(s,M);
   return true, A;
end function;

intrinsic ModularAbelianVariety(s::MonStgElt : 
    Cremona := false  ) -> ModAbVar
{Abelian variety defined by string s.}
   t, A := CreateFromLabel(s, Cremona, 0);
   require t : "No abelian variety with label",s;
   return A;
end intrinsic;

intrinsic ModularAbelianVariety(s::MonStgElt, Sign::RngIntElt : 
    Cremona := false ) -> ModAbVar
{Abelian variety defined by string s.  See the documentation 
for more details.}
   require Sign in {-1,0,1} : 
       "Sign must be -1, 0, or 1."; 
   t, A := CreateFromLabel(s, Cremona, Sign);
   require t : "No abelian variety with label",s;
   return A;
end intrinsic;

/***************************************************************************

  << Creation Functions >>

  Each of these intrinsics creates an object of type {\tt ModAbVar},
i.e., a modular abelian variety.

The {\tt Sign} parameter is an integer that is either $0$ or $\pm 1$.
The reason for this parameter is that certain computations can be
done more quickly, at the expense of spurious powers of~$2$ in very
places.  If $A$ is a modular abelian variety
defined over $\R$, then $H_1(A,\Z)$ is equipped with an action
of complex conjugation~$c$.   
By working in the $+1$ or $-1$ eigenspaces for~$c$, many algorithms
that involve linear algebra with $H_1(A,\Z)$ run an order of
magnitude faster, since the spaces we have to consider have half
the dimension.  However, working in these eigenspaces will
produce incorrect results, e.g., most rational numbers will
be off by a power of~$2$, and if a group should be $G\oplus G$,
it will be output as $G'$, where $G'$ is isomorphic to $G$ away
from~$2$.  

EXAMPLES

\beginex{creation-basic}%\>--------------
We first create $J_0(46)$.
\begincode
> J := JZero(46); J;
Modular abelian variety JZero(46) of dimension 5 and level 2*23 over Q
\endcode
Next, we create the weight~$4$ analogue of $J_0(46)$, which is
a Grothendieck motive over~$\Q$.
\begincode
> J := JZero(46,4); J;
Modular motive JZero(46,4) of dimension 16 and level 2*23 over Q
\endcode
Now we create each of the Jacobians intermediate between $J_1(22)$
and $J_0(22)$.
\begincode
> JH(23,1);
Modular abelian variety JZero(23) of dimension 2 and level 23 over Q
> JH(23,2);
Modular abelian variety J_H(23) of dimension 2 and level 23 over Q
> JH(23,11);
Modular abelian variety J_H(23) of dimension 12 and level 23 over Q
> JH(23,22);
Modular abelian variety Js(23) of dimension 12 and level 23 over Q
\endex
Thus $J_1(23)$ has dimension $12$:
\begincode
> J := JOne(23);
> J;
Modular abelian variety JOne(23) of dimension 12 and level 23 over Q
\endcode
We can also create an abelian variety attached to a Dirichlet
character.
> A := ModularAbelianVariety(b); A;
Modular abelian variety of dimension 6 and level 2^2*5^2 over Q
\endcode
\endex


\beginex{creation-modular}%\>--------------
We can also create spaces associated to spaces of modular symbols
and modular forms.
\begincode
> M := ModularSymbols(50);
> A := ModularAbelianVariety(M); A;
Modular abelian variety of dimension 2 and level 2*5^2 over Q
\endcode
Notice that ModularAbelianVariety automatically replaced
$M$ by its cuspidal subspace.
\begincode
> B := ModularAbelianVariety(CuspidalSubspace(M)); B;
Modular abelian variety of dimension 2 and level 2*5^2 over Q
> A eq B;
true
\endcode
We can also construct the above modular symbols spaces
from modular forms.
\begincode
> S := CuspForms(50);
> C := ModularAbelianVariety(S);
> A eq C;
true
\endcode
We can also construct the modular abelian variety attached
to a newform.
\begincode
> N := Newforms(S); N;
[* [*
q - q^2 + q^3 + q^4 - q^6 + 2*q^7 + O(q^8)
*], [*
q + q^2 - q^3 + q^4 - q^6 - 2*q^7 + O(q^8)
*] *]
> f := N[1][1];
> D := ModularAbelianVariety(f); D;
Modular abelian variety Af of dimension 1 and level 2*5^2 over Q
> EllipticCurve(D);
Elliptic Curve defined by y^2 + x*y + y = x^3 - x - 2 over 
Rational Field
\endcode
In this example, the abelian variety has dimension~$1$, so there
is a corresponding elliptic curve.

\endex

 ***************************************************************************/

function Copy_ModAbVar(A)
   B := HackobjCreateRaw(ModAbVar);
   if assigned A`base_ring then B`base_ring := A`base_ring; end if;
   if assigned A`associated_abvar_over_Q then B`associated_abvar_over_Q := A`associated_abvar_over_Q; end if;
   if assigned A`component_group_order then B`component_group_order := A`component_group_order; end if;
   if assigned A`conductor then B`conductor := A`conductor; end if;
   if assigned A`dimension then B`dimension := A`dimension; end if;
   if assigned A`endomorphism_ring then B`endomorphism_ring := A`endomorphism_ring; end if;
   if assigned A`endomorphism_algebra then B`endomorphism_algebra := A`endomorphism_algebra; end if;
   if assigned A`hecke_algebra then B`hecke_algebra := A`hecke_algebra; end if;
   if assigned A`field_of_definition then B`field_of_definition := A`field_of_definition; end if;
   if assigned A`frobpolys then B`frobpolys := A`frobpolys; end if;
   if assigned A`homology then B`homology := A`homology; end if;
   if assigned A`inner_twists then B`inner_twists := A`inner_twists; end if;
   if assigned A`weights then B`weights := A`weights; end if;
   if assigned A`num_rational_points then B`num_rational_points := A`num_rational_points; end if;
   if assigned A`cm_twists then B`cm_twists := A`cm_twists; end if;
   if assigned A`isogeny_decomp then B`isogeny_decomp := A`isogeny_decomp; end if;
   if assigned A`tamagawa_numbers then B`tamagawa_numbers := A`tamagawa_numbers; end if;
   if assigned A`decomposition then B`decomposition := A`decomposition; end if;

   if assigned A`is_simple then B`is_simple := A`is_simple; end if;
   if assigned A`is_selfdual then B`is_selfdual := A`is_selfdual; end if;
   if assigned A`is_Af then B`is_Af := A`is_Af; end if;
   if assigned A`is_quaternionic then B`is_quaternionic := A`is_quaternionic; end if;
   if assigned A`is_only_motivic then B`is_only_motivic := A`is_only_motivic; end if;
   if assigned A`newform_number then B`newform_number := A`newform_number; end if;
   if assigned A`newform then B`newform := A`newform; end if;
   if assigned A`modular_param then B`modular_param := Copy_MapModAbVar(A`modular_param); end if;
   if assigned A`modular_embedding then B`modular_embedding := Copy_MapModAbVar(A`modular_embedding); end if;
   if assigned A`modular_degree then B`modular_degree := A`modular_degree; end if;
   if assigned A`embeddings then B`embeddings := A`embeddings; end if;
   if assigned A`only_maximal_newform_blocks then B`only_maximal_newform_blocks := A`only_maximal_newform_blocks; end if;
   if assigned A`level then B`level := A`level; end if;
   if assigned A`point_precision then B`point_precision := A`point_precision; end if;
   if assigned A`name then B`name := A`name; end if;

   return B;
end function;

function InCreation(A)
   if not assigned A`modular_embedding or not assigned A`modular_param then
      return true;
   end if;   
   return false;
end function;

function Create_ModAbVar(name, H, K, R, J)
   Verbose("Create_ModAbVar", 
         Sprintf("Creating %o of dimension %o", 
         name eq "" select "ModAbVar" else name, Dimension(H)), "");
   assert Type(name) eq MonStgElt;
   assert Type(H) eq ModAbVarHomol;
   assert IsDefnField(K);
   assert Type(J) in {Tup, BoolElt};
   // If J is a tup, then it must be a three-tuple: 
   //            < J, matrix_embed, matrix_param >,
   //  where J is a modular abelian variety such that
   //  IsAttachedToModularSymbols(J) is true, 
   //  matrix_embed defines a map from H1(A,Q) to H1(J,Q), 
   //  matrix_param defines a map from H1(J,Q) to H1(A,Q), 
   //  and the product of matrix_embed and matrix_param is
   //  multiplication by an integer.   These maps don't have
   //  to preserve integral homology -- they will be scaled 
   //  to do so.
   // ** If exactly one of matrix_embed or matrix_param
   //    is equal to the boolean false, then it is computed
   //    from the other using inverse morphisms

   A := HackobjCreateRaw(ModAbVar);
   A`name := name;
   A`homology := H;
   A`field_of_definition := K;
   A`base_ring := R;
   A`point_precision := 10;   // a default.
   if Type(J) eq BoolElt then  // identity map is modular param.
      assert assigned H`modsym;
      A`modular_param := nIsogeny(A,1);
      A`modular_embedding := nIsogeny(A,1);
   else
      assert Type(J) eq Tup;
      assert #J eq 3;
      J, matrix_embed, matrix_param := Explode(J);
      assert Type(J) eq ModAbVar;
      J := ChangeRing(J,R);
      assert {Type(matrix_embed), Type(matrix_param)}
                      subset {BoolElt, AlgMatElt, ModMatFldElt};
      if Type(matrix_embed) eq BoolElt then
     assert Type(matrix_param) ne BoolElt;
         assert Nrows(matrix_param) eq Dimension(Homology(J));
         assert Ncols(matrix_param) eq Dimension(Homology(A));
         A`modular_param := Create_MapModAbVar_MultiplyToMorphism(J, A, matrix_param, K);
         A`modular_embedding := LeftInverseMorphism(A`modular_param);
      elif Type(matrix_param) eq BoolElt then
         assert Type(matrix_embed) ne BoolElt;
         assert Nrows(matrix_embed) eq Dimension(Homology(A));
         assert Ncols(matrix_embed) eq Dimension(Homology(J));
         A`modular_embedding := Create_MapModAbVar_MultiplyToMorphism(A, J, matrix_embed, K);
         A`modular_param := RightInverseMorphism(A`modular_embedding);
      else
         A`modular_embedding := Create_MapModAbVar_MultiplyToMorphism(A, J, matrix_embed, K);
         assert Nrows(matrix_param) eq Dimension(Homology(J));
         assert Ncols(matrix_param) eq Dimension(Homology(A));
         A`modular_param := Create_MapModAbVar_MultiplyToMorphism(J, A, matrix_param, K);
      end if;
      // optimize the modular parameterization somewhat (i.e., if one is 2*phi, replace it with phi).
      // We are not making the parameterizations "optimal" in sense of connected kernel.
      A`modular_embedding := DivideOutIntegers(A`modular_embedding);
      A`modular_param := DivideOutIntegers(A`modular_param);
   end if;
   return A;
end function;



function Create_ModAbVar_ModSym(name, X) 
   Verbose("Create_ModAbVar_ModSym", "",
         Sprintf("Creating %o from %o", name, X));
   // The modular abelian variety associated to the spaces in the 
   // sequence X of  modular symbols.  The modular symbols
   // spaces need *not* be of weight 2.  If any space has 
   // weight bigger than 2, there is an abelian variety, but it
   // is only defined over the real numbers.    
   if Type(X) eq ModSym then
      X := [X];
   end if;
   assert Type(X) eq SeqEnum and #X gt 0 and Type(X[1]) eq ModSym;
   sign := Sign(X[1]);
   K := QQ;  // min field of definition
   is_only_motivic := false;
   for i in [1..#X] do
      M := X[i];     
      assert not IsMultiChar(M);
      if not IsCuspidal(M) then
         M := CuspidalSubspace(M);
         X[i] := M;
      end if;
      assert Sign(M) eq sign ;
      eps := DirichletCharacter(M);
      assert BaseField(M) cmpeq BaseRing(Parent(MinimalBaseRingCharacter(eps)));
      // as a motive it *is* over Q, and lots of stuff makes sense, so we'll 
      // just formally work as in weight 2 case, with same constraints...
      if Weight(M) ne 2 then 
         is_only_motivic := true;
      end if;
   end for;
   H := Create_ModAbVarHomol(X);
   A := Create_ModAbVar(name, H, K, K, true);
   A`is_only_motivic := is_only_motivic;
   return A;
end function;


function Create_ModAbVar_AssociatedToModularSymbolsSubspace(name, M);
   assert Type(name) eq MonStgElt;
   assert Type(M) eq ModSym;
   M := CuspidalSubspace(M);

   S := CuspidalSubspace(AmbientSpace(M));
   N := Level(M);
   J := ModularAbelianVariety([S]);

   /* Write subspace with respect to basis for S and restrict down to Q. */
   VS := VectorSpace(S);
   W := VectorSpace(BaseField(S),Dimension(S));
   V := VectorSpaceWithBasis([W|Coordinates(VS, b) : b in Basis(VectorSpace(M))]);
   V := RestrictionOfScalars(V);

   H,embed := Create_ModAbVarHomol_Subspace(Homology(J),V);
   A := Create_ModAbVar(name, H, QQ, QQ, <J, embed, false>);
   return A;
end function;



/***************************************************************************
 << Invariants >>

The invariants commands can be used to obtain the base ring, dimension,
character of defining modular form, a field of definition, the level,
the sign, the weights, and a short name of a modular abelian variety.

EXAMPLES
We illustrate all the commands for $J_0(23)$.
\begincode
> A := JZero(23);
> BaseRing(A);
Rational Field
> Dimension(A);
2
> DirichletCharacter(A);
1
> FieldOfDefinition(A);
Rational Field
> Level(A);
23
> Sign(A);
0
> Weights(A);
{ 2 }
\endcode

This is an example of a nontrivial Dirichlet character.
\begincode
> eps := DirichletCharacter(JOne(23)(2)); eps;
$.1^2
> Order(eps);
11
\endcode


We illustrate the {\tt Weights} command in several cases.
\begincode
> Weights(JZero(11));
{ 2 }
> Weights(JZero(11,4));
{ 4 }
> Weights(JOne(13,3));
{ 3 }
> Weights(DirectSum(JZero(11),JOne(13,3)));
{ 2, 3 }
> Weights(DirectSum(JZero(11),JZero(13,3)));
{ 2 }
\endcode

We display a few fields of definition.
\begincode
> FieldOfDefinition(JOne(13));
Rational Field
> FieldOfDefinition(BaseExtend(JZero(11),QuadraticField(7)));
Rational Field
> FieldOfDefinition(ChangeRing(JZero(11),GF(7)));
Finite field of size 7
\endcode
In the following example we quotient $J_0(11)$ out by a $5$-torsion point.
The resulting abelian variety might not be defined over $\Q$, and the
{\tt FieldOfDefinition} command currently plays it safe and returns 
$\Qbar$.
\begincode
> A := JZero(11);
> G := nTorsionSubgroup(A,5);
> H := Subgroup([G.1]);
> H;
Finitely generated subgroup of abelian variety with invariants [ 5 ]
> FieldOfDefinition(A/H);
Algebraically closed field with no variables
\endcode


 ***************************************************************************/

intrinsic Weights(A::ModAbVar) -> Set
{The set of weights of A.  (The weight need not be unique since direct sums of modular
symbols spaces of different weights are allowed.)}
   if not assigned A`weights then
      A`weights := Set(&cat[
          [Weight(M) : M in ModularSymbols(Af[1])] : Af in Factorization(A)]);
   end if;
   return A`weights;
end intrinsic;


intrinsic Level(A::ModAbVar) -> RngIntElt
{An integer N so that A is a quotient of a power of J_1(N).  Note
that N need not be minimal.  It is determined by how A is explicitly
represented as a quotient of modular Jacobians.}
   if not assigned A`level then
      if IsAttachedToModularSymbols(A) then
         A`level := LCM([Level(M) : M in A`homology`modsym]);
      else
         A`level := Level(Domain(ModularParameterization(A)));
      end if;
   end if;
   return A`level;
end intrinsic;
      

intrinsic FieldOfDefinition(A::ModAbVar) -> Fld
{The best known field of definition of A.}
   return A`field_of_definition;
end intrinsic;

intrinsic DirichletCharacter(A::ModAbVar) -> GrpDrchElt
{If A = A_f is attached to a newform, then this returns the Nebentypus
character of f.  Note that since f is only well-defined up to
Gal(Qbar/Q) conjugacy, the character is also only well-defined up
to Gal(Qbar/Q) conjugacy.}
   t, Af := IsAttachedToNewform(A);
   require t:"Argument 1 must be of the form A_f for a newform f.";
   M := Homology(Af)`modsym[1];
   return DirichletCharacter(M);
end intrinsic;

intrinsic DirichletCharacters(A::ModAbVar) -> List
{List of all Dirichlet characters of spaces of modular symbols associated with
the modular symbols abelian variety that parameterizes A.}
   J := Domain(ModularParameterization(A));
   M := ModularSymbols(J);
   ans := [* *];
   for m in M do 
      Append(~ans,DirichletCharacter(m));
   end for;
   return ans;
end intrinsic;


intrinsic Dimension(A::ModAbVar) -> RngIntElt
{The dimension of A.}
   if not assigned A`dimension then
      H := Homology(A);
      A`dimension := HomologySignZeroDimension(H);
   end if;
   return A`dimension;
end intrinsic;

intrinsic BaseRing(A::ModAbVar) -> Rng
{The ring that A is defined over.}
   return A`base_ring;
end intrinsic;


function ShortDesc(A) 
   assert Type(A) eq ModAbVar;
   if A`name eq "" then
      return Sprintf("modular abelian variety of dimension %o",Dimension(A));
   else
      return A`name;
   end if;
end function;


intrinsic Sign(A::ModAbVar) -> RngIntElt
{The sign of A, which is either 0, -1, or +1.  If +1 or -1, this means
we only compute the corresponding complex-conjugation eigenspace of
the homology of A, so various computations will be off by a factor of 2.}
   return Homology(A)`sign;
end intrinsic;


function Name(A)
   // A short string that describes A.
   // We skip type checking, so this function can apply to objects
   // of other types (e.g., when imported into decomp.m).
   return A`name;
end function;

procedure SetName(A, name) 
   // Set the name of A.
   A`name := name;
end procedure;


/***************************************************************************

 << Conductor >>

Let $A$ be a modular abelian variety over $\Q$.  The conductor command
computes the conductor of $A$ by factoring $A$ into newform abelian
varieties and using that the conductor of $A_f$ is $N^d$, where
$N$ is the level of $f$ and $d$ is the dimension of $A_f$.

EXAMPLES

> Factorization(Conductor(JZero(33)));
[ <3, 1>, <11, 3> ]
> Factorization(Conductor(JZero(11)^5));
[ <11, 5> ]
> Factorization(Conductor(OldSubvariety(JZero(46))));
[ <23, 4> ]
> Factorization(Conductor(JOne(25)));
[ <5, 24> ]


 ***************************************************************************/

intrinsic Conductor(A::ModAbVar) -> RngIntElt
{The conductor of the abelian variety A.  We require that A
is defined over Q.  When A=A_f is attached to a newform of level N,
then the conductor of A is N^d, where d is the dimension of A.}
   // TODO: Extend this for abelian varieties not necessarily over Q.
   // Is this easy?
   require BaseRing(A) cmpeq QQ : "Argument 1 must be defined over Q.";
   if not assigned A`conductor then
      if not IsAttachedToNewform(A) then
         D := Factorization(A);
         A`conductor := &*[Conductor(X[1])^#X[2] : X in D];
      else
         A`conductor := Level(A)^Dimension(A);
      end if;
   end if;
   return A`conductor;
end intrinsic;

/***************************************************************************

 << Number of Points >>

Given an abelian variety $A$ over a field $K$ the {\tt
NumberOfRationalPoints} or {\tt \#} command compute a divisor and
multiple of $\#A(K)$.   When finite, the multiple of the number of
rational points is computed using reduction mod primes up to 100.
Currently the lower bound is nontrivial only when $A$ is a quotient
of $J_0(N)$.

EXAMPLES

> #JZero(11);
5 5
> #JZero(23);
11 11
> #JZero(37);
Infinity Infinity
> #JOne(13);
1 19
> #JOne(23);
1 408991
> Factorization(408991);
[ <11, 1>, <37181, 1> ]
> NumberOfRationalPoints(ModularAbelianVariety("43B"));
7 7


 ***************************************************************************/

intrinsic NumberOfRationalPoints(A::ModAbVar) -> RngIntElt, RngIntElt
{Divisor and multiple of the cardinality of A(K),  where A is a modular abelian 
 variety defined over a field K.  If K is an abelian 
number field, then we assume the Birch and Swinnerton-Dyer 
conjecture. }
   return ComputeSize(A);
end intrinsic;

intrinsic '#'(A::ModAbVar) -> RngIntElt
{Same as NumberOfRationalPoints.}
   require IsAbelianVariety(A) : "Argument 1 must be an abelian variety.";
   return ComputeSize(A);
end intrinsic;



/***************************************************************************

  << Inner Twists and Complex Multiplication >>

If $f$ is a newform then an {\em inner twist} of $f$ is a Dirichlet
character $\chi$ such that the twist of $f$ by $\chi$ equals a Galois
conjugate of $f$, at least at Fourier coefficients whose subscript is
coprime to some fixed integer.  A {\em CM twist} is a nontrivial
character $\chi$ such that $f$ twisted by $\chi$ equals $f$, at least
at Fourier coefficients whose subscript is coprime to some fixed
integer.  The commands below find the CM and inner twists of the newform
corresponding to a newform abelian variety.

The optional parameter {\tt Proof} to each command is by default {\tt
false}.  If {\tt true}, it uses a huge number of terms of
$q$-expansions to ensure that that inner twist is really an inner
twist.  If {\tt false}, it uses far less (and is hence very quick),
and in practice this should be OK.


EXAMPLES
We compute the inner twists for $J_1(13)$.
\begincode
> A := JOne(13); A;
Modular abelian variety JOne(13) of dimension 2 and level 13 over Q
> CMTwists(A);
[]
> A2 := BaseExtend(A,AlgebraicClosure(RationalField()));
> CMTwists(A2);
[]
> InnerTwists(A2);
[
    1,
    $.1^5
]
> Parent($1[2]);
Group of Dirichlet characters of modulus 13 over Cyclotomic Field
of order 6 and degree 2
\endcode
We compute the inner twists for the second newform factor
of $J_1(23)$.
\begincode
> A := Decomposition(JOne(23))[2]; A;
Modular abelian variety image(23A[2]) of dimension 10, level 23 
and conductor 23^10 over Q
> InnerTwists(BaseExtend(A,AlgebraicClosure(RationalField())));
[
    1,
    $.1^20
]
\endcode
The CM elliptic curve $J_0(32)$ has a nontrivial CM inner twist.
\begincode
> A := JZero(32);
> InnerTwists(BaseExtend(A,AlgebraicClosure(RationalField())));
[
    1,
    $.1
]
> CMTwists(BaseExtend(A,AlgebraicClosure(RationalField())));
[
    $.1
]
\endcode
Quotients of $J_0(N)$ can also have nontrivial inner twists, which are
not CM twists.
\begincode
> J := JZero(81);
> A := Decomposition(J)[1];
> InnerTwists(BaseExtend(A,AlgebraicClosure(RationalField())));
[
    1,
    $.1
]
> CMTwists(BaseExtend(A,AlgebraicClosure(RationalField())));
[]
> Newform(A);
q + a*q^2 + q^4 - a*q^5 + 2*q^7 + O(q^8)
\endcode

The following is an example of a 4-dimensional abelian variety $A=A_f$ in 
$J_0(512)$ that has four inner twists, none of which are CM twists.
One can use this fact to prove that if $a_p$ is a prime-indexed Fourier
coefficient of $f$, then $a_p^2\in \Z$.  Thus no single $a_p$ generates
the degree $4$ field generated by all $a_n$.
\begincode
> J := JZero(512, 2, +1);
> A := Decomposition(J)[7]; A;
Modular abelian variety 512G of dimension 4, level 2^9 and 
conductor 2^36 over Q with sign 1
> f := Newform(A); f;
 q + 1/12*(a^3 - 30*a)*q^3 + 1/12*(-a^3 + 42*a)*q^5 +
     1/6*(-a^2 + 18)*q^7 + O(q^8)
> Coefficient(f,3)^2;
6
> Coefficient(f,5)^2;
12
> Coefficient(f,7)^2;
8
> Abar := BaseExtend(JZero(512,2,+1)(7),AlgebraicClosure(RationalField()));
> InnerTwists(Abar);
[
    1,
    $.1,
    $.2,
    $.1*$.2
]
> CMTwists(Abar);
[]
\endcode

 ***************************************************************************/

intrinsic CMTwists(A::ModAbVar : Proof := false) -> SeqEnum
{The CM inner twists characters of the newform abelian variety A = A_f 
that are defined over the base ring of A.  For all CM twists, base
extend to AlgebraicClosure(RationalField()) first.}
   if not assigned A`cm_twists then
      Verbose("CMTwists", Sprintf("Computing CM twists of %o.", A),
              Sprintf("Proof = %o", Proof));
      t, A0 := CanChangeRing(A,QQ);
      require t : "Argument 1 must be base extension of abelian variety over Q.";
      t, A0 := IsAttachedToNewform(A0);
      require t : "Argument 1 must be attached to a newform.";
      dummy := InnerTwists(A : Proof := Proof);
      require assigned A`cm_twists : "Argument 1 must be attached to a newform.";
   end if;
   return A`cm_twists;
end intrinsic;

intrinsic InnerTwists(A::ModAbVar : Proof := false) -> SeqEnum
{The inner twists characters of the newform abelian variety A = A_f 
that are defined over the base ring of A. For all CM twists, base
extend to AlgebraicClosure(RationalField()) first.}

   K := BaseRing(A);
   require IsField(K) : "Argument 1 must be over a field.";
   require Characteristic(K) eq 0 : "Argument 1 must be over a field of characteristic 0.";
   if not assigned A`inner_twists then
      Verbose("InnerTwists", Sprintf("Computing inner twists of %o.", A),
              Sprintf("Proof = %o", Proof));
      if Type(K) eq FldRat then
         A`inner_twists := [];
         A`cm_twists := [];
         return [];
      end if;
      t, A0 := CanChangeRing(A,QQ);
      require t : "Argument 1 must be base extension of abelian variety over Q.";
      t, A0 := IsAttachedToNewform(A0);
      require t : "Argument 1 must be attached to a newform.";
       A`inner_twists, A`cm_twists := 
            Compute_Inner_Twists_Over_K(ModularSymbols(A0)[1],BaseRing(A), Proof);         
   end if;
   return A`inner_twists;
end intrinsic;



/***************************************************************************

  << Predicates >>

Most of the predicates below work in full generality.  The ones whose
domain of applicability is somewhat limited are {\tt IsIsomorphic},
{\tt IsQuaternionic}, and {\tt IsSelfDual}.  The {\tt IsIsomorphic}
command will at least definitely determine whether any two simple
modular abelian varieties over $\Q$ are isomorphic.  In theory, it is
possible to determine isomorphism for more general classes of modular
abelian varieties, but this has not been implemented.


EXAMPLES
We test whether a few objects are actually abelian varieties.
\begincode
> A := JZero(11);
> A11 := ChangeRing(A,GF(11));
> IsAbelianVariety(A11);
false
> AZ := ChangeRing(A,Integers());
> IsAbelianVariety(AZ);
false
> A3 := ChangeRing(A,pAdicRing(3));
> IsAbelianVariety(A3);
true
\endcode
A modular motive is sometimes also an abelian variety, but only over
the complex numbers. 
\begincode
> A := JZero(11,4); A;
Modular motive JZero(11,4) of dimension 2 and level 11 over Q
> IsAbelianVariety(A);
false
> IsOnlyMotivic(A);
true
> IsAbelianVariety(BaseExtend(A,ComplexField()));
true
\endcode

Abelian varieties $J_0(N)$ and $J_s(N)$ are attached to modular symbols, as are
newform abelian varieties.
\begincode
> J := JZero(37);
> IsAttachedToModularSymbols(J);
true
> A := Decomposition(J)[1];
> IsAttachedToModularSymbols(A);
false
> t, Af := IsAttachedToNewform(A);
> IsAttachedToModularSymbols(Af);
true
> IsIsomorphic(A,Af);
true Homomorphism from 37A to 37A given on integral homology by:
[1 0]
[0 1]
> IsAttachedToModularSymbols(Js(17));
true
> IsAttachedToModularSymbols(JOne(17));
false
\endcode

We test isogeny between a few abelian varieties.
\begincode
> IsIsogenous(JZero(11),JOne(11));
true
> IsIsogenous(JZero(11)*JZero(11),JZero(22));
true
> IsIsogenous(JZero(11)*JZero(11),JZero(33));
false
> IsIsogenous(JZero(11),JZero(14));
false
> IsIsogenous(JZero(11)^2,JZero(22));
true
> A := JZero(37)(2); B := JOne(13);
> IsIsogenous(A*B*A, A*A*B);
true
> A := JZero(43);
> G := RationalCuspidalSubgroup(A);
> IsIsogenous(A,A/G);
true
\endcode

Next we test isomorphism between some abelian varieties.  
\begincode
> A := JZero(43);
> G := RationalCuspidalSubgroup(A);
> B := A/G;
> CanDetermineIsomorphism(A,B);
true false
> IsIsomorphic(A,B);
false
> IsIsomorphic(JZero(11),JOne(11));
false
> IsIsomorphic(JZero(13),JOne(13));
false
> IsIsomorphic(Js(13),JOne(13));
true Homomorphism from Js(13) to JOne(13) given on integral 
homology by:
[ 0  0 -1  0]
[ 0  0  0 -1]
[ 1  0 -1  0]
[ 0  1  0 -1]
> CanDetermineIsomorphism(Js(17),JOne(17));
false All tests failed to decide whether A and B are isomorphic.
\endcode

We test whether certain Jacobians have simple factors
with multiplicity one.
\begincode
> HasMultiplicityOne(JZero(43));
true
> HasMultiplicityOne(JZero(33));
false
> Decomposition(JZero(33));
[
    Modular abelian variety 33A of dimension 1, level 3*11 and 
    conductor 3*11 over Q,
    Modular abelian variety N(11,33,1)(11A) of dimension 1, level
    3*11 and conductor 11 over Q,
    Modular abelian variety N(11,33,3)(11A) of dimension 1, level
    3*11 and conductor 11 over Q
]
> 
\endcode

In the next example we give an example of a non-quaternionic
surface.  
\begincode
> IsQuaternionic(JOne(13));
false
\endcode

Jacobians are self dual, and there is a surface of level $43$ that is self dual and
a surface of level $69$ that is not.
\begincode
> IsSelfDual(JOne(13));
true
> IsSelfDual(JZero(69)(2));
false
\endcode

The surface $A$ below is isomorphic to its dual.
The natural polarization has kernel that is the kernel
of multiplication by $2$.
\begincode
> A  := JZero(43)(2);
> A;
Modular abelian variety 43B of dimension 2, level 43 and 
conductor 43^2 over Q
> IsSelfDual(A);
true
> phi := ModularPolarization(A);
> Invariants(Kernel(phi));
[ 2, 2, 2, 2 ]
\endcode

We test a few abelian varieties for simplicity.
\begincode
> IsSimple(JOne(25));
false
> IsSimple(JOne(13));
true
> IsSimple(JZero(11)^10);
false
> IsSimple(NewSubvariety(JZero(100)));
true
\endcode


 ***************************************************************************/

intrinsic IsQuaternionic(A::ModAbVar) -> BoolElt
{True if and only if some simple factor of A over the base ring 
has quaternionic multiplication.}
/*TODO: This will currently always return false, since the
function decomp.m:Create_Af_and_Homs always sets it to
false, since I haven't programmed deciding whether multiplication
is quaternionic in general.
*/
   if not assigned A`is_quaternionic then
      A`is_quaternionic := Decide_If_Quaternionic(A);
   end if;
   return A`is_quaternionic;
end intrinsic;

intrinsic IsSimple(A::ModAbVar) -> BoolElt
{True if and only if A has no proper abelian subvarieties over
 BaseRing(A).}
   if not assigned A`is_simple then
      X := Factorization(A);
      A`is_simple := #X eq 1 and #X[1][2] eq 1;
   end if;
   return A`is_simple;
end intrinsic;

intrinsic IsAttachedToNewform(A::ModAbVar) -> BoolElt, ModAbVar, MapModAbVar
{True if A is isogenous to a newform abelian variety A_f.  
This intrinsic also returns the abelian variety A_f.
Third return argument is explicit isogeny from A_f to A.}
   if not assigned A`is_Af then
      if not assigned A`is_Af then
         D := Factorization(A);
         if #D eq 0 then
            A`is_Af := false;
         elif #D gt 1 or #D[1][2] gt 1 then
            A`is_Af := false;
         else
            A`is_Af := <D[1][1], D[1][2][1]>;
         end if;
      end if;
   end if;
   if Type(A`is_Af) eq BoolElt then
      if A`is_Af then
         return true,  A, IdentityMap(A);
      else
         return false, _, _;
      end if;
   end if;
   return true, A`is_Af[1], A`is_Af[2];
end intrinsic;

intrinsic CanDetermineIsomorphism(A::ModAbVar, B::ModAbVar) -> 
BoolElt, BoolElt, MapModAbVar, RngIntElt
{True if we can determine whether or not A and B are isomorphic.  
If we can determine isomorphism, also returns true if A and B are isomorphic 
and an explicit isomorphism, or false if they are not isomorphic. 
If we can not determine isomorphism, also returns the reason why we can not
as a string.  If A and B are simple and defined over Q, then it is always
possible to determine whether A and B are isomorphic.  If one of A or B
has simple factors of multiplicity one, then in principal it is possible,
but the algorithm has not been programmed.}

   print "Tseno's version";

   require Sign(A) eq Sign(B) : "Arguments 1 and 2 must have the same sign.";
   Verbose("CanDetermineIsomorphism", Sprintf("Isomorphism between %o and %o?",A, B),""); 
   // print "Isogenous?";
   if not IsIsogenous(A,B) then
      Verbose("CanDetermineIsomorphism", "","   * Not even isogenous.");
      return true, false, "", 1;
   end if;
   if A eq B then
      Verbose("CanDetermineIsomorphism", "","   * Equal.");
      return true, true, IdentityMap(A), 2;
   end if;
   // print "Simple?";
   if IsSimple(A) and IsSimple(B) then
      Verbose("CanDetermineIsomorphism", "","   * Both simple.");
      f := NaturalMap(B,A);
      deg := Degree(f);
      Verbose("CanDetermineIsomorphism", "",  Sprintf("* Natural map has degree",deg));
    
      is_field, K, K_to_End, End_to_K := IsField(BaseExtend(End(A),QQ));
      if not is_field then
         return false, false, "", 3; 
      //      "Isomorphism testing not programmed when endomorphism rings are not fields", _;
      end if;

      // The endomorphism algebra is a field. 
      if Degree(K) eq Dimension(A) then
         t, d := IsSquare(deg);
         if not t then   // Since deg(K) = Dim(A), the
                         // degrees of endomorphisms are perfect squares
                         // since they are the squares of norms.
                         // This uses that the dimension of the endomorphism
                         // algebra equals the dimension of A, and that
                         // the endomorphism algebra is a field.
                         // (See Page 126, Prop 12.12 of Milne's 
                         // "Abelian Varieties" in [Cornell-Silverman].)
            return true, false, "", 4;
         end if;
      end if;

      // There is an isomorphism if and only if Hf contains
      // an element of degree deg.  We try to decide this by mapping
      // to the field Hom_0(A,A) and solving a norm equation.
      H := Hom(A,B);    // all homomorphism from A to B.  
      Hf := Pullback(H, f);  // subspace of End(A,A).


      Verbose("CanDetermineIsomorphism", "", 
              Sprintf("Looking for an element of norm", d));
      // Find all elements with norm d in K.
      if Type(K) eq FldRat then
         X := [d, -d];
         t := true;
      else
         // These two lines are commented out, since they would only be
         // useful if End(A) were maximal, and it usually isn't.
         // Fortunately, MAGMA solves Diophantine norm equations in
         // orders, so we don't have a problem. 
         // O_K := MaximalOrder(K);
         // t, X := NormEquation(O_K,d : All := true);
         G := [End_to_K(f) : f in Generators(End(A))];
         // print "Computing order";
         O := Order(G);
         // print "Solving norm equation";
         t, X := NormEquation(O,d : All := true);
         t1, Y := NormEquation(O,-d : All := true);
         // print "Solving done";
         if t then
            X := [K!a : a in X];
         end if;
	 if t1 then
	    for y in Y do
		Append(~X, K!y);
	    end for; 
	 end if;
      	 // print "Appending done";
      end if;
      // This is changed because of +/- business.
      if not (t or t1) then
         return true, false, "", 5;
      end if;
      // Now lift each of the elements of X back to the endomorphism
      // ring of A, and check to see if any lie in Hf.  One does if
      // A is isomorphic to B.
      for x in X do
         // print x;
         g := K_to_End(x);
         // print g;
         if g in Hf then   // g has degree deg and g = psi*f for some
                           // psi in Hom(A,B) and f has degree deg, 
                           // so psi=g*f^(-1) is an explicit isomorphism.
            // print "Good solution";
            return true, true, g*f^(-1), 6;
         end if;
      end for;
      // print "None of the solution worked";
      return true, false, "", 6;
   end if;
   f := NaturalMap(A,B);
   if IsIsogeny(f) then
      G := ComponentGroupOfKernel(f);
      if #G eq 1 then
         return true, true, f;
      end if;
      I := Invariants(G);
      if Sign(A) ne 0 then
         I := [Integers() |OddPart(n) : n in I | OddPart(n) ne 1];
      end if;
      if #I eq 2*Dimension(A) and #Set(I) eq 1 then
         return true, true, (1/I[1])*f;
      end if;
      /* We generalize prop 12.12 from Milne's article again (see above)
         to deduce that the degree of f must be a square, so long as 
            1) each factor in a simple decomposition of A has multiplicity 1.  
            2) the endomorphism algebra of each simple factor is a field
               of dimension equal to the degree of the factor.
         The proof is as follows.  Suppose, for simplicity that 
                              g:A-->C x D
         is an isogeny and that C and D are
         non-isogenous simple abelian varieties.  Suppose also that
         End_0(C) and End_0(D) are fields of degree equal to the dimension
         of C and D, respectively.  Let h be the complementary morphism
         to g, so that g*h = h*g = [deg(g)].   If phi:C --> C is a surjective
         endomorphism, then deg(h*phi*g) = deg(phi)*deg(g)^(2*dim(A)).
         Since h*phi*g is an endomorphism of C x D, it is a cross product of
         an endomorphism of C and an endomorphism of D.  As such, its
         degree is the product of two squares (by Milne), hence a square.
         Since deg(g)^(2*dim(A)) is a square, it follows that phi must 
         have degree a square as well.
      */
      if HasMultiplicityOne(A) then
         small_fields := true;
         if Type(BaseRing(A)) ne FldRat then
            for F in Decomposition(A) do 
               is_field, K, K_to_End, End_to_K := IsField(BaseExtend(End(F[1]),QQ));
               if not is_field or is_field and Degree(K) ne Dimension(F[1]) then
                  small_fields := false;
                  break;
               end if;
            end for;
         end if;
         if small_fields then
            t, _ := IsSquare(Degree(f));
            if not t then
               return true, false, _;
            end if;
         end if;
      end if;
   end if;
   return false, "All tests failed to decide whether A and B are isomorphic.", _;
end intrinsic;

intrinsic HasMultiplicityOne(A::ModAbVar) -> BoolElt
{True if the simple factors of A appear with multiplicity one.}
   for F in Factorization(A) do
      if #F[2] gt 1 then
         return false;
      end if;
   end for;
   return true;
end intrinsic;

intrinsic IsIsomorphic(A::ModAbVar, B::ModAbVar) -> BoolElt, MapModAbVar, RngIntElt
{True if A and B are isomorphic.  If true, also returns an explicit
isomorphism.  This command will work if A and B are defined over Q and
the simple factors occur with multiplicity one, and may work
otherwise, but it may terminate with an error in the general case.
Use the command CanDetermineIsomorphism to avoid getting
an error.}

   t, iso, map, reas := CanDetermineIsomorphism(A,B);
   require t : iso;
   if not iso then
      return iso, "", reas;
   end if;
   return iso, map, reas;
end intrinsic;


intrinsic MinimalIsogeny(A::ModAbVar, B::ModAbVar) -> BoolElt,
MapModAbVar, RngIntElt {True if A and B are isogenous. In that case
returns the minimal degree of an isogeny between A and B and an
explicit isogeny of that degree. False if A and B are not isogenous.
If A and B are simple and defined over Q, then it is always possible
to the minimal degree of an isogeny.}

   require Sign(A) eq Sign(B) : "Arguments 1 and 2 must have the same sign.";
   Verbose("MinimalIsogeny", Sprintf("Minimal isogeny between %o and %o?",A, B),""); 
   
   if not IsIsogenous(A,B) then
      Verbose("MinimalIsogeny", "","   * Not even isogenous.");
      return false, "", 0;
   end if;
   if A eq B then
      Verbose("MinimalIsogeny", "","   * Equal.");
      return true, IdentityMap(A), 1;
   end if;
   if IsSimple(A) and IsSimple(B) then
      Verbose("MinimalIsogeny", "","   * Both simple.");
      f := NaturalMap(B,A);
      deg := Degree(f);
      a, b := SquareFree(deg);
      Verbose("MinimalIsogeny", "",  Sprintf("* Natural map has degree",deg));
    
      is_field, K, K_to_End, End_to_K := IsField(BaseExtend(End(A),QQ));
      if not is_field then
         return false, "", 0; 
      //      "Minimal isogeny testing not programmed when endomorphism rings are not fields", _;
      end if;

      // The endomorphism algebra is a field. 

      i := 0;
      while true do	
         i := i+1;

         // There is an isogeny of degree a*i^2 if and only if Hf contains
         // an element of degree a^2*b^2*i^2.  We try to decide this by mapping
         // to the field Hom_0(A,A) and solving a norm equation.
         H := Hom(A,B);    // all homomorphism from A to B.  
         Hf := Pullback(H, f);  // subspace of End(A,A).

         Verbose("MinimalIsogeny", "", 
                 Sprintf("Looking for an element of norm", a*b*i));
         // Find all elements with norm a*b*i in K.
         G := [End_to_K(f) : f in Generators(End(A))];
         O := Order(G);
         t, X := NormEquation(O,a*b*i : All := true);
         t1, Y := NormEquation(O,-a*b*i : All := true);
         if t then
            X := [K!a : a in X];
         end if;
   	 if t1 then
   	    for y in Y do
   		Append(~X, K!y);
   	    end for; 
   	 end if;

         if t or t1 then
            // Now lift each of the elements of X back to the endomorphism
            // ring of A, and check to see if any lie in Hf.  One does if
            // A is isogenous B with an isogeny of degree a*i^2.
            for x in X do
               g := K_to_End(x);
               if g in Hf then   // g has degree a^2*b^2*i^2 and g = psi*f for some
                                 // psi in Hom(A,B) and f has degree a*b^2, 
                                 // so psi=g*f^(-1) is an explicit isogeny of degree a*i^2.
                  return true, g*f^(-1), a * i^2;
               end if;
            end for;
         end if;
      end while;  
   end if;
end intrinsic;


intrinsic IsIsogenous(A::ModAbVar, B::ModAbVar) -> BoolElt
{True if A and B are isogenous.  If this can *not* be determined,
then an error message is displayed and the program terminates.
It is always possible to determine whether or not A and B
are isogenous when both are defined over Q.}
   require BaseRing(A) cmpeq BaseRing(B) : 
        "Arguments 1 and 2 must be over the same base field.";
   if Dimension(A) ne Dimension(B) then
      return false;
   end if;
   if not IsSimple(A) then
      if IsSimple(B) then
         return false;
      end if;
      // Neither isogeny simple.
      DA := Factorization(A);
      DB := Factorization(B);
      // Pair off factors.
      for C in DA do 
         fail := true;
         for i in [1..#DB] do 
            if IsIsogenous(DB[i][1],C[1]) then 
               if #DB[i][2] ne #C[2] then
                  return false;
               end if;
               Remove(~DB,i);
               fail := false;
               break;
            end if;
         end for;
         if fail then
            return false;
         end if;
      end for;
      return true;
   end if;
   if not IsSimple(B) then
      return false;
   end if;
   // Both are isogeny simple. 
   A := Factorization(A)[1][1];
   B := Factorization(B)[1][1];
   if BaseRing(A) cmpeq QQ then
      // Since base is Q, these are newform abelian varieties, 
      // so can decide isogeny using degeneracy maps.
      f := NaturalMap(A,B,1);   
      return not IsZero(f);
   end if;
   if IsAttachedToModularSymbols(A) and IsAttachedToModularSymbols(B) then
      return ModularSymbols(A)[1] eq ModularSymbols(B)[1];
   end if;
   require false : "Unable to determine whether or not A and B are isogenous.";
end intrinsic;

intrinsic IsAbelianVariety(A::ModAbVar) -> BoolElt
{True if A is an abelian variety, i.e., defined over a ring of 
 characteristic 0 in which the conductor is invertible, or 
 a finite field that does not divide the conductor
 of A.  For example, if A has positive dimension and is defined over Z,
 then Raynaud's theorem implies that A is not an abelian variety.}
   if Dimension(A) eq 0 then
      return true;
   end if;
   if IsOnlyMotivic(A) then
      return Type(BaseRing(A)) eq Type(ComplexField());
   end if;
   R := BaseRing(A);
   n := Characteristic(R);
   if n eq 0 then
      if Type(R) in {FldRat, FldNum, FldCyc, FldQuad} then
         return true;
      end if;
   end if;
   t, A0 := CanChangeRing(A,QQ);
   require t : "It must be possible to change the base ring of " * 
               "argument 1 to the RationalField().";
   N := Conductor(A0);
   t := IsInvertible(R!N);
   return t;
end intrinsic;


intrinsic IsSelfDual(A::ModAbVar) -> BoolElt
{True if A is known to be isomorphic to its dual.  There is
an error message if MAGMA is unable to decide.}
   if not assigned A`is_selfdual then
      if Dimension(A) eq 1 then
         A`is_selfdual := true;
      else
         t, dual := CanComputeDual(A);
         require t : "Unable to determine the dual of A ("*dual*").";
         t, ans := CanDetermineIsomorphism(A,dual);
         require t : "Unable to determine if A is isomorphic to its dual ("*ans*").";
         if t then
            A`is_selfdual := ans;
         end if;
      end if;
   end if;
   return A`is_selfdual;
end intrinsic;

intrinsic IsAttachedToModularSymbols(A::ModAbVar) -> BoolElt
{True if the underlying homology of A is being computed using a space
of modular symbols.  For example, this will be true for J_0(N) 
and for newform abelian varieties.}
   return assigned (A`homology)`modsym;   
end intrinsic;

intrinsic IsOnlyMotivic(A::ModAbVar) -> BoolElt
{True if any of the modular forms attached to A have weight bigger than 2.}
   if not assigned A`is_only_motivic then
      if Type(BaseRing(A)) eq FldPr then
         A`is_only_motivic := false;
      else
         Amb := Codomain(ModularEmbedding(A));
         // TODO: Polish is_only_motivic
         // This is potentially a BUG, though definitely not serious.
         // (It might just print "abelian variety" when it should
         // print "motivic".)
         // At worst it would say that something is not motivic when
         // it is --- since we never do anything such that this matters,
         // it isn't an issue for now.  This should be improved
         // for the final release.
         A`is_only_motivic := assigned Amb`is_only_motivic and
                           Amb`is_only_motivic;
      end if;
   end if;
   return A`is_only_motivic;
end intrinsic;


/***************************************************************************

  << Printing >>

 ***************************************************************************/

intrinsic HackobjPrintModAbVar(A::ModAbVar, level::MonStgElt)
{}

   if InCreation(A) then
      printf "Modular abelian variety.";
   else
      printf "Modular %o%o of dimension %o%o level %o %oover %o%o", 
         IsOnlyMotivic(A) select "motive" else "abelian variety",
         A`name ne "" select " "*A`name else "", 
         Dimension(A), 
         assigned A`conductor select "," else " and",
         FactorPrint(Level(A)), 
         assigned A`conductor select Sprintf("and conductor %o ",FactorPrint(Conductor(A))) else "",
         RingName(BaseRing(A)), 
         Sign(A) ne 0 select Sprintf(" with sign %o",Sign(A)) else "";
   end if;
end intrinsic;


/***************************************************************************

  << Equality and Inclusion Testing >> 

These functions test whether two abelian varieties are exactly equal
or if one is a subset of another.

EXAMPLES

This example illustrates that taking the direct product of abelian varieties
is not commutative, in the sense of equality (though it is in the sense
of isomorphism).
\begincode
> A := JZero(11);
> B := JZero(14);
> A*B eq A*B;
true
> A*B eq B*A;
false
> IsIsomorphic(A*B, B*A);
true Homomorphism N(1) from JZero(11) x JZero(14) to JZero(14) x JZero(11) 
given on integral homology by:
[0 0 1 0]
[0 0 0 1]
[1 0 0 0]
[0 1 0 0]
\endcode

The first inclusion below is as expected, but the second
non-inclusion might be surprising.  We do not consider $J_0(11)$ as a
subset of $J_0(22)$, even though there is an injective map from one to
the other, since to be a subset is much stronger than just the
existence of an inclusion map.
\begincode
> JZero(37) subset JZero(37);
true
> JZero(11) subset JZero(22);
false
> IsInjective(NaturalMap(JZero(11),JZero(22)));
true
\endcode

 ***************************************************************************/

intrinsic 'eq'(A::ModAbVar, B::ModAbVar) -> BoolElt
{True if A and B are equal.}
   if Dimension(A) ne Dimension(B) then
      return false;
   end if;
   if not (A`field_of_definition cmpeq B`field_of_definition) then
      return false;
   end if;
   if not (A`base_ring cmpeq B`base_ring) then
      return false;
   end if;
   if IsAttachedToModularSymbols(A) then
      if not IsAttachedToModularSymbols(B) then
         return false;
      end if;
      return A`homology eq B`homology;  
   end if;
   // Now over same field and ring, and neither is attached to modular symbols.
   // so we consider them the same if they are exactly the same model for a 
   // sub-abelian variety of something attached to modular symbols.
   if Homology(A) ne Homology(B) then   // check if lattices the same lattice.
      return false;
   end if;
   eA := ModularEmbedding(A);
   eB := ModularEmbedding(B);
   return Codomain(eA) eq Codomain(eB) and Matrix(eA) eq Matrix(eB);

end intrinsic;

intrinsic 'in'(A::ModAbVar, X::List) -> BoolElt
/* EXCLUDE FROM HANDBOOK */
{True if A is equal to one of the elements of the list X.}
   return &or [A eq B : B in X];
end intrinsic; 

intrinsic 'subset'(A::ModAbVar, B::ModAbVar) -> BoolElt
{True if A is a subset of B.}
   if A eq B then
      return true;
   end if;
   if Dimension(A) eq 0 then
      return true;
   end if;
   if Dimension(A) gt Dimension(B) then
      return false;
   end if;
   eA := ModularEmbedding(A);
   eB := ModularEmbedding(B);
   if Codomain(eA) ne Codomain(eB) then
      return false;
   end if;
   if not IsInjective(eA) then
      return false;
   end if;
   if not IsInjective(eB) then
      return false;
   end if;
   return RowSpace(Matrix(eA)) subset RowSpace(Matrix(eB));
end intrinsic;


/***************************************************************************

  << Modular Embedding and Parameterization >>

Every modular abelian variety $A$ is equipped with a modular parameterization
and a modular embedding.  The modular parameterization is a surjective
homomorphism from a modular symbols abelian variety, such as $J_0(N)$.
The modular embedding is a homomorphism to a modular symbols abelian
variety, which is only guaranteed to be {\em injective in the category
of abelian varieties up to isogeny}.   The structure of these two homomorphisms
is extremely important as it is completely defines $A$.

EXAMPLES

\begincode
> X := [JZero(11),ModularAbelianVariety("37B")];
> CommonModularStructure(X);
[*
Homomorphism from JZero(11) to JZero(11) x JZero(37) given on integral 
homology by:
[1 0 0 0 0 0]
[0 1 0 0 0 0],
Homomorphism from 37B to JZero(11) x JZero(37) given on integral 
homology by:
[0 0 1 1 1 0]
[0 0 0 0 0 1]
*] [*
Homomorphism from JZero(11) x JZero(37) to JZero(11) (not printing 6x2 
matrix),
Homomorphism from JZero(11) x JZero(37) to 37B (not printing 6x2 
matrix)
*]
\endcode

The next example illustrates that the modular ``embedding'' need
only be an embedding in the category of abelian varieties up to isogeny.
\begincode
> A := JZero(37)(1);
> x := A![1/2,1];
> B := A/Subgroup([x]);
> e := ModularEmbedding(B);
> e;
Homomorphism from modular abelian variety of dimension 1 to 
JZero(37)_Qbar given on integral homology by:
[ 1 -1  1  0]
[ 2 -2 -2  2]
> IsInjective(e);
false
\endcode
Moreover, the modular parameterization is surjective, but it need
be optimal (have connected kernel).
\begincode
> pi := ModularParameterization(B);
> IsSurjective(pi);
true
> ComponentGroupOfKernel(pi);
Finitely generated subgroup of abelian variety with invariants [ 2 ]
> IsOptimal(pi);
false
\endcode




 ***************************************************************************/


intrinsic ModularParameterization(A::ModAbVar) -> MapModAbVar
{A surjective morphism to A from an abelian variety attached to
modular symbols.}
   return A`modular_param;
end intrinsic;

intrinsic ModularEmbedding(A::ModAbVar) -> MapModAbVar
{A morphism with finite kernel from A to a modular 
abelian variety attached to modular symbols.  
This is only guaranteed to be an embedding in the category 
of abelian varieties up to isogeny.}
   return A`modular_embedding;
end intrinsic;


intrinsic CommonModularStructure(X::[ModAbVar]) -> List, List
{This intrinsic finds modular abelian varieties J_e and J_p associated to modular symbols 
and returns a list of finite-kernel maps from the abelian varieties in X to J_e 
and a list of modular paramaterizations from J_p to the abelian varieties in X.}
   require #X gt 0 : "Argument 1 must be nonempty.";

   Dp := [* Domain(ModularParameterization(X[1])) *]; 
   De := [* Codomain(ModularEmbedding(X[1])) *]; 
   for i in [2..#X] do 
      B := Domain(ModularParameterization(X[i]));
      if not (B in Dp) then
         Append(~Dp, B);
      end if;
      B := Codomain(ModularEmbedding(X[i]));
      if not (B in De) then
         Append(~De, B);
      end if;
   end for;
 

   _, embed_p, param_p := DirectSum([x : x in Dp]);
   _, embed_e, param_e := DirectSum([x : x in De]);

   param_X := [* *];
   embed_X := [* *];
   for A in X do
      for i in [1..#Dp] do
         if Dp[i] eq Domain(ModularParameterization(A)) then
            Append(~param_X, param_p[i]*ModularParameterization(A));
            break;
         end if;
      end for;
      for i in [1..#De] do 
         if De[i] eq Codomain(ModularEmbedding(A)) then
            Append(~embed_X, ModularEmbedding(A)*embed_e[i]);
            break;
         end if;
      end for;
   end for;

   return embed_X, param_X;
end intrinsic;

/***************************************************************************

  << Coercion >>

Coercion can be used to create points on modular abelian varieties from
vectors on a basis of integral homology, from other elements of modular
abelian varieties, or from modular symbols.  See the examples below,
especially the last one, to understand some of the subtleties of
coercision that arise because we view an abelian variety as a
vector space modulo a lattice, and the lattice can be embedded
in any way in $\Q^n$.

EXAMPLES
If you coerce a {\em sequence} of rationals or reals into an abelian variety
$A$, then MAGMA computes the corresponding linear combination of a basis of 
integral homology and returns the point it defines.   The sequence
must have length the rank of the integral homology.
\begincode
> JZero(11)![1/2,1/5];
Element of abelian variety defined by [1/2 1/5] modulo homology
\endcode
If you coerce exactly two cusps (or extended reals) into $A$,
then MAGMA computes the point corresponding to that modular symbol.
\begincode
> JZero(11)![Cusps()|1/2,1/5];
0
> JZero(11)![Sqrt(2),0];
Element of abelian variety defined by 
[1.414213562373095048801688724198 0] modulo homology
> JZero(11)![Cusps()|0,Infinity()];     // cusps
Element of abelian variety defined by [0 1/5] modulo homology
> JZero(11)![0,Infinity()];             // extended reals
Element of abelian variety defined by [0 1/5] modulo homology
\endcode
Coercion of modular symbols also works for higher weight. 
\begincode
> JZero(11,4)![0,Infinity()];
Element of abelian variety defined by [-4/61 5/61 1/61 -1/61] 
modulo homology
> R<x,y> := PolynomialRing(RationalField(),2);
> JZero(11,4)!<x^2,[0,Infinity()]>;
Element of abelian variety defined by [-4/61 5/61 1/61 -1/61] 
modulo homology
> JZero(11,4)!<y^2,[0,Infinity()]>;
Element of abelian variety defined by [44/61 -55/61 -11/61 11/61]
modulo homology
\endcode

You can also coerce elements from abelian subvarieties into
an ambient abelian variety.
\begincode
> J := JZero(37); A := Decomposition(J)[1];
> x := A![1/5,0];
> Parent(x);
Modular abelian variety 37A of dimension 1, level 37 and 
conductor 37 over Q
> x in J;
false
> y := J!x; y;
Element of abelian variety defined by [1/5 -1/5 1/5 0] modulo 
homology
> y in J;
true
> Parent(y);
Modular abelian variety JZero(37) of dimension 2 and level 37 over Q
\endcode

Coercion also provides an easy way to create the $0$ element.
\begincode
> JZero(37)!0;
0
\endcode

The following example illustrates the subtlety of coercion
when the element being coerced in is a vector instead of a
sequence.   We create the quotient of $J_0(11)$ by a cyclic
subgroup of order $10$.  The lattice that defines $J_0(11)$
is $\Z\times \Z$, but the lattice that defines this quotient
is $(1/10)\Z\times \Z$.  Thus the natural quotient map on
vector spaces is defined by the identity matrix.  On the other
hand, the matrix
of the quotient with respect to a basis for integral homology
has determinant $10$.
\begincode
> A := JZero(11);
> x := A![1/10,0]; x;
Element of abelian variety defined by [1/10 0] modulo homology
> Order(x);
10
> B,pi := A/Subgroup([x]);
> B;
Modular abelian variety of dimension 1 and level 11 over Qbar
> pi;
Homomorphism from JZero(11)_Qbar to modular abelian variety of 
dimension 1 given on integral homology by:
[10  0]
[ 0  1]
> Matrix(pi);
[1 0]
[0 1]
> IntegralMatrix(pi);
[10  0]
[ 0  1]
> base := Basis(IntegralHomology(B)); base;
[
    (1/10    0),
    (0 1)
]
\endcode
If we coerce in the sequence [1/10,0] we get the point in $B$
that is represented by $1/10$th of the first generator for 
homology.   If we coerce in the vector $(1/10,0)$, we instead
get the element of $B$ represented by that element of the
rational homology, which is $0$, since the lattice that defines
$B$ is embedded in such a way that it contains $(1/10,0)$. 
\begincode
> y := B![1/10,0]; y;
Element of abelian variety defined by [1/10 0] modulo homology
> Order(y);
10
> z := B!base[1]; z;
0
\endcode

 ***************************************************************************/

intrinsic HackobjCoerceModAbVar(A::ModAbVar, x::.) -> BoolElt, ModAbVarElt
{Coerce x into A.  The argument x can be an element of a modular abelian
variety, the integer 0, a sequence obtained by Eltseq'ing an element of A (i.e.,
a linear combination of *integral* homology),  a vector on the basis for 
*rational* homology, or a tuple of the form <P(X,Y),[u,v]> that defines
a modular symbol.}
   require Characteristic(BaseRing(A)) eq 0 : 
         "Creation of points on abelian varieties over " * 
         "finite fields not allowed.";
    
   H := Homology(A);
   case Type(x):
      when Tup:
         t, z := IsCoercible(H,x);
         if not t then
            return false, "Could not coerce modular symbol into abelian variety.";
         end if;
         return true, A!z;
      when ModAbVarElt:
         if Parent(x) eq A then
            return true, x;
         end if;
         e := ModularEmbedding(Parent(x)); 
         f := ModularEmbedding(A);
         if Codomain(f) eq Codomain(e) then
            y := e(x);
            if Element(y) in RowSpace(Matrix(f)) then
               pi := ModularParameterization(A);
               return true, pi(y);
            end if;
         end if;   
         return false, "Illegal coercion."; 
      when SeqEnum, ModTupFldElt, ModTupRngElt, LatElt:
         lattice := Type(x) eq SeqEnum;
         if Type(x) ne SeqEnum then
            w := x;
            x := Eltseq(x);
         end if;
         if #x gt 0 and Type(x[1]) in {SetCspElt, ExtReElt} then
            t, z := IsCoercible(H, x);   // as modular symbol.
            if not t then
               return false, "Unable to coerce modular symbol into space.";
            end if;
            return true, A!z;
         end if;
       
         if #x ne Dimension(Homology(A)) then
            return false, Sprintf("Argument 2 must be a sequence of %o elements.",
                      Dimension(Homology(A)));   
         end if;
         if #x eq 0 then
            return A!0;
         end if;
         if lattice then           
            B := Basis(IntegralHomology(A));
            V := VectorSpace(FieldOfFractions(Parent(x[1])),#x);
            w := &+[V|x[i]*B[i]  : i in [1..#x]];
         end if;

         t, v := IsCoercible(VectorSpace(H), w);

         K := Qbar;
         if not t then
            t, v := IsCoercible(RealVectorSpace(H), w);
            K := CC;
         end if;
         if not t then
            return false, "Illegal coercion.";
         end if;
         z := Create_ModAbVarElt(v, A, K);
         return true, z;

      when RngIntElt, FldRatElt:
         if x ne 0 then 
            return false, "Argument 2 must be 0 if it is an integer.";
         end if;
         return true, Create_ModAbVarElt_Zero(A);
      else 
         return false, "Illegal coercion";
   end case;      
end intrinsic;




/***************************************************************************

  << Modular Symbols to Homology >>

Modular symbols determine elements of the rational homology of
$J_0(N)$, $J_1(N)$, etc., and hence of arbitrary modular abelian
varieties, by using the modular parameterization.  The commands below
convert from modular symbols, which are represented in various ways,
to vectors on the basis for rational or integral homology.  

EXAMPLES

\begincode
> A := JZero(11);
> x := ModularSymbolToIntegralHomology(A,[0,Infinity()]); x;
(  0 1/5)
> z := A!x; z;
Element of abelian variety defined by [0 1/5] modulo homology
> Order(z);
5
> A := JZero(47);
> x := ModularSymbolToIntegralHomology(A,[0,Infinity()]); x;
(-1/23  3/23 -4/23  4/23 -3/23  1/23  2/23  8/23)
> z := A!x;
> Order(z);
23
\endcode

\begincode
> J := JZero(11,4);
> IntegralHomology(J);
Lattice of rank 4 and degree 4
Basis:
( 3  1 -1 -1)
( 1 -1 -3  1)
( 2 -2  2  2)
( 2 -2  2 -2)
Basis Denominator: 8
> ModularSymbolToIntegralHomology(J,[0,Infinity()]);
(-4/61  5/61  1/61 -1/61)
\endcode
Notice that for weight greater than $2$ the homomogenous polynomial
part of the modular symbol may be omitted.  If so, it defaults to $x^{k-2}$.
\begincode
> R<x,y> := PolynomialRing(RationalField(),2);
> ModularSymbolToIntegralHomology(J,<x^2,[0,Infinity()]>);
(-4/61  5/61  1/61 -1/61)
> ModularSymbolToIntegralHomology(J,<y^2,[0,Infinity()]>);
( 44/61 -55/61 -11/61  11/61)
\endcode
The result of coercion to rational homology is different because it is
written in terms of the basis for rational homology instead of the
basis for integral homology, and in this example the two basis differ.
\begincode
> ModularSymbolToRationalHomology(J,[0,Infinity()]);
( -7/488  -9/488 -11/488  13/488)
\endcode

Coercion is also a way to define torsion points on abelian varieties. 
\begincode
> JZero(37)![1/5,0,0,0];
Element of abelian variety defined by [1/5 0 0 0] modulo homology
\endcode

***************************************************************************/

function MStoHom_Coercion(A, x)
   if not IsAttachedToModularSymbols(A) then  
              // first reduce to modular symbols case.
      pi := ModularParameterization(A);
      y := MStoHom_Coercion(Domain(pi),x);
      if Type(y) eq BoolElt then
         return y;
      end if;
      return y*Matrix(pi);
   end if;
   t,y := IsCoercible(Homology(A),x);
   if not t then
      return false;
   end if;
   return y;
end function;

intrinsic ModularSymbolToRationalHomology(A::ModAbVar, x::ModSymElt) 
     -> ModTupFldElt
{}
   ans := MStoHom_Coercion(A, x);
   require Type(ans) ne BoolElt : "Argument 2 not coercible into argument 1.";
   return ans;
end intrinsic;

intrinsic ModularSymbolToRationalHomology(A::ModAbVar, x::SeqEnum) -> ModTupFldElt
{}
   ans := MStoHom_Coercion(A, x);
   require Type(ans) ne BoolElt : "Argument 2 not coercible into argument 1.";
   return ans;
end intrinsic;

intrinsic ModularSymbolToRationalHomology(A::ModAbVar, x::Tup) -> ModTupFldElt
{The element of rational homology naturally associated to the (formal)
modular symbol s=P(X,Y)\{alpha,beta\}, where alpha, beta are in P^1(Q)
and P is a homogeneous polynomial of degree 2.
The returned vector is written with respect to the basis of 
rational homology.  This intrinsic takes its input as a pair <P(X,Y),[c,d]>, 
where c, d are in Cusps().}
   ans := MStoHom_Coercion(A, x);
   require Type(ans) ne BoolElt : "Argument 2 not coercible into argument 1.";
   return ans;
end intrinsic;

intrinsic ModularSymbolToIntegralHomology(A::ModAbVar, x::SeqEnum) -> ModTupFldElt
{}
   ans := MStoHom_Coercion(A, x);
   require Type(ans) ne BoolElt : "Argument 2 not coercible into argument 1.";
   H := Homology(A);
   C := BasisChange_Rational_to_Lattice(H);
   return ans*C;
end intrinsic;

intrinsic ModularSymbolToIntegralHomology(A::ModAbVar, x::Tup) -> ModTupFldElt
{The element of integral homology naturally associated to the (formal)
modular symbol s=P(X,Y)\{alpha,beta\}, where alpha, beta are in P^1(Q)
and P is a homogeneous polynomial of degree 2.
The returned vector is written with respect to the basis of 
integral homology.
This intrinsic takes its input a sequence [a,b] or a 
pair <P(X,Y),[a,b]>,  where a, b are in Cusps().}
   ans := MStoHom_Coercion(A, x);
   require Type(ans) ne BoolElt : "Argument 2 not coercible into argument 1.";
   H := Homology(A);
   C := BasisChange_Rational_to_Lattice(H);
   return ans*C;
end intrinsic;


/***************************************************************************

  << Embeddings >>

The {\tt Embeddings} command contains a list of embeddings (up to isogeny)
from $A$ to other abelian varieties.  The {\tt AssertEmbedding} command
allows you to add an embedding to the beginning of the list.   The embeddings
are used for computing intersections, sums, etc., with the embedding at
the front of the list having highest priority.

EXAMPLES

Every modular abelian variety comes equipped with at least one
embedding, the ``modular embedding''.
\begincode
> Embeddings(JZero(11));
[*
Homomorphism from JZero(11) to JZero(11) given on integral homology by:
[1 0]
[0 1]
*]
> A := JZero(37)(1);
> Embeddings(A);
[*
Homomorphism from 37A to JZero(37) given on integral homology by:
[ 1 -1  1  0]
[ 1 -1 -1  1]
*]
\endcode
We add another embedding to the list of embeddings for $A$.
\begincode
> phi := NaturalMap(A,JZero(37*2));
> AssertEmbedding(~A,phi);
> Embeddings(A);
[*
Homomorphism N(1) from 37A to JZero(74) given on integral homology 
by:
[-1  1 -1  0  2 -1  1  0 -2  1 -2  2 -1  2 -1  1]
[-1  0  0  0  2 -1  1 -2  0  0 -1  1 -1  2  1 -1],
Homomorphism from 37A to JZero(37) given on integral homology by:
[ 1 -1  1  0]
[ 1 -1 -1  1]
*]
\endcode
The following intersection would not make sense if we
hadn't chosen an embedding of $A$ into $J_0(74)$.
\begincode
> B := Codomain(phi)(1); B;
Modular abelian variety 74A of dimension 2, level 2*37 and 
conductor 2^2*37^2 over Q
> #(A meet B);
1
\endcode



 ***************************************************************************/

intrinsic Embeddings(A::ModAbVar) -> List
{A list of morphisms from A into abelian varieties, 
which are used in making sense of intersections, 
sums, etc.  The embeddings at the beginning of the
list take precedence over those that occur later.  Note
that these maps might not really be injective; e.g.,
the modular embedding, which need only be injective on
homology, is always at the end of this list.}
   if not assigned A`embeddings then
      A`embeddings := [* ModularEmbedding(A) *];
   end if;
   return A`embeddings;
end intrinsic;

intrinsic AssertEmbedding(~A::ModAbVar, phi::MapModAbVar)
{Place phi at the beginning of Embeddings(A).  The morphism
phi must have finite kernel.}
   require Domain(phi) eq A : "The domain of argument 2 must be A.";
   require Nullity(phi) eq 0 : "Argument 2 must have finite kernel.";
   if not assigned A`embeddings then
      dummy := Embeddings(A);
   end if;
   if #A`embeddings gt 0 and A`embeddings[1] eq phi then
      return ;
   end if;
   X := [* phi *];
   for f in A`embeddings do
      Append(~X, f);
   end for;
   A`embeddings := X;
end intrinsic;

function AllCharactersTrivial(A)
   return &and[IsTrivial(DirichletCharacter(M)) : 
              M in ModularSymbols(Homology(Codomain(ModularEmbedding(A))))];
end function;

function Can_Compute_Inner_Twists(A)
   return IsAttachedToModularSymbols(A) and #ModularSymbols(A) eq 1 and 
     IsNew(ModularSymbols(A)[1]) and #NewformDecomposition(ModularSymbols(A)[1]) eq 1 ;
end function;


/***************************************************************************

  << Base Change >>

The {\tt BaseExtend} and {\tt ChangeRing} commands allow you to change
the base ring of a modular abelian variety.  The {\tt BaseExtend}
command is the same {\tt ChangeRing}, but is more restrictive.  For
example, if $A$ is over $\Q$ then {\tt BaseExtend(A,GF(2))} is not
allowed, but {\tt ChangeRing(A,GF(2))} is.  

Abelian varieties can have their base ring set to a finite field, but
there is very little that is implemented for abelian varieties over
finite fields.  Computing the number of points is implemented, but
creation of actual points or homomorphisms is not.

EXAMPLES
We consider $J_(13)$ over many fields and rings. 
\begincode
> A := JOne(13);
> BaseExtend(A,CyclotomicField(7));
Modular abelian variety JOne(13) of dimension 2 and level 13 over 
Q(zeta_7)
> BaseExtend(A,AlgebraicClosure(RationalField()));
Modular abelian variety JOne(13)_Qbar of dimension 2 and level 13 
over Qbar
> BaseExtend(A,RealField());
Modular abelian variety JOne(13)_R of dimension 2 and level 13 over R
> BaseExtend(A,ComplexField());
Modular abelian variety JOne(13)_C of dimension 2 and level 13 over C
> ChangeRing(A,GF(3));
Modular abelian variety JOne(13)_GF(3) of dimension 2 and level 13 
over GF(3)
> #ChangeRing(A,GF(3));
19 19
> B := ChangeRing(A,GF(13)); B;
Modular abelian variety JOne(13)_GF(13) of dimension 2 and level 13
over GF(13)
> IsAbelianVariety(B);
false
> ChangeRing(A,Integers());
Modular abelian variety JOne(13)_Z of dimension 2 and level 13 over
Z
> ChangeRing(A,PolynomialRing(RationalField(),10));
Modular abelian variety JOne(13) of dimension 2 and level 13 over 
Polynomial ring of rank 10 over Rational Field
Lexicographical Order
Variables: $.1, $.2, $.3, $.4, $.5, $.6, $.7, $.8, $.9, $.10
\endcode

 ***************************************************************************/

intrinsic BaseExtend(A::ModAbVar, R::Rng) -> ModAbVar
{Extend the base ring of A to R, if possible.  This is merely
a more restrictive version of ChangeRing.}
/* todo: probably need better tests here. */
   S := BaseRing(A);
   require not (IsField(S) and (Characteristic(S) ne Characteristic(R) or 
          Type(R) in {RngInt, RngOrd})) :
      "The base ring of argument 1 and argument 2 are incompatible."*
      " Try using ChangeRing instead.";
   return ChangeRing(A,R);
end intrinsic;

intrinsic ChangeRing(A::ModAbVar, R::Rng) -> ModAbVar
{Change the base ring of A to R, if possible.}
   if R cmpeq BaseRing(A) then
      return A;
   end if;
   if Type(R) eq FldRat and assigned A`associated_abvar_over_Q then
      return A`associated_abvar_over_Q;
   end if;
   B := Copy_ModAbVar(A);
   B`base_ring := R;
   delete B`endomorphism_ring;
   delete B`endomorphism_algebra;
   delete B`hecke_algebra;
   delete B`isogeny_decomp;
   delete B`is_simple;
   delete B`is_selfdual;
   delete B`conductor;
   delete B`inner_twists;
   delete B`cm_twists;
   delete B`decomposition;
   delete B`tamagawa_numbers;
   delete B`embeddings;  // would get infinite loops if try to base extend these...
   if not assigned B`associated_abvar_over_Q then
      if Type(BaseRing(A)) eq FldRat then
          B`associated_abvar_over_Q := A;
      end if;
   end if;
   if Characteristic(R) ne 0 then
      if Type(FieldOfDefinition(A)) eq FldRat then
         B`field_of_definition := GF(Characteristic(R));
      end if;
      B`field_of_definition := R;
   end if;
   Rname := RingName(R);
   if B`name ne "" and Rname ne "" and #Rname lt 8 then
      B`name := Sprintf("%o_%o",B`name,Rname);   
   end if;
   // base extend modular parameterizations, etc.
   if IsAttachedToModularSymbols(A) then
      B`modular_param := IdentityMap(B);
      B`modular_embedding := IdentityMap(B);
   else
      J := ChangeRing(Domain(ModularParameterization(A)),R);
      B`modular_param`domain := J;
      B`modular_param`codomain := B;
      B`modular_param`parent := Hom(J,B);
      B`modular_embedding`domain := B;
      B`modular_embedding`codomain := B`modular_param`domain;
      B`modular_embedding`parent := Hom(B,J);
   end if;
   // It is tempting to base extend all of them, instead of throwing 
   // them away.  However, one could easily get in an infinite loop 
   // by doing so.
   B`embeddings := [* B`modular_embedding *];
   return B;
end intrinsic;

intrinsic CanChangeRing(A::ModAbVar, R::Rng) -> BoolElt, ModAbVar
{True if it is possible to change the base ring of A to R, and A over R 
when possible.}
   t := true;
   /* todo: actually add some conditions! */
   return t, ChangeRing(A,R);
end intrinsic;



function HasOnlyTrivialCharacters(A)
   J := Codomain(ModularEmbedding(A));
   for M in ModularSymbols(J) do
      if not IsTrivial(DirichletCharacter(M)) then
         return false;
      end if; 
   end for;
   return true;
end function;



DEFAULT_COLUMNS := 70;

function Format(s, start, only_first)
/* Reformat string s so each new line starts with start, and
   so the lines have length GetColumns(). */
   cols := GetColumns()-4;
   if GetVerbose("ModAbVar") in {1,3} and #s gt cols * 4 then
      s := &*[s[i] : i in [1..cols*4]] * "...(truncated)";
   end if;
   if GetVerbose("ModAbVar") in {2,4} and #s gt cols * 10 then
      s := &*[s[i] : i in [1..cols*10]] * "...(truncated)";
   end if;
   if cols lt 1 then
      cols := DEFAULT_COLUMNS;
   end if;
   ans := "";
   IDX := [i : i in [1..#s] | s[i] ne "\n"];
   if #IDX eq 0 then
      return "";
   end if;
   s := &*[s[i] : i in IDX];
   while #s gt 0 do 
      ans := ans * start;
      if only_first then
         start := &*[" " : i in [1..#start]];
      end if;
      /*if #s lt cols - #start then
         ans := ans * s;
         return ans;
      end if;*/
      j := cols-#start;
      while j ge 1 and j le #s and not (s[j] in {" ", "\t", "\n"}) do 
         j := j - 1;
      end while;
      if j eq 0 then
         j := cols-#start;
      end if;
      if j ge 1 then
         add := &*[s[i] : i in [1..j] | i le #s];
      else   
         add := "";
      end if;
      ans := ans * add;
      if #start + #add lt cols then
         ans := ans * &*[" " : i in [#start+#add+1..cols]] ;
      end if;
      //ans := ans * " |";
      if j+1 le #s then
         s := &*[s[i] : i in [j+1..#s]];
      else
         s := "";
      end if;
      if #s gt 0 then
         ans := ans * "\n";
      end if;
   end while;
   return ans;
end function;

function VerboseHeader(cmd, level)
   LEN:=45;
   if #cmd gt LEN then
      cmd := &*[cmd[i] : i in [1..LEN]] * "...";
   else
      cmd := cmd * &*[" " : i in [#cmd+1..LEN+3]];
   end if;
   return Sprintf("*  ModAbVar-Verbose %o: %o (time %o)\n", level, cmd, Cputime());
end function;


function VerboseFooter()
   cols := GetColumns();
   if cols eq 0 then
      cols := DEFAULT_COLUMNS;
   end if;
   return &*["-" : i in [1..cols-2]];
end function;

procedure Verbose(cmd, msg1, msg2)
   file := false;
   level := GetVerbose("ModAbVar");
   if level eq 0 then
      return;
   end if;
   msg1 := Sprintf("%o",msg1);
   if level eq 3 then
      level := 1;
      file := true;
   elif level ge 4 then
      level := 2;
      file := true;
   end if;
   if level gt 1 then
      msg2 := Sprintf("%o",msg2);
      msg1 := msg1 * msg2;
   end if;
   out := VerboseHeader(cmd, 
           level eq 1 or #msg2 eq 0 select "1" else "2");
   out := out* Format(msg1, "|  ", false);
   out := out* VerboseFooter();
   if file then
      Write("ModAbVar-verbose.log",out);
   else
      print out;
   end if;
end procedure;
