Sharedwww / Tables / modsymold.gpOpen in CoCalc
\\ Weight 2 modular symbols in GP/PARI pre 2.0
\\  (by Joseph L. Wetherell)


\\ the M-symbol n(c:d) is represented by [c,d,n]

\\ We identify the symbols (c:d)=(-c:d)=(c:-d)=(-c:-d)
\\                           = -(d:c) = -(d:-c) = -(-d:c) = -(-d:-c)
\\ This gives all of the 2-term and neg-eigenspace identities,
\\ and reduces the size of the spaces.
{mnormal(cond,v, 
	tmp)=
  v = lift(v*mod(1,cond));
  if(v[1]>cond/2, v[1]=cond-v[1] ,);
  if(v[2]>cond/2, v[2]=cond-v[2] ,);
  if(v[1]>v[2],   tmp=v[1]; v[1]=v[2]; v[2]=tmp; v[3]=-v[3] ,);
  v;
}

{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,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,1];
      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(0, n=v,curn-1, ((ret[n][1]*d-c*ret[n][2])%cond)==0 ),
              ret[curn] = [c,d,1];
              curn = curn + 1;
            ,);
	  ,);
        );
    ,);
  );

\\ normalize all of the m-symbols
  for(n=1,num, 
    ret[n] = mnormal(cond,ret[n]);
    ret[n] = ret[n][2]*cond+ret[n][1];
  );
\\ remove duplicates
  ret = set(ret);
  
  ret = vector(length(ret), n, [ret[n]%cond, ret[n]\cond, 1] );
  for(n=1,length(ret),
    while(gcd(ret[n][1],ret[n][2])<>1,
       ret[n] = ret[n] + [cond,0,0];
    );
  );
  ret;
}


{hashmsymbols(cond,msymbols, 
	ret,tmp) = 
  ret = matrix(cond\2+1,cond\2+1,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,msymbols[k][3]];
        tmp = mnormal(cond,tmp);
        ret[tmp[1]+1,tmp[2]+1] = tmp[3]*k;
      );
    ,);
  );
  ret;
}

{dohash(cond,hash,v,  
	tmp)=  
  v=mnormal(cond,v);
  tmp = hash[v[1]+1,v[2]+1];
  [abs(tmp), sign(tmp)*v[3]];
}
doT1(v) = [v[1]+v[2],v[1],v[3]]
doT2(v) = [v[2],v[1]+v[2],v[3]]

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

  for(n=1,nsym,
    n1 = msymbols[n];
    n2 = doT1(msymbols[n]);
    n3 = doT2(msymbols[n]);
    n1 = dohash(cond,hash,n1);  s1=n1[2]; n1=n1[1];
    n2 = dohash(cond,hash,n2);  s2=n2[2]; n2=n2[1];
    n3 = dohash(cond,hash,n3);  s3=n3[2]; n3=n3[1];
    r = s1*relat[n1,] + s2*relat[n2,] + s3*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;
    ,);
  );

  n1 = vector(nsym,n,mask1[n]*n);
  n2 = 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/p1,q1)), 1);
  s2 = if(q2>2, lift(mod(1/p2,q2)), 1);

  mod(s1*q2-s2*q1, gcd(q1*q2,cond))==0;
}

{getcusp(cond,v,
	ret,negv) = 
  ret = 0;
  negv = [-v[1],v[2]];
  for(n=1,length(cusps),
    if(testcuspeq(cond,v,cusps[n]), ret=n, );
    if(testcuspeq(cond,negv,cusps[n]), ret=n, );
  );
  if(ret==0,
    cusps = concat(cusps,[v]);
    ret = length(cusps);
  ,);
  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,
	rels,m,k,kinv,ksz)=
  rels = vector(length(freebasis),n, delta(cond,msymbols[freebasis[n]]) );
  for(n=1,length(rels),
    rels[n] = [getcusp(cond,rels[n][1]),getcusp(cond,rels[n][2])];
  );
  m = matrix(length(cusps),length(freebasis),n1,n2,0);
  for(n=1,length(rels),
    m[rels[n][1],n] = m[rels[n][1],n] + 1;
    m[rels[n][2],n] = m[rels[n][2],n] - 1;
  );
  k = ker(m);
  ksz = matsize(k);
  kinvr = matinvr(supplement(k));
  [k, matextract(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 = cf(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];
    hashval = dohash(cond,hash,mods); 
	s1=hashval[2]; n1=hashval[1];
    sym[n1] = sym[n1]+s1;
  );
  ret=symtofree*sym;
  ret;
}

{tpmatrix(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);
  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;
}