Documentation for the modular symbol programs in the eclib library (guide to source code)

%  created 20070125
%  Time-stamp: <2012-05-07 17:30:02 jec>

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Note that until March/April 2012 my source code was split into four
subdirectories (procs, qcurves, qrank and g0n) only the last of which
contained modular symbol code for computing elliptic curves.  The
notes here are mostly about the latter.  In the current distribution,
there are library headers and implementation files in libsrc/eclib and
libsrc/ respectively, test programs in tests/ and user programs in progs/.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Level 0:  functionality provided from other libraries

Basic arithmetic, including polynomials and PARI/NTL interfacing (PARI
is only used for integer factorization), linear algebra (including
sparse LA), various utilities some of which are independently useful
(reduction of cubics and quartics, solving conics, ...).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Level 1:  low-level classes probably not useful for users to use
directly.

NB The efficiency of the implementations here is of crucial
importance!

----------------------------------------------------------------------

file: moddata
class: level

Holds the current level N ("modulus"), its prime factors and divisors,
some "global" flags and a reduce() function for reduction modulo N.

----------------------------------------------------------------------

file: moddata
class: moddata (derived from level)

Holds precomputed arithmetic data for the level, including lookup
tables of gcds with N and inverse modulo N

Test program: modtest (also tests symb class)

----------------------------------------------------------------------

file: symb
class: symb

Class for (c:d) or M-symbols.  Also holds a pointer to a moddata.
Functions for equality testing, normalization.

----------------------------------------------------------------------

file: symb
class: modsym

Class for {a,b} or modular symbols, where a,b are rational.
Includes constructor from a symb.

----------------------------------------------------------------------

file: symb
class: symblist

Class for an array of symbs with hash function for fast index lookup.

----------------------------------------------------------------------

file: symb
class: symbdata (derived from moddata)

Class for full M-symbol data.  Includes a symblist of "specials"
(M-symbols not of the form (c:1) or (1:d)),  bijections from the set
of M-symbols to/from indices in the range 0..nsymb-1, and
implementation of matrix operations (rof, rsof, sof, tof) applied to
symb indices.

----------------------------------------------------------------------

file: oldforms
class: oldforms

Class for retrieving oldforms from stored data files for all levels
properly dividing the current level.  Computes multiplicities of
these. Used for one purpose: in searching for newforms recursively, at
any point we will have restricted to some subspace cut out by
eigenvalues [a2, a3, a5, ...]  and the oldforms member function
dimoldpart() will tell us how much of the current dimension is
accounted for by oldforms.  This enables us to abandon the current
recursive branch if it is all old, and for the same reason results in
oldforms being discarded and only newforms kept.  [NB for this to work
it is crucial that before running the newforms search, date for the
dividing levels is available.  This is handled automatically by the
newforms class functions.]

----------------------------------------------------------------------

file: cusp
class: cusplist

Class to hold an array of inequivalent cusps (=rationals), insert and
locate cusps in the list, and test for cusp equivalence (using a
moddata pointer to know the level and also whether or not plusflag is
set, in which case the equivalence test is affected).

----------------------------------------------------------------------

file: homspace
class: mat22

Class to hold a 2x2 matrix, with functions to allow them to act on
rationals and symbs

----------------------------------------------------------------------

file: homspace
class: matop

Class to hold a formal sum of mat22s.  Constructors include:  Hecke
operators T_p, W_q, Heilbronn matrices

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Level 3: high-level classes useful for users to use directly.  I would
expect these to be wrapped for Sage (eventually)

----------------------------------------------------------------------

file: homspace
class: homspace (derived from symbdata)

The main class for holding a "modular symbols space".
The constructor does a lot of work:
  homspace(long n, // the level
	   int hp, // plus-space flag (0 or 1 or -1)
	   int hcusp, // cuspidal flag (0 or 1)
	   int verbose // verbosity (0 : no output
	               //	     1 : basic short output
	               //            2 : lots of detail)
	   );

The plus-space flag switches between working on H_1^+ (by quotienting
out extra relations (c:d)=(-c:d)), H_1^-,  or the full H_1.

