Sharedwww / Tables / model.gpOpen in CoCalc
\\ model.gp -- find models for modular curves.
\\ 
\\ As usual, it is assumed that a matrix E of new eigenforms exists
\\ on the PARI stack before this program is run. 
\\ The format of E is
\\     E[N,i] = [f(x), [a1(x), a2(x), ..., an(x)]]. 
\\ where each ai is a polynomial in Q[x] of deg < d = deg f 
\\ and f is an irreducible polynomial in Z[x].  The ai 
\\ should be thought of as expressions in a root of f,
\\ and we assume that they are integral over Z. 
\\
\\ bg.gp is assumed "available".
\\
\\ William A. Stein ([email protected]), Nov 12, 1998.

{nrows(M)=matsize(M)[1];}
{ncols(M)=matsize(M)[2];}
{find(v,a,i=0)=for(i=1,ncols(v),if(v[i]==a,return(i)));return(0);}

\\ intbasis1: -- int basis for space spanned by a single conjugacy class.
\\ INPUT  : N,i which defines an E[N,i] = [f,[a1,...,an]]
\\ OUTPUT : vector of q-expansions.
{intbasisclass(N,i, 
  M,K,v,j,k)=
  K=nfinit(E[N,i][1]);
  M=matrix(poldegree(E[N,i][1]),ncols(E[N,i][2]),j,k,0);
  for(j=1,ncols(E[N,i][2]),M[,j]=nfalgtobasis(K,E[N,i][2][j]));
  vector(nrows(M),j,sum(k=1,ncols(M),M[j,k]*q^k)+O(q^(ncols(M)+1)));
}

{intbasisnew(N,
  i,j,k,n,last,v,w)=
  print1("Integral basis of newforms at level ", N);
  n=0;last=0;
\\  if(nrows(E)<N,print("Need data for level ",N," to proceed!"); return(););
  for(i=1,ncols(E),
    if(E[N,i]!=0,
      n+=poldegree(E[N,i][1]),last=i-1;break())
  );
\\  if(n==0,print("Need data for level ",N," to proceed!"); return(););
  v=vector(n,i,0);
  j=1;
  for(i=1,last,
    w=intbasisclass(N,i);
    for(k=1,ncols(w), v[j]=w[k]; j++)
  );
  print(" -- ",ncols(v)," newforms.");
  v;
}


\\ intbasis:
\\ Find a Q-basis for S_2(Gamma_0(N)) consisting of forms with q-expansion
\\ coefficients in Z. 

{intbasis(N,
  fac,i,j,k,l,m,v,d,t,Snew)=
  fac=factor(N);
  v=vector(nrows(fac)*N*20,i,0);   \\ way over estimate.
  j=1;
  fordiv(N,d,
    Snew=intbasisnew(N/d);
    fordiv(d,t,
      for(k=1,ncols(Snew),v[j]=subst(Snew[k],q,q^t);j++);
    );
  );
  print("Finished computing integral basis at level ",N," dim = ", j-1);
  vector(j-1,i,v[i]);               \\ truncate away extra zeros.
}


\\ reducedintbasis:
\\ Compute a "reduced" integral basis for S_2(Gamma_0(N)),
\\ that is a Z-basis for the free abelian group of modular forms
\\ whose q-expansion coefficients lie in Z.

\\ save for later -- not sure how to do this in pari.

