\\ \\ a few GP/PARI scripts for unimodular lattices (GP/PARI version 2.13) \\ \\ by Bill Allombert and Gaetan Chenevier, 13/11/2020 \\ \\ most of the gp scripts below are already documented in the files: \\ \\ -- http://gaetan.chenevier.perso.math.cnrs.fr/unimodular_lattices/tools_unimodular_lattices.gp \\ \\ -- http://gaetan.chenevier.perso.math.cnrs.fr/unimodular_lattices/tools_unimodular_lattices_no_roots.gp \\ \\ However, several functions have been optimized below (or at least, improved). \\ \\ See also the following file for a short tutorial, whose aim is to study L28: \\ http://gaetan.chenevier.perso.math.cnrs.fr/unimodular_lattices/tuto_L28.gp \\ \\ \\ functions to install \\ install("F2m_rank","lG"); install("zm_to_ZM","G"); install("ZM_to_zm","G"); install("ZM_to_F2m","G"); install("F2m_mul","GG"); install("F2m_transpose","G"); install("F2m_to_Flm","G"); install("zm_mul","GG"); install("Flm_mul","GGU"); install("Flm_det","GU"); install("FpM_det","GG"); install(hash_GEN,uG); \\ \\ d-neighbors of I_n \\ \\ nei(v), with v=[d, iso, t, ***] : returns a gram matrix of the d-neighbor \\ of the standard lattice I_n associated with : \\ -- the Z/d-isotropic line generated by the column matrix iso , \\ -- the parameter t , necessarily 0 for d odd, and in {0,1} for d even. \\ nei(v)=I_nei(v[2],v[1],v[3]); I_nei(v,d,t=0)={ my(n, a, N, nei); n=#v~; if( d%2==1, a=(v~*v)[1,1]*(d+1)/(2*d), a=(v~*v)[1,1]/(2*d)+t*d/2 ); N=matconcat( [matid(n), v/d] ); N[1,1]=d; N[1,n+1]=N[1,n+1]-a; for(i=2,n,N[1,i]=-v[i,1]); nei=1/d*mathnf(d*N); return(nei~*nei); }; \\ \\ Gram matrices \\ \\ lllgram(gram) : returns an lll-reduced form of the Gram matrix gram . \\ lllgram(gram)={ my(basis); basis=qflllgram(gram); return(basis~*gram*basis); }; concatgram(g)={ my(res=matid(0),m,n,p,q); for(u=1,#g, [m,n]=matsize(res); [p,q]=matsize(g[u]); res=matrix(m+p,n+q,i,j, if( (j<=n) && (i<=m), res[i,j],0)+if( (j>n) && (i>m), g[u][i-m,j-n],0) ) ); return(res); }; \\ is_exceptional(gram) returns 1 iff the rank n unimodular lattice L with Gram \\ matrix gram (assumed to have no norm 1 vector) has a characteristic vector \\ with norm <8 \\ (our convention is that for n = 0 mod 8, this means that L is even) is_exceptional(gram,verb=0)= { my(n=#gram, w, N); w=gram^-1*matrix(n,1,i,j,gram[i,i]%2); if( n%8==0, return( w ==matrix(n,1,i,j,0) )); N=mathnf( matconcat([w, 2*matid(n)]) ); if(verb,return(qfminim(N~*gram*N,n%8)[1])); return( qfminim(N~*gram*N,n%8)[1] <> 0); }; \\ \\ mass formulae \\ \\ even_unimodular_mass_formula(n) : returns the mass of the rank n even \\ unimodular lattices (following Serre) \\ odd_unimodular_mass_formula(n) : returns the mass of the rank n odd \\ unimodular lattices (following Conway-Sloane) \\ tot_mass(n) : returns their sum (and 1 for n=0 as a convention) \\ no1_mass(n) : returns the mass of the rank n unimodular lattices with no \\ vector of norm 1 (and 1 for n=0 as a convention) \\ we have sum no1_mass(n) x^n * exp(x/2) = sum tot_mass(n) x^n \\ even_unimodular_mass_formula(n)= { if(n%8==0, return( abs( bernfrac(n/2)/n*prod( j=1,n/2-1,bernfrac(2*j)/(4*j) ) ) ), return(0) ); }; odd_unimodular_mass_formula(n)= { my(k,E); if(n==1, return(1/2) ); if(n%8==1 || n%8==7, k=(n-1)/2; return((2^k+1)/2*abs(prod(j=1,k,bernfrac(2*j)/(4*j)))) ); if(n%8==3 || n%8==5, k=(n-1)/2; return((2^k-1)/2*abs(prod(j=1,k,bernfrac(2*j)/(4*j)))) ); if(n%8==4, k=n/2; return((1-2^(-k))*(1-2^(1-k))*1/2*1/k!*abs(bernfrac(k)*prod(j=1,k-1,bernfrac(2*j)))) ); if(n%8==0, k=n/2; return((1-2^(-k))*(1+2^(1-k))/(2*k!)*abs(bernfrac(k)*prod(j=1,k-1,bernfrac(2*j)))) ); if(n%8==2 || n%8==6, k=n/2; E=vector((k+1)/2,u,1); for(u=1,(k-1)/2,E[u+1]=-sum(v=0,u-1,binomial(2*u,2*v)*E[v+1])); E=E[#E]; return(1/(2^(2*k+1)*(k-1)!)*abs(E*prod(j=1,k-1,bernfrac(2*j)))) ); }; tot_mass(n)= { if(n==0, return(1), return(odd_unimodular_mass_formula(n) +if(n%8==0,even_unimodular_mass_formula(n),0)) ); }; no1_mass(n)= { return(sum(i=0,n, (-1)^i*tot_mass(n-i)*2^(-i)*1/i! )); }; \\ \\ root systems \\ \\ root_system(gram) : returns the isomorphism class of the root system \\ (as a string) of the Gram matrix gram . \\ name_to_surname(rs) : transforms an isomorphism class of root system of type \\ string in a similar object of type vector. \\ surname_to_name(rs) : inverse transformations. \\ root_system(gram,opt=0)= { my(un, unstr, ba, rho, heights, ind_basis, basis, rank, vec_nei, bi); if(opt==0, un=qfminim(gram,1)[3]; if(un, unstr=Str("I",#un), unstr=""); if(#un==#gram,return(unstr)); if(#un, ba=matkerint(un~*gram); gram=ba~*gram*ba ), un=0; unstr="" ); \\ Rplus = the sytem of positive roots chosen by Pari my(Rplus=qfminim(gram,2)[3]); if(#Rplus==0, return(unstr)); if(un, unstr=concat(unstr," ")); \\ rho = Weyl vector of Rplus (half the sum of positive roots) rho=1/2*Rplus*matrix(#Rplus,1,i,j,1); \\ heights = column matrix of heights of the positive roots heights=Rplus~*(gram*rho); \\ ind_basis = indices of the simple roots ind_basis=select(x->heights[x,1]==1, vector(#Rplus,u,u)); rank=#ind_basis; \\ new indexation of the simple roots, from 1 to rank basis=vector(rank,i,Rplus[,ind_basis[i]]); \\ vec_nei[i] is the list of neighbors of the i-th simple root \\ in the Dynkin diagram of the root system vec_nei=vector(rank,u,List()); for(i=1,rank, bi=basis[i]~*gram; for(j=i+1,rank, if(bi*basis[j],listput(vec_nei[i],j); listput(vec_nei[j],i)) ) ); return(Str(unstr,print_name(vec_to_mset(tree_comp(vec_nei))))); }; \\ \\ (naive) scripts for connected components of a disjoint union of trees \\ tree_comp(T)= { if(#T==0, return( List() ) ); my(S, list_comp, new_comp); S=[1..#T]; list_comp=List(); while(#S, new_comp=connected_comp(S[1],T); listput(list_comp,concat(typ_comp(T,new_comp),#new_comp)); S=setminus(S,new_comp) ); return(list_comp); }; connected_comp(i,T,j=0)= { my(vois,comp); vois=setminus(Set(T[i]),Set(j)); comp=Set(i); if(#vois==0, return(comp)); for(k=1,#vois, comp=setunion(comp,connected_comp(vois[k],T,i)) ); return(comp); }; typ_comp(T,comp)= { my(v,s,typ="A"); for(i=1,#comp, v=T[comp[i]]; if(#v==3, s=sum(j=1,3,#T[v[j]]); if(s>4,typ="E",typ="D") ) ); return(typ); }; \\ sorting functions vec_to_mset(v)= { my(res); res=vecsort(v,cmp); res=vector(#res,i,[res[i],1]); res=gather(res); return(res); }; gather(v)= { my(x,n,i,res); if( #v==0, return(v) ); i=1; x=v[1][1]; n=v[1][2]; res=vector(#v); for(k=2, #v, if( v[k][1]==x, n=n+v[k][2], res[i]=[x,n]; i=i+1; x=v[k][1]; n=v[k][2] ) ); res[i]=[x,n]; return( vector(i,u,res[u]) ); }; \\ \\ naming root systems \\ \\ surname = string, such as "2A1 A2 D4 E8") \\ name = vector of [ mult, letter, rank] (like [ [2, "A", 1], [1, "A", 2], [1, "D", 4], [1, "E", 8] ]) \\ name_to_surname(rs)= { my(res=[], m); m=0; res=vector( sum(i=1, #rs, rs[i][1]) ); for(i=1, #rs, for(j=1,rs[i][1], m++; res[m]=concat(rs[i][2], rs[i][3]) ) ); return(print_name(vec_to_mset(res))); }; surname_to_name(str)= { my(res=[],mult,letter,number); str=List(Vec(str)); while(#str, mult=""; while(str[1]<>"I" && str[1]<>"A" && str[1]<>"D" && str[1]<>"E", mult=concat(mult,str[1]); listpop(str,1) ); if(mult=="",mult=1, mult=eval(mult)); letter=str[1]; listpop(str,1); number=""; while(#str && str[1]<>" ", number=concat(number,str[1]); listpop(str,1) ); if(#str && str[1]==" ",listpop(str,1)); number=eval(number); res=concat(res,[[mult,letter,number]]) ); return(res); }; print_name(v)= { my(name,space); name=""; for(i=1,#v, if(i==1, space="", space=" " ); name=concat(name, if(v[i][2]==1, Str(space,v[i][1]), Str(space,v[i][2],v[i][1]) ) ) ); return(name); }; \\ \\ imax(R) : maximum n such that A_n-1 embeds in root system R \\ imax(surname)= { my(name, vec, typ, rk); if(#surname==0,vec=[1], name=surname_to_name(surname); vec=vector(#name); for(i=1,#name, typ=name[i][2]; rk=name[i][3]; if(typ=="A",vec[i]=rk+1, typ=="D",vec[i]=rk, typ=="E" && rk==8, vec[i]=9, typ=="E" && rk==7, vec[i]=8, typ=="E" && rk==6, vec[i]=6 ) ) ); return(vecmax(vec)); }; \\ \\ order of Weyl group \\ \\ order_weyl(rs) = order of the Weyl group of the root system rs \\ order_pmweyl(rs,rk) = order of the subgroup of O(R^rk) generated by \\ the Weyl group of rs (a root system in R^rk) and -id. \\ order_weyl_irr(letter,n)= { if(letter=="A", return((n+1)!), letter=="D", return(2^(n-1)*n!), letter=="E", if( n==8, return(696729600), n==7, return(2903040), n==6, return(51840) ) ); }; order_weyl(rs)= { if(type(rs)=="t_STR", rs=surname_to_name(rs)); return( prod(i=1,#rs,order_weyl_irr(rs[i][2],rs[i][3])^(rs[i][1]))); }; test_two(typ,n)= { if(typ=="A" && n==1, return(1), typ=="A" && n>1, return(0), typ=="D" && (n%2)==1, return(0), typ=="D" && (n%2)==0, return(1), typ=="E" && n==8, return(1), typ=="E" && n==7, return(1), typ=="E" && n==6, return(0) ); }; is_minus1_weyl(rs,rk)= { my(rk_rs,is2=1); if(type(rs)=="t_STR", rs=surname_to_name(rs)); rk_rs=sum(i=1,#rs,rs[i][1]*rs[i][3]); if(rk_rs>rk,return("error")); if(rk_rs && rk_rs<>rk, is2=0, for(i=1,#rs,is2=is2*test_two(rs[i][2],rs[i][3])) ); if(rk_rs==0,is2=0); return(is2); }; order_pmweyl(rs,rk)= { rs=surname_to_name(rs); return(order_weyl(rs)*(2-is_minus1_weyl(rs,rk))); }; \\ \\ order of isometry group \\ reduce(q,B,t=20000)= { my(V=qfminim(q,B),n=V[1]/2,M=V[3],d=#q); if(n==0,return(0)); if(F2m_rank(ZM_to_F2m(M))!=d,return(0)); my(W=qfminim(q,B-1)[3],m=#W,kstart); if (m==0,kstart=n,kstart=1); for(k=kstart, #q, for(i=1, t, my(N=Mat(vector(d,i,if(i>k,W[,random(m)+1],M[,random(n)+1])))); my(D=FpM_det(N,1009)); if((D==1 || D==1008) && abs(matdet(N))==1, my(Ni=ZM_to_zm(N^-1)); my(P=zm_to_ZM(zm_mul(Ni,ZM_to_zm(M)))); my(R=N~*q*N); return([P,R])))); }; diag(M)=vector(#M,i,M[i,i]); qfreduce(q)= { my(n=vecmax(diag(q))); my(T=qflllgram(q),Q=T~*q*T); my(N=vecmax(diag(Q))); if (n>N,q=Q;n=N); for(i=1,n,my(D=reduce(q,i));if(D,return(D))); [qfminim(q,n),q]; }; order_aut(gram,fl=[5,2])= { my(un, ba, Rplus, rho, w, a, b, c, heights, ind_basis, rank, basis, vec_nei, bi); un=qfminim(gram,1)[3]; a=2^#un*((#un)!); if(#un==#gram,return(a)); if(#un, ba=matkerint(un~*gram); gram=ba~*gram*ba); my(init); [init,gram]=qfreduce(gram); Rplus= qfminim(gram,2)[3]; if( #Rplus==0, return(a*qfauto(qfisominit(gram,fl,init))[1]), rho=Rplus*matrix(#Rplus,1,i,j,1); w=gram*rho; b=qfauto(qfisominit([gram, w*w~],fl,init))[1]*1/2; heights=Rplus~*w; ind_basis=select(x->heights[x,1]==2,vector(#Rplus,u,u)); rank=#ind_basis; basis=vector(rank,i,Rplus[,ind_basis[i]]); vec_nei=vector(rank,u,List()); for(i=1,rank, bi=basis[i]~*gram; for(j=i+1,rank, if(bi*basis[j], listput(vec_nei[i],j); listput(vec_nei[j],i)) ) ); c=order_weyl(print_name(vec_to_mset(tree_comp(vec_nei)))); return(a*b*c); ); }; \\ \\ the invariant BV \\ graph(L)= { my(S); S=qfminim(L,3)[3]; if(#S==0,return(0)); S=ZM_to_F2m(S); L=ZM_to_F2m(L); return(F2m_to_Flm(F2m_mul(F2m_mul(F2m_transpose(S),L),S))); }; BV(L,p=1009)= { my( G=graph(L) ); if(G===0,return(G)); my(S,M); S=Flm_mul(G,G,p); M=vector(#S,i,vec_to_mset2(Col(S[,i]))); return(Set(M)); }; vec_to_mset2(v)= { Vec(ZM_to_zm(matreduce(v))); }; \\ hashed version of BV HBV(L)=hash_GEN(BV(L)); \\ \\ checking functions \\ \\ checking neighbor form : isotropy and ordering check_nf(v)= { my( [d,iso,s]= v[1..3] ); if( [d,iso,s] <> sort_iso([d,iso,s]), return(0) ); my( dd = (2-d%2)*d ); return( sum(i=1,#iso~, Mod(iso[i,1]^2,dd) )==Mod(0,dd) ); }; \\ check integrality and unimodularity of a gram matrix install("RgM_is_ZM","iG"); check_int_un_dim(gram,dim)= { if(#gram<>dim, return(0)); if(RgM_is_ZM(gram)==0, return(0)); return(matdet(gram)==1); }; \\ check distinct invariants check_distinct_invariants(v,depth)= { my( invar=fun_invar(depth) ); my( res=vector(#v) ); my(ind=0); print1(ind,"/",#res," ",Strchr(13)); parfor(i=1,#v, invar(v[i]), c, res[i]=c; ind++; print1(ind,"/",#res," ",Strchr(13)) ); return( #res==#Set(res) ) }; \\ check root systems check_root_systems(v)= { return( parapply(x->(root_system(nei(x))==x[5]),v) ); }; \\ check order automorphism groups check_order_aut(v)= { res=vector(#v); my(ind=0); print1(ind,"/",#res," ",Strchr(13)); parfor(i=1,#res, order_aut(nei(v[i]))==1/v[i][4], c, res[i]=c; ind++; print1(ind,"/",#res," ",Strchr(13)) ); return(res); };