If the cupsidal flag is 1 then the space computed is the kernel of
the boundary map to the cusps, and all operators are restricted to
that subspace.  However the "restricted" operators which compute
directly the action of Tp etc on a subspace only work properly (1) for
subspaces of the dual and (2) for the full (including non-cuspidal)
homology.

Functionality (most useful class functions only):

  long h1cuspdim() // returns dim of cuspidal subspace
  long h1dim()     // returns dim of whole space
  long h1ncusps()  // returns number of cusp classes

  svec schain(const symb& s) // returns coordinates of symb s
  svec schaincd(long c, long d) // returns coords of symb (c:d)
  svec schain(long nn, long dd) // returns coords of modsym {0,nn/dd}
  svec schain(const rational& r) // returns coords of modsym {0,r}

  vec cycle(long n, long d) // as schain(n,d) but cuspidal only,
  result is a coordinate vector w.r.t. the basis for cuspidal homology

  vec cycle(const rational& r) // as schain(r) but cuspidal
  vec cycle(const modsym& m) // as schain(b)-schain(a) but cuspidal, m={a,b}

  mat heckeop(long p, int dual, int display=0) // p'th Hecke operator
  (Tp or W_p according as p does not / does divide level)
  NB (1) if cuspidal==1, restricts to cuspidal subspace
     (2) if dual==1, transposes

  mat heckeop_restricted(long p, const subspace& s, int dual, int
  display=0)
  // as previous but restricted to *dual* subspace s, so dual==1 is
  // mandatory else results are not guaranteed to be meaningful

  smat s_heckeop(long p, int dual, int display=0)
  smat s_heckeop_restricted(long p, const ssubspace& s, int dual, int
  display=0)
  // Exactly as previous two but returns a sparse matrix

  mat newheckeop(long p, int dual, int display=0) // Tp computation
  using Heilbronn matrices.  Faster than heckeop(), but not currently
  implemented in a "restricted" version, hence not used in the main
  newform-finding programs. *But* for a one-off calculation of a Hecke
  op on the whole space you should use this one!

  mat wop(long q, int dual, int display=0) // returns W_q matrix
  smat s_wop(long q, int dual, int display=0) // returns W_q sparse matrix
  mat fricke(int dual, int display=0) // returns Fricke involution matrix
  mat conj(int dual,int display=0) // returns conjugation matrix
  (=identity if plusflag==1!)
  mat conj_restricted(const subspace& s, int dual,int display=0) // ditto, restricted to subspace s
  smat s_conj(int dual, int display=0) // as above, returning sparse mat
  smat s_conj_restricted(const ssubspace& s, int dual, int display=0)  // as above, returning sparse mat

  vec maninvector(long p) // for good p, returns sum_{a mod p}{0,a/p}
  vec manintwist(long p) // for good p, returns sum_{a mod p}chi(a){0,a/p}

----------------------------------------------------------------------

file: newforms
classes: newform and newforms

Overview:  newform holds an individual newform, with a pointer to the
parent newforms for "global" data;  newforms (which is derived from
level) contains a vast anout of stuff, explained below, including an
array of individual newform-s.

----------------------------------------------------------------------

class: newform

For more details see newforms.h

Main data fields:

  newforms *nf;  // the "parent"
  int plusflag;   // 1 for old-style newform, 0 for old-style h1newform
  vec bplus,bminus; // DUAL eigenvectors (bminus only used if plusfalg==0)
  scalar type;            // 2 for rectangular, 1 for triangular
			  //  period lattice
  vector<long> aplist, aqlist;
// aplist is a vector of Fourier coefficients, indexed by all primes
// aqlist is a vector of the Wq eigenvalues
// Bother are needed since (a) the q'th Fourier coefficient is 0 when
// q^2|N so does not determine the eigenvalue; (b) the range of the
// aplist may not go ar enough to include all bad primes, but we really
// want to have all aq (to get the sign of the functional equation,
// analytic rank parity etc.)

  long ap0;     // Eigenvalue of first "good" p
  long sfe;     // sign of functional equation
  rational loverp;  // L(f,1)/x where x = least real part of a period
                    // =np0/dp0