\\    if(sum(i=1,s,e[i])==d,  \\ monomial of degree d,

\\ monomials
\\ returns list of exponents of monomials of degree d in s variables.
{monomials(s,d,
  v,w,i,j,k,n) =
  if(s<1||d<1,return([]));
  if(s+d-1==d,v=vector(1,i,0),       \\ bug in PARI binomial function.
    v = vector(binomial(s+d-1,d),i,vector(s,n,d)));
  k=1;
  v[k]=vector(s,n,if(n==1,d,0)); k++;
  for(i=0,d,
    w = monomials(s-1,d-i);
    for(j=1,ncols(w),v[k]=vector(s,n,if(n==1,i,w[j][n-1]));k++)
  );
  v;
}

\\ findrels:
\\ This function finds all linear degree d relations satisfied by a 
\\ a list {f1, ..., fs} of power series.
\\
\\ INPUT:  * A list f1,...,fs of q-expansions a1*q+a2*q^2+...+O(q^(prec+1)).
\\         * An integer n.
\\         * prec = precision of the fi's. 
\\
\\ OUTPUT: * A matrix of integers.  The rows of the matrix represent
\\           linear relations satisfied by f1^d, f1^{d-1}f2, ..., fs^d.
\\           The rows are linear independent. 
\\
\\ AUTHOR:  William Stein ([email protected]), Nov 1998. 

{findrelmat(v, d, prec, 
   M, i, j, k, mon, g) = 
  \\ A: construct M.
  \\    iterate through the set of monomials in s variables of degree d. 
  print("Computing monomials of degree ", d,".");
  mon = monomials(ncols(v),d);
\\  prec=poldegree(truncate(v[1]));
\\  for(i=2,ncols(v),prec=min(prec,poldegree(truncate(v[i]))));
  print("Building ", prec,"x", ncols(mon), 
        " matrix of monials in modular forms.");
  M = matrix(prec,ncols(mon),i,k,0);
  j=1;
  for(i=1,ncols(mon),
    g = 1;
    for(k=1,ncols(v), 
      g*=(v[k]+O(q^(prec+1)))^mon[i][k]);         \\ g <--  "f1^{e1}*f2^e2*..." 
    for(k=1,prec,
      M[k,j]=polcoeff(g,k-1));
    j++;
  );
  M;
}


{findrels(v, d, prec, 
    M)=
  M = findrelmat(v,d,prec);
  print("Computing LLL reduced relations using matkerint.");
  R=matkerint(M);    \\ LLL reduced kernel.
  print("Done computing relations.");
  R;
}



\\ findindeprels:
\\ Find all _independent_ relations of degree <= d. 
\\ 
\\ Method: Compute relations of degree 2.  These automatically give
\\         rise to several relations of degree 3.  Compute the 




{itoa(n)=if(n==1,return("a"));if(n==2,return("b"));if(n==3,return("c"));
if(n==4,return("d"));if(n==5,return("e"));if(n==6,return("f"));
if(n==7,return("g"));if(n==8,return("h"));if(n==9,return("i"));if(n==10,return("j"));
if(n==11,return("k"));if(n==12,return("l"));if(n==13,return("m"));if(n==14,return("n"));
if(n==15,return("o"));if(n==16,return("p"));if(n==17,return("q"));if(n==18,return("r"));
if(n==19,return("s"));if(n==20,return("t"));if(n==21,return("u"));if(n==22,return("v"));
if(n==23,return("w"));if(n==24,return("x"));if(n==25,return("y"));if(n==26,return("z"));
}  

{printbasis(v, prec=10, 
  i)=
  for(i=1,ncols(v),
    if(i<=26,print1(itoa(i)),print1("f",i,));
    print(" = ",v[i]+O(q^prec)));
}

show_basis(N)=printbasis(intbasis(N),10);

{printmon(mon,
  i,first)=
  first=1;
  for(i=1,ncols(mon),
    if(mon[i]!=0,
      if(!first,print1("*"));
      first=0;
      if(i<=26,print1(itoa(i)),
        print1("f",i));
      if(mon[i]>1,
        print1("^",mon[i])
      )
    );
  );
}

{printrel(r, s, d,  
  i, j, k, mon, rel,first) = 
  mon = monomials(s,d);
  if(nrows(r)==0, print("no relations of degree ", d); return());
  for(k=1,ncols(r), 
    first=1;
    for(i=1,ncols(mon),
      if(r[i,k]!=0,
        if(!first && r[i,k]>0,print1(" + "));
        first=0;
        if(r[i,k]<0, print1(" - "); r[i,k]=-r[i,k]);
        if(abs(r[i,k])!=1,
          print1(r[i,k],"*"));
        printmon(mon[i]);
      );
    );
    print(" = 0");
  );

}


{show_rels(N,d,prec,
  v,r)=
  v=intbasis(N);
  r=findrels(v,d,prec);
  printrel(r, ncols(v), d);
}

{initleveldata()=
E=matrix(400,40,x,y,0);
print("reading levels to 100");
read("1-100.gp");
print("reading levels to 200");
read("101-200.gp");
print("reading levels to 300");
read("201-300.gp");
\\print("reading levels to 400");
\\read("301-400.gp");
0;
}


\\ Let r be a list of relations of degree d.
\\ This function returns the list of relations of degree d+1 induced by the 
\\ relations r.  For example, if r "is"
\\        3*a^2+b^2
\\ Then this function would return the 2 relations
\\      3*a^3+a*b^2 and 3*a^2*b+b^3

{inducedrels(r,s,d,
  ind,i,j,k,l,mon_d,mon_dp1,mon,xmon,v)=
  mon_d = monomials(s,d);      \\ monomials of degree d
  mon_dp1 = monomials(s,d+1);  \\ monomials of degree d+1.
  \\ v=vector(ncols(mon_d),i,0);
  ind=matrix(ncols(mon_dp1),s*ncols(r), i,j,0); 
  j=1;
  for(i=1,s,
    for(k=1,ncols(r),
      for(l=1,nrows(r),
        if(r[l,k]!=0,
          mon=mon_d[l];
          mon[i]++;                \\ multiply by ith variable
          xmon=find(mon_dp1,mon);  \\ find in the degree d+1 monomial list
          if(xmon==0,print("inducedrels: ERROR IN PROGRAMMING!"));
          ind[xmon,j]=r[l,k];
        );
      );
      j++;
    );
  );
  ind;
}

{numcubicrels(N,
  g)=
  g=g0(N);
  g*(g+7)*(g-4)/6 + 5;
}

{show_inducedrels(N,d,prec,
  v,r,r2)=
  v=intbasis(N);
  r=findrels(v,d,prec);
  print("THE DEGREE ", d, " RELATIONS:");
  printrel(r, ncols(v), d);
  print("THE INDUCED DEGREE ", d+1, " RELATIONS:");
  r2=inducedrels(r,ncols(v),d);
\\  printrel(r2, ncols(v), d+1);
  print("computing rank of an ",matsize(r2)," matrix");
  print("THERE RANK OF INDUCED RELS IS ", matrank(r2));
  print("THE NUMBER OF CUBIC RELS IS   ", numcubicrels(N));
}


print("initleveldata() = load data for levels <= 300");
print("show_basis(N)   = show Q-basis at level N");
print("show_rels(N,d,prec) = compute the relations of degree d at level N, ");
print("                      to precision prec<=500.");
print("show_inducedrels(N,d,prec) = ");

{systematic(N,
   m,prec)
  while(N<=220,
    prec=floor(mu0(N)/3.0)+15;
    show_inducedrels(N,2,prec);
    N+=2;
  );
}