CoCalc
Sharedwww / Tables / modsym.gpOpen in CoCalc
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\
\\ The program modsym.gp 
\\
\\
\\ HISTORY: 
\\
\\    *) Original program written for PARI 1.x by Joseph L. Wetherell
\\
\\    *) July 1999: William Stein updated the program for PARI 2.x.
\\
\\    *) March 2001: Dominique Bernardi and Bernadette Perrin-Riou made 
\\                   the modifications described below:
\\
\\       Les modifications sont de type : dans le programme modsym.gp, 
\\       les relations binaires n'�taient pas toutes prises en compte,
\\       par exemple, certaines relations du type 2*x = 0 (voir par 
\\       exemple le cas N=6, N=106).   Nous avons donc d� enlever 
\\       l'astuce qui permettait de r�duire  les matrices de taille 
\\       N � une taille N/2, ce qui ralentit bien s�r le programme.
\\       Nous avons d'autre part rajout� un signe +1 ou -1 afin 
\\       de calculer aussi l'application de l'espace des formes 
\\       modulaires dans la partie + ou - de l'homologie.
\\       Enfin, nous avons rajout� deux proc�dures ellsym et 
\\       xpm qui permettent de calculer les composantes x^\pm E 
\\       du symbole modulaire associ�e � une courbe elliptique 
\\       donn�e dans GP/Pari par ellinit ou plut�t � la forme 
\\       modulaire asssoci�e.  Le but est de calculer des valeurs 
\\       de fonctions L p-adiques associ�es � une courbe elliptique 
\\       et � un discriminant de corps quadratiques (positif ou n�gatif).
\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

talk=1;    \\ verbosity flag.

{mnormal(cond,v)=
   lift(v*Mod(1,cond));
}

{generatemsymbols(cond,
         num,ret,curn)=
\\ num = [\Gamma:\Gamma_0(N)] = N * Prod_{p|cond} (1+p^-1)
   num = cond;
   v = factor(cond);
   for(n=1,length(v~),
     num = num*(1+1/v[n,1]);
   );

\\ initialize M-symbol list
   ret = vector(num,c,0);
   curn = 1;

\\ generate M-symbols in three lists:
\\ list 1: (c:1) for 0 <= c < N
   for(c=0, cond-1,
     ret[curn] = [c,1];
     curn = curn + 1;
   );

\\ list 2: (1:d) for 0 <= d < N and gcd(d,N)>1
   for(d=0, cond-1,
     if(gcd(d,cond)>1,
       ret[curn] = [1,d];
       curn = curn + 1;
     ,);
   );

\\ list 3: (c:d) with c|N, c<>1,N,  gcd(c,d)=1, gcd(d,N)>1
   v = curn;
   for(d=2, cond-1,
     if(gcd(d,cond)>1,
         fordiv(cond,c,
           if(c<>1 && c<>cond && gcd(c,d)==1,
             if(0 == sum(n=v,curn-1, ((ret[n][1]*d-c*ret[n][2])%cond)==0 ),
               ret[curn] = [c,d];
               curn = curn + 1;
             ,);
           ,);
         );
     ,);
   );

if (curn -num - 1, print("pas le bon nombre de symboles"));
   ret;
}

{hashmsymbols(cond,msymbols,
         ret,tmp) =
   ret = matrix(cond,cond,c1,d1,0);
   for(n=1,cond-1,
     if(gcd(n,cond)==1,
       for(k=1,length(msymbols),
         tmp = [msymbols[k][1]*n,msymbols[k][2]*n];
         tmp = mnormal(cond,tmp);
         ret[tmp[1]+1,tmp[2]+1] = k;
       );
     ,);
   );
   ret;
}

{dohash(cond,hash,v,
         tmp)=
   v=mnormal(cond,v);
   hash[v[1]+1,v[2]+1];
}

doV(v) = [-v[1],v[2]]
doS(v) = [-v[2],v[1]]
doT1(v) = [-v[1]-v[2],v[1]]
doT2(v) = [v[2],-v[2]-v[1]]