// The bext few are technical, needed for computing periods.
  long cuspidalfactor;  // pdot =cuspidalfactor*np0
  long pdot,np0,dp0,qdot;  // np0=1+p0-ap0, pdot = maninvector(p0).bplus,
                           //                    = cuspidalfactor*dp0
                           // qdot={q,infinity}.bplus
  long lplus, lminus;  // primes = +1, -1 mod 4
  long mplus, mminus;  // mplus*x=sqrt(lplus)*L(fxlplus,1)
                       // mminus*yi=sqrt(-lminus)*L(fxlminus,1)
  long a,b,c,d,dotplus,dotminus; // matrix for period integration
  // Either type=1, lattice=[2x,x+yi]
  // Or     type=2, lattice=[x,yi]
  // & integral over [a,b;Nc,d] is dotplus*x+dotminus*yi

  long degphi;             // degree of Weil parametrization
  // Not currently being set or used, since our algorithm was too slow
  // for large levels and we just use MW's programs instead (externally).

  vec coords;              // vector components of each freegen (will
			   // become one column of homspace::projcoord)
  // Used in mass-production of ap after newforms have been found

Member function for newform class:

Constructors:
  newform(const vector<int>& data, const vector<long>& aq, const vector<long>& ap, newforms* nfs);
// To be used when constructing from stored data, so no basis vector
//  -- since we do not compute the H1 space unless we have to
  newform(const vec& vplus, const vec& vminus, const vector<long>& ap, newforms* nfs,long ind=-1);
// To be used when constructing from the recursive search, including
  basis vector(s) (vminus only relevant when plusflag==0);  all the
  technical data fields will be computed in this constructor, using
  the parent newforms's homspace.
// New constructor(s) needed, from an elliptic curve and/or the data
  [level, aplist, aqlist] which will need to find the basis vector(s)
  from splitting off a 1D eigenspace from the appropriate homspace and
  then use the first constructor above.

  void add_more_ap(int nap); // Extends the aplist (if necessary)
  until it contains atleast nap terms.  Only used in the second
  constructor after all the newforms for a given level have been
  found, so that they all have apslists of the same length.

Format of stored data:  newforms are stored in files in a directory
whose name is set via the environment variable NF_DIR which defaults
to "newforms".  Filenames within that dir are the level (in decimal)
preceded by 'x', e.g. newforms/x11 holds the data for level 11. The
filename is constructed by the fiunction nf_filename(), in moddata.h/cc.

The filename is constructed by the function
char* nf_filename(long n, char c)  // use c='x'
where it is the callers responsibility to delete[] the result after
use.

The data files are written by newforms::output_to_file().

The file format is binary (for speed of input/output and compact size
-- my definitive newforms dirsctory is 8GB for levels to 130000).  The
following description is taken from newforms.h:

/* Data stored in a newform (and in data files newforms/x$N):
   (Numbers refer to lines of data file)
Items 1-18 are "int" while the ap and aq are "short"
3.  sfe : sign of functional equation (=-product of aq)
4.  ap0 : p0'th Hecke eigenvalue, p0=smallest good prime
5.  np0 : np0=1+p0-ap0
6.  dp0 : dp0/np0=L/P=L(f,1)/2x
7.  lplus : prime =1 (mod 4) with L(f,lplus,1) nonzero
8.  mplus : L(f,lplus,1)*sqrt(l)=mplus*x
9.  lminus : prime =3 (mod 4) with L(f,lminus,1) nonzero
10. mminus : L(f,lminus,1)*sqrt(-l)=mminus*yi
11-14. a, b, c, d : entries of a matrix M=[a,b;N*c,d] in Gamma_0(N) s.t.
15. dotplus       : the integral of f over {0,M(0)} is
16. dotminus      : dotplus*x+dotminus*yi
17. type : type 1 if period lattice = [2x,x+yi], type 2 if [x,yi]
18. degphi : degree of modular parametrization
aq : list of Wq-eigenvalues at bad primes
ap : list of Tp- & Wq-eigenvalues at all primes
*/

There are some shell scripts for viewing the stored data:

nnf N : returns number of newforms at level N (on file -- no computation!)
nap N : returns number of ap stored at level N (on file -- no computation!)
showdata N : shows just the technical data, not all the ap
showeigs N : shows all the aq and ap
shownf N : shows all data for level N (pipe through more or less!)

NB At present these scripts assume NF_DIR="newforms"

----------------------------------------------------------------------

class: newforms (derived from level and splitter_base)

For more details see newforms.h

Parent classes: level provides... the level, bad primes.
       		splitter_base provides facilities for (a) recursive
       		searching for Hecke eigenspaces from scratch using
       		only the ability to compute a sequence of operators
       		(optionally, restricted to the current space), with a
       		known finite set of possible eigenvalues for each; and
       		(b) splitting of the eigenspace for any given (valid)
       		sequence of eigenvalues.

		The reason this is done in a virtual base class way
		over in ../procs is that I use the *exact* same
		functionality of the code for modular symbols over
		imaginary quadratic fields.

Constructor:   newforms(long n, int plus, int cuspidalflag, int disp)
// just sets the flags, and the of and h1 pointers to 0
// the real work is done by calling either createfromscratch() or
// createfromdata(), see below.

Main data fields:

  int verbose; // controls output verbosity level (*none* if 0)
  long maxdepth, // bound on recursion depth in search
       cuspidal, // flags whether or not to use cuspidal homology
       plusflag; // flags whether or not to work in + space
  int basisflag;  // flag to determine how much work ::use() does.
  // use() is what the base class does with eigenspaces when found.

  // flag==0 for a recursive search, where use() will call the newform
  // constructor;

  // flag==1 for recovering newforms in a plusspace from existing data,
  // where all use() has to do is set the newforms's basis vector, as the
  // newform will already have been constructed with all the other data
  // needed.

  vec mvp; // the "manin vector" sum_{a mod p} {0,a/p} for p=p0, the
      	   //  smallest good prime
  map<long,vec> mvlplusvecs, mvlminusvecs;
  	   // arrays of quadratic-twisted manin vectors, indexed by
      	   // good primes l where l=1(4) for plusvecs and l=3(4) for
      	   // minusvecs.  Each newform has one of each (determined by
      	   // its own lplus and lminus primes);  but we store these
      	   // vectors in the newforms class for efficiency when more than
      	   // one newform shares the same lplus or lminus!

  oldforms* of; // pointer to oldforms whcih provde the dimoldpart()
  	    	// facility to enable old spaces to be discarded
		// during the recursive search

  homspace* h1; // pointer to homspace for the level.  We make this a
  	    	// member rather than deriving the newforms class from
  	    	// homspace, which would be more natural, as for some
  	    	// functionality (e.g. reading the newforms data from
  	    	// files and recomputing the elliptic curves) we do 
  	    	// not need the homology information, so don't want to
  	    	// compute it.  Basically, createfromscratch() make
  	    	// the homspace (by calling ::makeh1()) and
  	    	// createfromdata() does not.

  int j0; 
  long nq, dq; // data used for ap computation
  std::set<long> jlist; 
  long n1ds, // the number of newforms
       j1ds; // index used when recreating newforms
  vector<newform> nflist; // the actual newform-s.


  // These just call the similar function in the associated homspace 
  // Here i is the index of the prime to use
  mat opmat(int i, int d, int v=0) 
  mat opmat_restricted(int i, const subspace& s, int d, int v=0) 
  smat s_opmat(int i, int d, int v=0) 
  smat s_opmat_restricted(int i, const ssubspace& s, int d, int v=0) 

  // The next 3 functions are to provide functionality needed by the
  // base splitter class

  // the dimension of the underlying homspace
  long matdim(void)  {return h1->dimension;} 
  // The list of possible eigenvalues for the i'th operator
  vector<long> eigrange(int i)
  // Given an initial sequence of eigenvalues [a2,a3,...], uses the
  // oldforms member to return the corresponding dimension of oldforms
  long dimoldpart(const vector<long> l);

  void display(void)  // Output to stdout
  void output_to_file(int binflag=1) // output to data file in NF_DIR

  // add newform with basis b1, eiglist l to current list (b2 not used):
  void use(const vec& b1, const vec& b2, const vector<long> l); 

  // find newforms using homology; ntp is number of eigenvalues to use
  // for oldforms, *not* the number computed via homology (use addap()
  // for that):
  void createfromscratch(long ntp);

  // read newforms from file, if it exists, otherwise (perhaps) revert
  // to createfromscratch:
  void createfromdata(long ntp, int create_from_scratch_if_absent=1);

  // Construct bases (homology eigenvectors) from eigenvalue lists:
  void makebases();

  // computes vector of ap, one for for each newform	
  vector<long> apvec(long p);  

  // compute ap for all primes up to the last'th, insert them into each
  // newforms own aplist:
  void addap(long last); // adds ap for primes up to the last'th prime

  // Sort newforms 
  void sort(int oldorder=0);
  
  // next three functions() implemented in periods.cc

  // Given newform with no intdata, compute least real (part of)
  // period -- unless sfe=-1 and n=square, in which case return 0
  int get_real_period(long i, bigfloat& x, int verbose=0) const;

  // Given all data, compute the periods as a Cperiods
  Cperiods getperiods(long i, int method=-1, int verbose=0);

  // Given all data & Cperiods, compute the curve
  Curve getcurve(long i, int method, bigfloat& rperiod, int verbose=0);

  // next two implemented in pcprocs.cc

  // Computes x0, y0 (real & imag parts of periods) & a matrix which
  // gives these scaled by dotplus & dotminus.  rp_known is set if we
  // know x0 to be the least real part of a period (usually true)
  int find_matrix(long i, long dmax, int&rp_known, bigfloat&x0, bigfloat&y0);

  // Given an imaginary period y1, finds a prime lminus =3(mod 4) and
  // <=lmax for which L(f,lminus,1) is nonzero and hence a multiple
  // mminus of y1.
  // if lmax==0 it carries on until a suitable lminus is found
  int find_lminus(long i, long lmax, const bigfloat& y1);