{msymbolrelations(cond,msymbols,hash,signe,
                 mask1,n1,n2,n3,s1,s2,s3,relat,nsym,freebasis,symtofree,
         tempr) =
   nsym  = length(msymbols);
   relat = matid(nsym);
   mask1 = vector(nsym,n,1);


   for(n=1,nsym,
     n1 = msymbols[n];
     n2 = doS(msymbols[n]);
     n1 = dohash(cond,hash,n1);
     n2 = dohash(cond,hash,n2);
     r = relat[n1,] + relat[n2,];
     minc = vecmax(abs(r));
     for(k=1,nsym,
       if(r[k] && abs(minc)>=abs(r[k]),
         minc = r[k];
         mincloc = k;
       ,);
     );
     if(minc,
       mask1[mincloc] = 0;
       r = r/minc;
       relat = relat - relat[,mincloc]*r;
     ,);
   );

   for(n=1,nsym,
     n1 = msymbols[n];
     n2 = doT1(msymbols[n]);
     n3 = doT2(msymbols[n]);
     n1 = dohash(cond,hash,n1);
     n2 = dohash(cond,hash,n2);
     n3 = dohash(cond,hash,n3);
     r = relat[n1,] + relat[n2,] + relat[n3,];
     minc = vecmax(abs(r));
     for(k=1,nsym,
       if(r[k] && abs(minc)>=abs(r[k]),
         minc = r[k];
         mincloc = k;
       ,);
     );
     if(minc,
       mask1[mincloc] = 0;
       r = r/minc;
       relat = relat - relat[,mincloc]*r;
     ,);
   );

        for(n=1,nsym,
     n1 = msymbols[n];
     n2 = doV(msymbols[n]);
     n1 = dohash(cond,hash,n1);
     n2 = dohash(cond,hash,n2);
     r = relat[n1,] -signe*relat[n2,];
     minc = vecmax(abs(r));
     for(k=1,nsym,
       if(r[k] && abs(minc)>=abs(r[k]),
         minc = r[k];
         mincloc = k;
       ,);
     );
     if(minc,
       mask1[mincloc] = 0;
       r = r/minc;
       relat = relat - relat[,mincloc]*r;
     ,);
   );

   n1 = vector(nsym,n,mask1[n]*n);
   n2 = eval(Set(simplify(n1)));
   freebasis = vector(length(n2)-1,n,n2[n+1]);
   symtofree = matrix(nsym,length(freebasis),n1,n2,relat[n1,freebasis[n2]]);
   [freebasis,symtofree~];
}

cusps;
{initcusps() = cusps=[];}

{testcuspeq(cond,cusp1,cusp2,
         p1,q1,p2,q2,s1,s2,n)=
   p1 = cusp1[1];
   q1 = cusp1[2];
   p2 = cusp2[1];
   q2 = cusp2[2];
   s1 = if(q1>2, lift(Mod(1,q1)/Mod(p1,q1)), 1);
   s2 = if(q2>2, lift(Mod(1,q2)/Mod(p2,q2)), 1);

   foo=Mod(s1*q2-s2*q1, gcd(q1*q2,cond))==0;

   foo;
}

{getcusp(cond,v,signe,
         ret,negv) =
   ret = 0;
   negv = [-v[1],v[2]];
   for(n=1,length(cusps),
     if(testcuspeq(cond,v,cusps[n][1]), ret=[n,cusps[n][2]],
     if(testcuspeq(cond,negv,cusps[n][1]), ret=[n,signe*cusps[n][2]]));
   );
   if(ret==0,
     cusps = concat(cusps, [[v, 1-((signe < 0) && testcuspeq(cond,v,negv))]]);
     ret = length(cusps);
     ret = [ret,cusps[ret][2]];
   ,);
   ret;
}

{delta(cond,msymbol,
         v,c1,c2,ret)=
   v = bezout(msymbol[1],msymbol[2]);
   if(v[3]<>1, print("msymbol not coprime:", msymbol);1/0 ,);
   c1 = [v[2],msymbol[1]];
   c2 = [-v[1],msymbol[2]];
   ret=[c2/gcd(c2[1],c2[2]), c1/gcd(c1[1],c1[2])];
   ret;
}

{kerdelta(cond,msymbols,hash,freebasis,signe,
         rels,trels,m,k,kinv,ksz)=
   initcusps();
   rels = vector(length(freebasis),n, delta(cond,msymbols[freebasis[n]]) );
   for(n=1,length(rels),
     rels[n] = [getcusp(cond,rels[n][1],signe),getcusp(cond,rels[n][2],signe)];
   );
   m = matrix(length(cusps),length(freebasis),n1,n2,0);
   for(n=1,length(rels),
     m[rels[n][1][1],n] = m[rels[n][1][1],n] + rels[n][1][2];
     m[rels[n][2][1],n] = m[rels[n][2][1],n] - rels[n][2][2];
   );
   k = matker(m);
   ksz = matsize(k);
   kinvr = (matsupplement(k))^(-1);
   [k, vecextract(kinvr,2^ksz[2]-1,2^ksz[1]-1)];
}