----------------------------------------------------------------------

file: periods
class: character

Precomputes values of a quadratic character modulo a prime (or trivial)

----------------------------------------------------------------------

file: periods
class summer 

Base class for several others all of which compute sums of the form
sum_{n=1}^{\infty} a_n f(n) 
for various real or complex-values f(n), where the a_n are the
coefficients of a newform.  Using the usual Buhler-Gross-Zagier
recursion.

Data fields: see periods.h

----------------------------------------------------------------------

file: periods
class periods_via_lfchi (derinved from summer)

Class to compute the periods of a newform using the values of
L(f,chi,1) for suitable quadratic characters chi.

Constructor:
  periods_via_lfchi (const level* iN, const newform* f); 

Member functions:
  void compute(void); // does the work
  bigfloat rper(void) // return the real period
  bigfloat iper(void) // return the imaginary period
  Cperiods getperiods() // return the period lattice

----------------------------------------------------------------------

file: periods
class: periods_direct (derived from summer)

Class to compute the periods of a newform by integrating over {0,g(0)}
for a matrix g in Gamma_0(N);  from the real and imaginary parts of
the result and the integer scaling factors stored in the newform class
(and the lattice type) we can recover the full period lattice.

Constructor:
  periods_direct (const level* iN, const newform* f); 

Member functions:
  void compute(void); // Computes integral over {0,g(0)} where g is
       		      // the matrix in Gamma_0(N) stored in the newform
  void compute(long a, long b, long c, long d);  
// Computes integral over {0,g(0)} where g is [a,b;N*c,d]
  bigfloat rper(void) // return the real period
  bigfloat iper(void) // return the imaginary period
  Cperiods getperiods() // return the period lattice

----------------------------------------------------------------------

file: periods
class: part_period  (derived from summer)

Class to compute the integral of a newform over {z,\infty} for any z
in the upper half plane -- for example, to compute Heegner points.

Constructor:
  part_period (const level* iN, const newform* f); 

Member functions:
  void compute(const bigcomplex& z0); // do the work
  bigcomplex getperiod() 	      // get the answer

----------------------------------------------------------------------

file: periods
class: ldash1   (derived from summer)

Class to compute the analytic rank r and value L^(r)(f,1) for a
newform f.  The parity of r comes from the newform's s.f.e., and
whether or not r=0 is also determined from the newform's data.  This
gives an initial guess of r\in{0,1,2}.  If the value is very small
then r is increased by 2 and the value recomputed.  (Tested ok for
lots of rank 3 curves and the rank 4 curve at N=234446.)

Constructors:
  ldash1 (const level* iN, const newform* f); 
  ldash1 (const newforms* nf, long i);  // the i'th newform

Member functions:
  void compute(void); // does the work
  long rank() 	      // returns the rank
  bigfloat value()    // returns the value

----------------------------------------------------------------------

file: periods
class: lfchi (derived from summer)

Class to compute L(f,chi,1) for a newform f and quadratic character
chi (modulo an odd prime l not dividing the level)

Constructor:
  lfchi (const level* iN, const newform* f);

(Does nothing specific to the prime l)

Member functions:
  void compute(long ell); // does the work for the prime l
  bigfloat value(void) 	  // gives the value
  bigfloat scaled_value(void) // gives the value * sqrt(l)


----------------------------------------------------------------------

file: pcprocs

Contains two new functions plus implementations of newforms class functions
newforms::find_matrix() and newforms::find_lminus() described above

// returns a rational approximation a/b to real x
void ratapprox(bigfloat x, long& a, long& b);

int get_curve(long n, long fac, long maxnx, long maxny,
	      const bigfloat& x0, const bigfloat& y0, 
	      long& nx, long& ny, int& type, int detail=0);

Input 

n: target conductor of E
fac: known factor of c6(E)
x0, y0: reals known to be multiples of least real/imag periods of E
maxnx, maxny: bounds on these multiples
detail: verbose flag

Output: nx, ny, type

If type==1 then [2*x,x+y*i] or if type==2 then [x,y*i], where x=x0/nx,
y=y0/ny, is the period lattice of an elliptic curve of conductor n
(whose c6 is divisible by fac).


----------------------------------------------------------------------

files: nfd.h, nfd.cc, tnfd.cc
class: nfd

Class which implements d-dimensional newforms (i.e. with Hecke
eigenvalues of degree d over Q).

Program tnfd is a test program for this.  Experimental, not currently
working.


----------------------------------------------------------------------

files: fixc6.h, fixc6.cc, fixc4.data, fixc6.data
class: fixc6

These can be completely ignored, since I gave up using them after
N=130000 and just use multiprecision floating point always, even
though it is slower, with some cleverness to automatically increase
precision as needed.  Future plans include re-introducing the standard
double precision code (which is all still there) and using that where
possible and sensible, for speed and flexibility.

// Background: for curves with large c6 (or c4), the value cannot be computed
// with sufficient precision without using multiprecision floating
// point arithmetic, which slows down all the other cases
// unnecessarily.  To avoid this, after first computing the "correct"
// c6 value once in multiprecision, we add that value to a table
// indexed by (level, form#), and look up in the table when we need
// the value.

// The fixc6 class has two static data members, of type map<
// pair<long,int>, bigint> such that an entry ((N,i),c6) or ((N,i),c4)
// says that the c6 or c4 value for form i at level N is c6 or c4.  Of
// course, for most (N,i) pairs this is blank -- and we must avoid
// inserting wrong dummy entries of the form ((N,i),0).

// April 2005: added facility for fixing c4 as well as c6, but the
// class name is unchanged

NB For levels up to 130000 these files are quite big:
   443   1329  10319 fixc4.data
 25808  77424 679655 fixc6.data
 26251  78753 689974 total

That is, out of 570226 (optimal) curves, 25808 require this "fixing"
of c6 and 443 also for c4.

NB This is *not used* when the programs are comiled with
multiprecision floating point (i.e. NTL not NTL_INTS) so can be
ignored for Sage.  But the programs which produce tables of curves
from stored data *do* recompute the periods of the newforms and hence
get the curves from their periods, and this is *much* faster using
standard C doubles!