{modulartofree(cond,hash,symtofree,v,
         q,hashval,n1,s1,sym,mods,ret)=
   if(v[2]==0, v=[1,cond]);
   q = contfrac(v[1]/v[2]);
   q[1]=1;
   for(n=3,length(q),
     q[n]=q[n]*q[n-1]+q[n-2];
   );
   sym = vector(matsize(symtofree)[2],n,0)~;
   for(n=2,length(q),
     mods = [(-1)^(n)*q[n],q[n-1],1];
     n1 = dohash(cond,hash,mods);
     sym[n1] = sym[n1]+1;
   );
   ret=symtofree*sym;
   ret;
}

{tpmatrix(cond,msymbols,hash,relations,hdata,p,
         freebasis,symtofree,hbase,hdual,
         sz,ret,col,sym,cusp,a)=
   if(!isprime(p) || cond%p == 0,
      print("Error, this program can only compute T_p for p prime and
coprime to N.");
      return([]);
   );
   freebasis = relations[1];
   symtofree = relations[2];
   hbase = hdata[1];
   hdual = hdata[2];

   sz = matsize(hbase);
   ret = matrix(sz[2],sz[2],n1,n2,0);

   for(c=1,sz[2],
     col = vector(sz[1],n,0)~;
     for(r=1,sz[1],

       a = hbase[r,c];
       if(a,
         sym = msymbols[freebasis[r]];
                 cusp = delta(cond,sym);
         col = col -
           a*modulartofree(cond,hash,symtofree,cusp[1]*[p,0;0,1]);
         col = col +
           a*modulartofree(cond,hash,symtofree,cusp[2]*[p,0;0,1]);
         for(n=0,p-1,
           col = col -
              a*modulartofree(cond,hash,symtofree,cusp[1]*[1,0;n,p]);
           col = col +
              a*modulartofree(cond,hash,symtofree,cusp[2]*[1,0;n,p]);
         );
       ,);
     );
     ret[,c] = hdual*col;
   );

   ret;
}

{gammap(cond,msymbols,hash,relations,hdata,p,
         freebasis,symtofree,hbase,hdual,
         sz,ret,col,sym,cusp,a)=
   freebasis = relations[1];
   symtofree = relations[2];
   hbase = hdata[1];
   hdual = hdata[2];

   sz = matsize(hbase);

   col = vector(sz[1],n,0)~;
   for(n=1,p-1,
     col = col + modulartofree(cond,hash,symtofree,[n,p]);
   );

   hdual*col;
}

defaultmodsymspace=0;
{modsym(N,signe,
        symbols, hash, relations, hdata,
        t1,t2,t3,t4)=
    gettime();
    if(talk,print1("1. Generating M-symbols \t\t("));
    symbols   = generatemsymbols(N);
    if(talk,print(t1=gettime()," ms)");
            print1("2. Hashing M-symbols \t\t\t("));
    gettime();
    hash      = hashmsymbols(N,symbols);
    if(talk,print(t2=gettime()," ms)");
            print1("3. Quotienting out by relations  \t("));
    gettime();
    relations = msymbolrelations(N,symbols,hash,signe);
    if(talk,print(t3=gettime()," ms)");
            print1("4. Computing the kernel of delta \t("));
    gettime();
    hdata     = kerdelta(N,symbols,hash,relations[1],signe);
    if(talk,print(t4=gettime()," ms)"));
    if(talk,print("Total time.......................\t(",
            t1+t2+t3+t4," ms)"));

    defaultmodsymspace=[N,symbols,hash,relations,hdata,signe];
}

{T(p, M=0)=
    if(type(M)=="t_INT",M=defaultmodsymspace;
            if(M==0,print("T: Error -- no default space.")));
    if(matsize(M)[2] != 6,
      print("T: invalid modular symbols object."));
    tpmatrix(M[1],M[2],M[3],M[4],M[5],p,M[6]);
}

{ellsym (ell,signe,
   modele, cond, M, r, temp, res)=
   modele = ellglobalred(ell);
   cond = modele [1];
   M = modsym (cond, signe);
   r = matsize (M[5][1]) [2]; print("r = ", r);
   res = if (r > 1, temp = [;];
     forprime (p = 2, 10000,
       if (matrank (temp) == r-1, break,);
       print ("p = ", p);
       if (cond % p, temp = concat (temp, T (p, M) - ellap(ell, p)*matid(r)),)
     );
     mattranspose (matker (mattranspose (temp))) * M[5][2],
     M[5][2]
   );
   defaultellsyms = [cond, M[3], M[4][2], res, signe, ell]
}

{xpm (a, b, parms)=
   if (parms,, parms = defaultellsyms);
   (parms [4] * modulartofree (parms [1], parms [2], parms [3], 
     [a,b], parms [5]))[1];
}


{f(N,p)=
	modsym(N,1); print(factor(charpoly(T(p),x)));
	modsym(N,-1); print(factor(charpoly(T(p),x)));
}