Annexe au memoire "Formes automorphes et voisins de Kneser des reseaux de Niemeier"

Liste de programmes pour PARI/GP
--------------------------------

Version:

GP/PARI CALCULATOR Version 2.5.0 (released)
i386 running darwin (x86-64 kernel) 64-bit version
compiled: Jul  4 2011, gcc-4.2.1 (Apple Inc. build 5666) (dot 3)
(readline v6.2 enabled, extended help enabled)


(a) Justification de la table V.3
---------------------------------


G=[[1,-1,0,0,0,0,0,0], [0,1,-1,0,0,0,0,0], [0,0,1,-1,0,0,0,0], [0,0,0,1,-1,0,0,0], [0,0,0,0,1,-1,0,0], [0,0,0,0,0,1,-1,0], [0,0,0,0,0,0,1,-1], [0,0,0,0,0,0,1,1],[1/2,1/2,1/2,1/2,1/2,1/2,1/2,1/2]];

L=1/2*mathnf(2*matrix(8,9,i,j,G[j][i]));

RP=L*qfminim(mattranspose(L)*L)[3];

GRAM=mattranspose(RP)*RP; 

a1=RP;

a1a1=[]; for(i=1,#RP,for(j=i+1,#RP,if(GRAM[i,j]==0,a1a1=concat(a1a1,[[i,j]]))));

a2=[]; for(i=1,#RP,for(j=i+1,#RP,if(GRAM[i,j]==-1,a2=concat(a2,[[i,j]]))));

a1a1a1=[]; for(i=1,#a1a1,for(j=a1a1[i][2]+1,#RP,if(vecextract(GRAM,a1a1[i],[j])==[0;0],a1a1a1=concat(a1a1a1,[concat(a1a1[i],j)]))));

a1a2=[]; for(i=1,#a2,for(j=1,#RP,if(vecextract(GRAM,a2[i],[j])==[0;0],a1a2=concat(a1a2,[concat(a2[i],[j])]))));

a3=[]; for(i=1,#a1a1,for(j=1,#RP,if(vecextract(GRAM,a1a1[i],[j])==[-1;-1],a3=concat(a3,[concat(a1a1[i],[j])]))));

d4=[]; for(i=1,#a1a1a1,for(j=1,#RP,if(vecextract(GRAM,a1a1a1[i],[j])==[-1;-1;-1],d4=concat(d4,[concat(a1a1a1[i],j)]))));

a4=[]; a2double=[]; for(i=1,#RP,for(j=1,#RP,if(GRAM[i,j]==-1,a2double=concat(a2double,[[i,j]])))); for(i=1,#a1a1,for(j=1,#a2double,if(vecextract(GRAM,a1a1[i],a2double[j])==[-1,0;0,-1],a4=concat(a4,[concat(a1a1[i],a2double[j])]))));

a1a1a2=[]; for(i=1,#a2,for(j=1,#a1a1,if(vecextract(GRAM,a2[i],a1a1[j])==[0, 0;0, 0],a1a1a2=concat(a1a1a2,[concat(a2[i],a1a1[j])]))));

a2a2=[]; for(i=1,#a2,for(j=i+1,#a2,if(vecextract(GRAM,a2[i],a2[j])==[0, 0;0, 0],a2a2=concat(a2a2,[concat(a2[i],a2[j])]))));

a1a1a1a1=[]; for(i=1,#a1a1a1,for(j=a1a1a1[i][3]+1,#RP,if(vecextract(GRAM,a1a1a1[i],[j])==[0;0;0],a1a1a1a1=concat(a1a1a1a1,[concat(a1a1a1[i],j)]))));

a1a3=[]; for(i=1,#a3,for(j=1,#RP,if(vecextract(GRAM,a3[i],[j])==[0;0;0],a1a3=concat(a1a3,[concat([j],a3[i])]))));

% L est une matrice de Gram du reseau E_8 standard, note E_8 ci-dessous.
% RP est un systeme positif de l'ensemble des racines de E_8. Il s'agit donc de 120 racines numerotees de 1 a 120, a savoir vecextract(RP,[j]) pour j=1,...,120 (vecteurs colonnes).
% Pour chacun des systemes de racines R de rang <= 4 ci-dessus, on enumere les sous-reseaux de E_8 isometriques a Q( R), 
% en retournant une Z-base de chacun de ces reseaux qui est  constituee d'elements de RP.
% Concretement, le programme retourne un vecteur v, tel que pour i=1,...,#v, v[i] est lui-meme un vecteur, dont les coordonnees sont les indices des elements de RP formant la base choisie du i-eme reseau. 
% Le calcul se fait par induction sur le rang de R.
% Exemple : mattranspose(vecextract(RP,a2[1]))*vecextract(RP,a2[1]) retourne la matrice de Cartan de A2.
% Temps approximatifs, avec un processeur 2.13 GHz INtel Core 2 Duo : 
% a1a1a1 (4 min), a1a2 (2 min), d4 (15 sec.), a4 (2 min), a1a1a2 (5h30), a2a2 (13 min), a1a1a1a1 (1h), a1a3 (1h30)


lag=[1,I,0,0,0,0,0,0;0,0,1,I,0,0,0,0;0,0,0,0,1,I,0,0 ;0,0,0,0,0,0,1,I];

f(X)=vector(9,u,[u-1,sum(i=1,#X,matdet(lag*vecextract(RP,X[i]))^(u-1))]);

% Exemple : f(d4) renvoie [[0, 3150], [1, 0], [2, 0], [3, 0], [4, 4800], [5, 0], [6,-4800], [7, 0], [8, 43200]]




(suite) preuve du corollaire IX.5.15
------------------------------------

orthocomp(v)={A=[]; for(i=1,#RP,if(mattranspose(vecextract(RP,v))*vecextract(RP,[i])==0,A=concat(A,[i]))); rho=1/2*sum(i=1,#A,vecextract(RP,[A[i]])); B=[]; for(i=1,#A,if(mattranspose(rho)*vecextract(RP,[A[i]])==1,B=concat(B,[A[i]]))); B};

e6=[];
for(i=1,#a2,e6=concat(e6,[orthocomp(a2[i])]));

d6=[];
for(i=1,#a1a1,d6=concat(d6,[orthocomp(a1a1[i])]));

esp=[1,0,0,0,0,0,0,0;0,1,0,0,0,0,0,0;0,0,1,0,0,0,0,0;0,0,0,1,0,0,0,0;0,0,0,0,1,0,0,0;0,0,0,0,0,1,0,0];

g(X)=sum(i=1,#X,matdet(esp*vecextract(RP,X[i]))^2);

(g(d6)*2*6!*2^5)/(g(e6)*2*51840)

%idem que ci-dessus, avec e6 et d6, qui sont les orthogonaux respectifs des %a_2 et a1a1 (seuls systemes de racines ADE de rang 6 et d'indice de connexion 3 et 4)



(b) justification de la table C.2
----------------------------------

allocatemem(2^30);

A1A1=concat(vector(#a1a1,u,[vecextract(RP,[a1a1[u][1]])~,vecextract(RP,[a1a1[u][2]])~]),vector(#a1a1,u,[-vecextract(RP,[a1a1[u][1]])~,vecextract(RP,[a1a1[u][2]])~]));

A2=concat(concat(vector(#a2,u,[vecextract(RP,[a2[u][1]])~,vecextract(RP,[a2[u][2]])~]),vector(#a2,u,[vecextract(RP,[a2[u][1]])~+vecextract(RP,[a2[u][2]])~,-vecextract(RP,[a2[u][1]])~])),vector(#a2,u,[vecextract(RP,[a2[u][2]])~,-vecextract(RP,[a2[u][1]])~-vecextract(RP,[a2[u][2]])~]));

vnorm(d)={A=[]; B=L*qfminim(mattranspose(L)*L,d)[3]; for(i=1,matsize(B)[2],if(vecextract(B,[i])~*vecextract(B,[i])==d,A=concat(A,[vecextract(B,[i])]))); 
matrix(8,#A,i,j,A[j][i,1])};

N=vnorm(4);

{D7=[];  B=mattranspose(RP)*N;  for(i=1,#RP,for(j=1,#N,if((B[i,j]==1)||(B[i,j]==-1),D7=concat(D7,[[i,[-B[i,j],j]]]))));};

{D8=[];  B=mattranspose(RP)*N;  for(i=1,#RP,for(j=1,#N,if(B[i,j]==0,D8=concat(D8,[[i,[1,j]],[i,[-1,j]]]))));};

{D12=[];  B=mattranspose(N)*N;  for(i=1,#N,for(j=1,#N,if((B[i,j]==2)||(B[i,j]==-2),D12=concat(D12,[[i,[-B[i,j]/2,j]]]))));}

S=vnorm(6);

{D11=[];  B=mattranspose(RP)*S;  for(i=1,#RP,for(j=1,#S,if((B[i,j]==1)||(B[i,j]==-1),D11=concat(D11,[[i,[-B[i,j],j]]]))));};

{D12p=[];  B=mattranspose(RP)*S;  for(i=1,#RP,for(j=1,#S,if(B[i,j]==0,D12p=concat(D12p,[[i,[1,j]],[i,[-1,j]]]))));};

U=[2,I,I,I,I,0,0,0];

V=[0,0,0,I,-I,I,I,2];

P(x,y,j,k,u,v,w)=(((x*u~)*(v*y~)-(y*u~)*(v*x~))[1,1])^(k-4)*((w*x~)[1]*X+(w*y~)[1]*Y)^j;

sym(Q)=Q+subst(subst(subst(Q,X,Z),Y,X),Z,Y);

coeff(j,k,u,v,w,Q)=sym(sum(s=1,#Q,P(Q[s][1],Q[s][2],j,k,u,v,w)));

coeffd7(j,k,u,v,w)=sum(s=1,#D7,P(vecextract(RP,[D7[s][1]])~,D7[s][2][1]*vecextract(N,[D7[s][2][2]])~,j,k,u,v,w));

coeffd8(j,k,u,v,w)=sum(s=1,#D8,P(vecextract(RP,[D8[s][1]])~,D8[s][2][1]*vecextract(N,[D8[s][2][2]])~,j,k,u,v,w));

coeffd11(j,k,u,v,w)=sum(s=1,#D11,P(vecextract(RP,[D11[s][1]])~,D11[s][2][1]*vecextract(S,[D11[s][2][2]])~,j,k,u,v,w));

coeffd12(j,k,u,v,w)=sum(s=1,#D12,P(vecextract(N,[D12[s][1]])~,D12[s][2][1]*vecextract(N,[D12[s][2][2]])~,j,k,u,v,w));

coeffd12p(j,k,u,v,w)=sum(s=1,#D12p,P(vecextract(RP,[D12p[s][1]])~,D12p[s][2][1]*vecextract(S,[D12p[s][2][2]])~,j,k,u,v,w));

tau2(j,k)=coeffd12(j,k,U,V,U+V)/coeff(j,k,U,V,U+V,A2);

% utilise RP, L, a1a1 et a2 du (a)
% A1A1 liste les couples de racines orthogonales de E_8, modulo (x,y) -> (y,x) et (x,y) -> (-x,-y) (un quart des cas)
% A2 liste les couples de racines de E_8 de matrice de Gram la matrice de Cartan de A2, modulo les memes tranformations que ci-dessus.
% Si Q=A1A1 ou A2, coeff(j,k,u,v,w) donne la moitie du coefficient de Fourier de la forme f_{j,k,u,v,w} en la matrice de Cartan de Q,
% on exploite dans les deux cas les symetries (X,Y) -> (Y,X) et (-X,-Y).
% N est le vecteur constitue des elements x de E_8 tels que x.x=4, modulo x -> -x (il y en a 1080).  
% D7 liste les couples (x,y) d'elements de E_8 de matrice de Gram la matrice [2, -1; -1; 4], modulo (x,y) -> (-x,-y) (le calcul prend environs 13 minutes).
% pour des raisons de stockage, nous donnons plutot l'indice de x comme colonne de RP, et celui de y comme colonne de N (avec un signe supplementaire)
% D8 liste les couples (x,y) d'elements de E_8 de matrice de Gram la matrice [2, 0; 0, 4], modulo (x,y) -> (-x,-y) (12 min)
% pour des raisons de stockage, nous donnons plutot l'indice de x et y comme colonnes de N (avec un signe supplementaire pour y)
% idem pour D11, D12 et D12p, pour les matrices de Gram respectives [2, -1; -1, 6], [2, 0; 0, 6] et [4, -2; -2, 4] (resp. 1h05, 4h25 et 1h de calcul)
% coeffd7(j,k,u,v,w) (resp. coeffd8(j,k,u,v,w)) donne la moitie du coefficient de Fourier de la forme f_{j,k,u,v,w} en la matrice de Gram des couples d'elements de D7 (resp. D8)
% Exemple : 
% coeff(4,10,U,V,U+V,A2)
% -17625600*X^4 + 35251200*Y*X^3 - 52876800*Y^2*X^2 + 35251200*Y^3*X - 17625600*Y^4
% coeff(4,10,U,V,U+V,A2)/-17625600
% X^4 - 2*Y*X^3 + 3*Y^2*X^2 - 2*Y^3*X + Y^4
% factor(subst(%,Y,1))
% [X^2 - X + 1 2]
% coeff(4,10,U,V,U+V,A1A1)/-17625600
% -2*X^4 + 18*Y^2*X^2 - 2*Y^4
% coeffd7(4,10,U,V,U+V)/-17625600
% -16*X^4 + 32*Y*X^3 - 168*Y^2*X^2 + 152*Y^3*X + 112*Y^4
% coeffd8(4,10,U,V,U+V)/-17625600
% 36*X^4 - 252*Y^2*X^2 - 252*Y^4

calclambda(j,k,p,q)=2^(k-2)*(subst(subst(p,X,2*x-y),Y,y)+subst(subst(p,X,-x+2*y),Y,x)+subst(subst(p,X,x+y),Y,-x+y))/q+2^(j+2*k-4);

%Exemple : 
% -calclambda(6,8,-16*Y^2*(X-Y)*(X+Y)*(17*X^2+13*Y^2),y^2*x^2*(y-x)^2)*2;
% 409600


(c) Formules explicites de Mestre-Riemann-Weil
----------------------------------------------

g(t)=(1-abs(t))*cos(Pi*t)+sin(Pi*abs(t))/Pi;
F(t,l)=g(t/l)/cosh(t/2);
phi(x)=1/2*psi((x+1)/2)-1/2*psi(x/2);
R(x)=2*Pi^2*exp(-x)/((x^2+Pi^2)^2);
nn=100;
R1(w,l)=2*l*sum(n=0,nn,R(2*l*(1/2+w/4+n)));
R2(l)=2*l*sum(n=1,nn,(-1)^(n+1)*n*R(l*n));
R3(l)=4*l*sum(n=0,nn,(-1)^n*R(l*(n+1/2)));
J(w,l)=log(Pi)-real(psi(1/2+w/4+I*Pi/(2*l)))+1/Pi*imag(psi(1/2+w/4+I*Pi/(2*l)))-1/(2*l)*real(psi'(1/2+w/4+I*Pi/(2*l)))+R1(w,l);
diff1eps(l)=1+2*Pi/l*imag(phi(1+I*Pi/l))+2*Pi/l^2*imag(phi'(1+I*Pi/l))+R2(l);
triv(l)=1/2*(J(0,l)+diff1eps(l));
eps(l)=1/2*(J(0,l)-diff1eps(l));
Phi0(l)=l*8/Pi^2;
Phi12(l)=4*real(phi(1/2+I*Pi/l))-4/Pi*imag(phi(1/2+I*Pi/l))+4/l*real(phi'(1/2+I*Pi/l))+R3(l);

N=100;
isprimepower(n)=if(matsize(factor(n))[1]>1,0,factor(n));
primepowers(n)={ my(v); v=[]; 
for(i=2,n,if(isprimepower(i),v=concat(v,[[isprimepower(i)[1,1], 
isprimepower(i)[1,2]]]))); v};
pippower(n)=sum(i=1,n,#primepowers(i)-#primepowers(i-1));
liste=primepowers(N);
E4=1+240*sum(k=1,N,sigma(k,3)*q^k+O(q^(N+1)));
E6=1-504*sum(k=1,N,sigma(k,5)*q^k+O(q^(N+1)));
D11=q*prod(k=1,N-1,1-q^k+O(q^N))^24;
D15=E4*D11;
D17=E6*D11;
D19=E4^2*D11;
D21=E4*E6*D11;
hecke(F,p,w)=x^2-polcoeff(F,p)/p^(w/2)*x+1;
sym2hecke(F,p,w)=x^3+(polcoeff(F,p)^2/p^w-1)*(x-x^2)-1;
r1=[[1,0], vector(#liste,u,1)];
rD11=[[0,0,11], vector(#liste,u,polsym(hecke(D11,liste[u][1],11),liste[u][2]+1)[liste[u][2]+1])];
rD15=[[0,0,15], vector(#liste,u,polsym(hecke(D15,liste[u][1],15),liste[u][2]+1)[liste[u][2]+1])];
rD17=[[0,0,17], vector(#liste,u,polsym(hecke(D17,liste[u][1],17),liste[u][2]+1)[liste[u][2]+1])];
rD19=[[0,0,19], vector(#liste,u,polsym(hecke(D19,liste[u][1],19),liste[u][2]+1)[liste[u][2]+1])];
rD21=[[0,0,21], vector(#liste,u,polsym(hecke(D21,liste[u][1],21),liste[u][2]+1)[liste[u][2]+1])];
rSym2D11=[[0,1,22], vector(#liste,u,polsym(sym2hecke(D11,liste[u][1],11),liste[u][2]+1)[liste[u][2]+1])];
rr=[r1, rD11, rD15, rD17, rD19, rD21, rSym2D11];


biloo(p,q,l)=(p[1]*q[1]+p[2]*q[2])*triv(l)+(p[1]*q[2]+p[2]*q[1])*eps(l)+(p[1]+p[2])*sum(i=3,#q,J(q[i],l))+(q[1]+q[2])*sum(i=3,#p,J(p[i],l))+sum(i=3,#p,sum(j=3,#q,J(p[i]+q[j],l)+J(abs(p[i]-q[j]),l)));
mult(p,l)=if(biloo(p,p,l)>0,Phi0(l)/biloo(p,p,l),"vide");
epsil(p,q)=I^(p[1]*q[2]+p[2]*q[1]+(p[1]+p[2])*sum(i=3,#q,q[i]+1)+(q[1]+q[2])*sum(i=3,#p,p[i]+1)+sum(i=3,#p,sum(j=3,#q,2*max(p[i],q[j])+2)));
vecepsil(p)=vector(#rr,i,epsil(p,rr[i][1]));
C(p,l)=Phi0(l)-biloo(p,p,l);
t(p,q,l)=-Phi12(l)/2*(1-epsil(p,q))/2-biloo(p,q,l)+sqrt(C(p,l)*C(q,l));
tvec(p,l)=vector(#rr,i,t(p,rr[i][1],l));


% p est un vecteur de longueur >=2 reprˇsentant un ˇlˇment de K_infty, l'element V = a 1 + b eps_C/R + I_{w_1} + I_{w_2} + .... + I_{w_k} etant reprˇsentˇ par p = [a,b,w_1,w_2,....,w_k], 
% exemple : I_{11} est reprˇsentˇ par p = [0,0,11],  et eps_C/R + I_22 + I_10 + I_0 par [0,1,20,0,10]

%%% Dˇmonstration des cas n=2 ou w=0 du thˇor¸me IX.3.2
% \p 2;
% vector(21,w,J(w,log(2))) renvoie [0.85, 0.73, 0.61, 0.51, 0.41, 0.32, 0.23, 0.15, 0.074, 0.0028, -0.065, -0.13, -0.19, -0.25, -0.30, -0.35, -0.40, -0.45, -0.50, -0.54, -0.59], cela justifie la table IX.3
% [triv(log(2)), eps(log(2))] renvoie [0.56, 0.42], cela dˇmontre le thˇor¸me IX.3.2 pour w=0.
% Phi12(log(2))/2 renvoie 0.28, cela montre m(I_{13})=0.
% vector(6,i,mult([0,0,2*i+9],log(8))) renvoie [1.2, 1.3, 1.5, 1.7, 1.9, 2.1], i.e. m(I_w)<=1 si w<=19 et m(I_{21})<=2. 
% a(l)=1/2*Phi12(l)+J(21,l);   
% b(l)=sqrt((Phi0(l)-triv(l))*(Phi0(l)/2-biloo([0,0,21],[0,0,21],l))); 
% a(log(6)) renvoie 0.17 et b(log(6)) renvoie 0.13, cela montre m(I_{21}) <2. 


gramoowodd(w,l)=matrix((w+1)/2,(w+1)/2,u,v,biloo([0,0,2*u-1],[0,0,2*v-1],l));
triodd(M)={my(RES); RES=[]; for(j=1,matsize(M)[2],if(sum(i=1,matsize(M)[1],M[i,j]+abs(M[i,j]))*sum(i=1,matsize(M)[1],M[i,j]-abs(M[i,j]))==0, RES=concat(RES,[j]))); N=abs(vecextract(M,RES))};
convertodd(M)={my(u); u=vector(#M,j,[0,0]); for(k=1,#M,for(i=1,#(M~),u[k]=concat(u[k],vector(abs(M[i,k]),j,2*i-1)))); u};
n1mat(M)=vecmax(vector(matsize(M)[2],j,sum(i=1,matsize(M)[1],M[i,j])));

% Dˇmonstration du cas w impair <=19 du thˇor¸me IX.3.2
% \p 21
% B=gramoowodd(19,log(9));
% qfsign(B) renvoie [10,0] : B est dˇfinie positive.
% [vecmin(abs(B)),vecmax(abs(B))] renvoie [0.0118399918510793214589, 3.47690749750616665860].
% A=round(10^4*B)/10^4;
% qfsign(A) renvoie [10,1]
% n1mat(A^-1)*1.0 renvoie la norme de A^-1 subordonnˇe a la norme sup, on trouve 3.22005138412652261362
% Cela dˇmontre le lemme IX.3.19.
% 10^4*(1-n1mat(A^-1)*#A*10^-4)^-1*Phi0(log(9)) renvoie 17867.5660731641503170 < 17868
% M=qfminim(10^4*A,17867)[3];
% #M renvoie 24 : 24 couples de solutions opposˇes.
% convertodd(triodd(M)) renvoie un vecteur constituˇ des solutions sous formes de vecteurs "p", en ne retenant que les effectifs, i.e. le vecteur de taille 11 suivant :
% [[0, 0, 15], [0, 0, 13], [0, 0, 13, 19], [0, 0, 11], [0, 0, 11, 19], [0, 0, 9], [0, 0, 9, 19], [0, 0, 7, 19], [0, 0, 5, 19], [0, 0, 17], [0, 0, 19]].  
% cela dˇmontre le lemme IX.3.21
% vector(5,i,mult([0,0,19,2*i+3],log(9))) renvoie (prˇcision \p 3) [1.08, 1.22, 1.26, 1.19, 1.05] : c'est le lemme IX.3.23
% matrix(5,6,i,j,t([0,0,19,15-2*i],rr[j][1],log(5))) donne la table IX.4 (prˇcisions \p 4)

gramooweven(w,l)=matrix(w/2+2,w/2+2,u,v,if((u<(w/2+1))&&(v<(w/2+1)),biloo([0,0,2*u],[0,0,2*v],l),if((u>w/2)&&(v>w/2),if(u==v,triv(l),eps(l)),J(min(2*u,2*v),l))));
trieven(M)={my(RES); RES=[]; for(j=1,matsize(M)[2],if((sum(i=1,matsize(M)[1],M[i,j]+abs(M[i,j]))*sum(i=1,matsize(M)[1],M[i,j]-abs(M[i,j]))==0)&&(sum(i=1,matsize(M)[1]-2,Mod(M[i,j],2))==Mod(M[matsize(M)[1],j],2)), RES=concat(RES,[j]))); abs(vecextract(M,RES))};
converteven(M)={my(u); u=vector(#M,j,[M[#(M~)-1,j],M[#(M~),j]]); for(k=1,#M,for(i=1,#(M~)-2,u[k]=concat(u[k],vector(abs(M[i,k]),j,2*i)))); u};

% Dˇmonstration du cas w pair <=20 du thˇor¸me IX.3.2
% \p 21
% A=round(10^4*gramooweven(20,log(9)))/10^4;
% [vecmin(abs(A)),vecmax(abs(A))]*1.0 renvoie [0.0172, 3.03]
% qfsign(A) renvoie [12,0]
% 1.0*n1mat(A^-1)*#A*10^-4 renvoie 0.0137574147093330776015 < 1 : Gram(20,log(9)) est dˇfinie positive
% 10^4*(1-n1mat(A^-1)*#A*10^-4)^(-1)*Phi0(log(9)) renvoie 18058.4694454765267122 < 18059
% M=qfminim(10^4*A,18058)[3]; renvoie 73 couples solutions
% poids=converteven(trieven(M)) renvoie les solutions effectives de determinant 1, i.e. la liste des 15 reprˇsentations suivantes 
%[[0, 0, 4, 20], [0, 0, 6, 20], [0, 1, 18], [0, 1, 20], [1, 0, 16, 20], [1, 0, 14, 20], [1, 0, 12, 20], [1, 0, 10, 20], [1, 0], [0, 0, 14, 20], [0, 0, 12, 20], [0, 0, 10, 18], [0, 0, 10, 20], [0, 0, 8, 18], [0, 0, 8, 20]]
% vecmax(vector(#poids,i,mult(poids[i],log(9)))) renvoie 1.58 (\p 3), cela termine la dˇmonstration du lemme IX.3.27.
% [C([0,0,18,10],log(10)),C([1,0,20,10],log(16)),C([1,0,20,16],log(16))] renvoie  [-0.00012, -0.13, -0.012].
% for(i=1,#poids,print([poids[i],vector(3,u,t(poids[i],rr[u][1],log(5)))])) justifie la table IX.5

% vˇrification du lemme 3.32

bilf(v,w,l)=sum(u=1,pippower(round(exp(l))-1),log(liste[u][1])*F(liste[u][2]*log(liste[u][1]),l)*v[u]*w[u]/(liste[u][1]^(liste[u][2]/2)));
Ci(p,r,l)=(epsil(p,r[1])-1)/4*Phi12(l)-biloo(p,r[1],l);
s(p,r,l)=Ci(p,r,l)+sqrt(C(p,l)*bilf(r[2],r[2],l));
svec(p,l)=vector(#rr,j,s(p,rr[j],l));
CI(p,E,l)=vector(#E,u,Ci(p,rr[E[u]],l));
gram(E,l)=matrix(#E,#E,i,j,bilf(rr[E[i]][2],rr[E[j]][2],l));
geom(p,E,l)={my(c,g); c=CI(p,E,l); g=gram(E,l); [(g^-1*c~)~, -c*g^-1*c~+C(p,l)]};
geom2(p,E,l)={my(c,g,h); c=CI(p,E,l); g=gram(E,l); h=[g[2,2], -g[1,2]; -g[2,1], g[1,1]];  if(#E==2,[(h*c~)~, -c*h*c~+matdet(g)*C(p,l)])};
positivite(q)=if(vecmax([vecmax(q[1]),q[2]])<0,0,1);
matgeom2(p,l)=matrix(#rr,#rr,i,j,if(i<>j,positivite(geom2(p,[i,j],l)),1));

% geom(p,E,l) est le test du scholie IX.3.31 avec p <-> V_0 et E une partie de {1,...,7} (indice des elements de R).
% geom2(p,E,l) renvoie det(g)*geom(p,E,l) si #E=2 (comme ca il n'est pas necessaire d'inverser la matrice g, noter que det(g)>0), rien si #E<>2.
% mat2geom(p,l) renvoie une matrice dont le coefficient (i,j) avec i>j est 0 si et seulement si le critere geometrique (ii) est viole par p,rr[i],rr[j] pour lambda=l.
% geom2([0,0,21,17,7],[3,4],log(14)) renvoie [[-0.345, -0.510], -0.00769]
% geom2([0,0,22,12],[2,3],log(5)) renvoie [[-0.0429, -0.0234], -0.00175]
% geom2([1,0,22,12],[1,2],log(5)) renvoie [[-0.0442, -0.0549], -0.00377]
% geom2([0,1,22,12,10],[5,7],log(38)) renvoie [[-0.261, -0.332], -1.62]
% geom2([0,0,22,4],[1,6],log(7)) renvoie [[-0.138, -0.142], -0.0216]
% geom2([1,0,22,16],[1,4],log(8)) renvoie [[-0.255, -0.194], -0.00564]
% geom2([0,1,22,16,10],[2,3],log(9)) renvoie [[-0.277, -0.135], -0.00245]
% geom2([0,0,22,20,14,4],[1,6],log(40)) renvoie [[-0.880, -1.84], -0.0214]
% tous ces resultats ont toutes leurs coordonnees <0, cela justife la table IX.6

triselodd(M,v)={my(RES); RES=[]; for(j=1,matsize(M)[2],if((sum(i=1,matsize(M)[1],M[i,j]+abs(M[i,j]))*sum(i=1,matsize(M)[1],M[i,j]-abs(M[i,j]))==0)&&(M[(v+1)/2,j]), RES=concat(RES,[j]))); abs(vecextract(M,RES))};
dim(p)=p[1]+p[2]+2*(#p-2);
tripardim(p,d)=select(x->dim(x)==d,p);

% Dˇmonstration du cas w = 21 du thˇor¸me IX.3.2
% \p 21
% for(i=14,40,print([i,#qfminim(gramoowodd(21,log(i)),Phi0(log(i)),2000,2)[3]])) explique le choix l=log(28)
% A=round(10^4*gramoowodd(21,log(28)))/10^4;
% [vecmin(abs(A)),vecmax(abs(A))]*1.0 renvoie [0.033, 4.0] (\p 2)
% n1mat(A^-1)*1.0 renvoie 6.43356770633294756205
% 10^4*(1-n1mat(A^-1)*#A*10^-4)^-1*Phi0(log(28)) renvoie 27202.3413238888608080 < 27203
% M=qfminim(10^4*A,27202)[3];
% #M renvoie 73
% poids=convertodd(triselodd(M,21)); 
% poids est un vecteur constitue des solutions sous formes de vecteurs "p", en ne retenant que les effectifs contenant I_{21}, i.e. le vecteur de taille 29 suivant : 
% [[0, 0, 3, 21], [0, 0, 3, 19, 21], [0, 0, 3, 17, 21], [0, 0, 3, 13, 21], [0, 0, 3, 13, 19, 21], [0, 0, 3, 15, 21], [0, 0, 5, 21], [0, 0, 5, 19, 21], [0, 0, 5, 17, 21], [0, 0, 5, 11, 19, 21], [0, 0, 5, 13, 21], [0, 0, 5, 13, 19, 21], [0, 0, 5, 15, 21], [0, 0, 5, 15, 19, 21], [0, 0, 7, 15, 21], [0, 0, 9, 15, 21], [0, 0, 15, 21], [0, 0, 7, 13, 21], [0, 0, 13, 21], [0, 0, 11, 17, 21], [0, 0, 11, 21], [0, 0, 7, 17, 21], [0, 0, 7, 19, 21], [0, 0, 7, 21], [0, 0, 9, 17, 21], [0, 0, 9, 19, 21], [0, 0, 9, 21], [0, 0, 17, 21], [0, 0, 21]]
% En plus de la liste des 26 ˇlˇments du lemme IX.3.34, on trouve donc [0,0,21] (normal) mais aussi 
% select(x->(mult(x,log(28))<1),poids) qui renvoie [[0, 0, 5, 11, 19, 21], [0, 0, 5, 15, 19, 21]]
% cela dˇmontre le lemme IX.3.34
% vector(#poids,u,vecmin(tvec(poids[u],log(27)))) renvoie [-1.22, -1.39, -0.243, -0.575, -0.222, -1.12, 0.680, -0.227, -0.324, -1.01, -0.947, -0.534, -0.0634, -0.877, -0.209, -0.603, -0.254, -0.507, 0.595, -0.839, -0.133, 0.0593, -0.229, -0.0970, -0.205, -0.454, 0.636, -0.783, 0.295]
% select(x->(vecmin(tvec(x,log(27)))>0.05),poids) renvoie [[0, 0, 5, 21], [0, 0, 13, 21], [0, 0, 7, 17, 21], [0, 0, 9, 21], [0, 0, 21]]
% cela dˇmontre le lemme IX.3.35



triseleven(M,v)={my(RES); RES=[]; for(j=1,matsize(M)[2],if((sum(i=1,matsize(M)[1],M[i,j]+abs(M[i,j]))*sum(i=1,matsize(M)[1],M[i,j]-abs(M[i,j]))==0)&&(sum(i=1,matsize(M)[1]-2,Mod(M[i,j],2))==Mod(M[matsize(M)[1],j],2))&&(M[v/2,j]), RES=concat(RES,[j]))); abs(vecextract(M,RES))};
remult(p,l)=if(biloo(p,p,l)>0,Phi0(l)/biloo(p,p,l),10);
mopt(p)=floor(vecmin(vector(10,u,remult(p,log(2+5*u)))));

% Dˇmonstration du cas w = 22 du thˇor¸me IX.3.2
% for(i=70,90,print([i,#qfminim(gramooweven(22,log(i)),Phi0(log(i)),2000,2)[3]])) explique le choix l=log(80)
% A=round(10^6*gramooweven(22,log(80)))/10^6;
% [vecmin(abs(A)),vecmax(abs(A))]*1.0 renvoie [0.009718, 3.665269]
% n1mat(A^-1)*#A*10^-6*1.0 renvoie 0.000267806759662115755539
% 10^6*(1-n1mat(A^-1)*#A*10^-6)^-1*Phi0(log(80)) renvoie 3552888.49057336309452 < 3552889
% M=qfminim(10^6*A,3552888)[3];
% #M renvoie 701
% poids=converteven(triseleven(M,22));
% #poids renvoie 158
% poidsms1=select(x->(remult(x,log(77))>1.95),poids); 
% poidsms1 renvoie [[0, 0, 10, 22], [0, 0, 12, 22], [0, 0, 8, 22]];
% #select(x->(remult(x,log(77))<2.2),poidsms1) renvoie 3; cela demontre la premiere assertion du lemme IX.3.36.
% poidstri=select(x->(vecmin(tvec(x,log(77)))>-0.001),poids);
% poidstri renvoie 8 un vecteur de taille 8, a savoir [[0, 0, 4, 22], [0, 0, 4, 14, 20, 22], [0, 1, 10, 16, 22], [0, 1, 10, 20, 22], [1, 0, 12, 22], [0, 0, 12, 22], [0, 0, 8, 22], [1, 0, 16, 22]];
% cela termine la dˇmonstration du lemme IX.3.36
% t([0,0,22,8],[0,0,11],log(5)) renvoie -0.0244 (sous \p 3)

t2(p,q,l)=-biloo(p,q,l)+sqrt((Phi0(l)/2-biloo(p,p,l))*C(q,l));

% t2([0,0,22,8],[0,0,21,9],log(22)) renvoie -0.14 (sous \p 2)
% t2([0,0,22,10],[0,0,21,9],log(22)) renvoie -0.029
% t2([0,0,22,12],[0,0,21,13],log(22)) renvoie -0.23
% cela termine la dˇmonstration du thˇor¸me IX.3.2

RPi=[[1,0], [0,0,11], [0,0,15], [0,0,17], [0,0,19], [0,0,21], [0,0,19,7], [0,0,21,5], [0,0,21,9], [0,0,21,13], [0,1,22]];
t12(p,l)=-Phi12(l)-biloo(p,[1,0],l)+sqrt(C([1,0],l)*C(p,l));

% vector(#RPi-1,u,t12(RPi[u+1],log(4))) renvoie [-1.1, -0.64, -0.49, -0.35, -0.23, -0.79, -0.86, -0.35, -0.050, -0.82], et prouve la proposition IX.3.38

% vˇrification de la proposition IX.3.39
% A=round(10^4*gramoowodd(23,9.74))/10^4;
% n1mat(A^-1)*#A*10^-4*1.0 renvoie 0.0614 < 1
% qfsign(A) renvoie [12,0], ce qui conclut.

(d) invariants sous SO(L), L=E_8+E_8 ou E_16
--------------------------------------------

A=[[x^8 + 6*x^7 + 14*x^6 + 14*x^5 - 14*x^3 - 14*x^2 - 6*x - 1, 1/5806080], [x^8 + 4*x^7 + 6*x^6 + 4*x^5 - 4*x^3 - 6*x^2 - 4*x - 1, 1/15360], [x^8 + 3*x^7 + 5*x^6 + 4*x^5 - 4*x^3 - 5*x^2 - 3*x - 1, 1/2592], [x^8 + 3*x^7 + 4*x^6 + 3*x^5 - 3*x^3 - 4*x^2 - 3*x - 1, 1/576], [x^8 + 3*x^7 + 2*x^6 - x^5 + x^3 - 2*x^2 - 3*x - 1, 1/8640], [x^8 + 2*x^7 + 2*x^6 + 2*x^5 - 2*x^3 - 2*x^2 - 2*x - 1, 1/1536], [x^8 + 2*x^7 + 2*x^6 + x^5 - x^3 - 2*x^2 - 2*x - 1, 1/60], [x^8 + 2*x^7 + x^6 - x^2 - 2*x - 1, 1/144], [x^8 + 2*x^7 - 2*x^5 + 2*x^3 - 2*x - 1, 1/384], [x^8 + 2*x^7 - 2*x^6 - 6*x^5 + 6*x^3 + 2*x^2 - 2*x - 1, 1/18432], [x^8 + x^7 + x^6 + x^5 - x^3 - x^2 - x - 1, 1/40], [x^8 + x^7 + x^6 - x^2 - x - 1, 1/288], [x^8 + x^7 + x^5 - x^3 - x - 1, 1/192], [x^8 + x^7 - x - 1, 1/28], [x^8 + x^7 - x^5 + x^3 - x - 1, 1/48], [x^8 + x^7 - x^6 - x^5 + x^3 + x^2 - x - 1, 1/120], [x^8 + x^7 - x^6 - 2*x^5 + 2*x^3 + x^2 - x - 1, 1/48], [x^8 + x^7 - 2*x^6 - 3*x^5 + 3*x^3 + 2*x^2 - x - 1, 1/576], [x^8 + 2*x^6 - 2*x^2 - 1, 1/768], [x^8 + x^6 - x^2 - 1, 1/72], [x^8 - 1, 5/64], [x^8 - x^6 + 2*x^5 - 2*x^3 + x^2 - 1, 1/432], [x^8 - x^6 + x^5 - x^3 + x^2 - 1, 1/36], [x^8 - x^6 + x^2 - 1, 1/24], [x^8 - x^6 - x^5 + x^3 + x^2 - 1, 1/36], [x^8 - x^6 - 2*x^5 + 2*x^3 + x^2 - 1, 1/432], [x^8 - 2*x^6 + 2*x^2 - 1, 19/4608], [x^8 - x^7 + x^6 - x^2 + x - 1, 1/288], [x^8 - x^7 + x^6 - x^5 + x^3 - x^2 + x - 1, 1/40], [x^8 - x^7 + x^5 - x^3 + x - 1, 1/48], [x^8 - x^7 + x - 1, 1/28], [x^8 - x^7 - x^5 + x^3 + x - 1, 1/192], [x^8 - x^7 - x^6 + 2*x^5 - 2*x^3 + x^2 + x - 1, 1/48], [x^8 - x^7 - x^6 + x^5 - x^3 + x^2 + x - 1, 1/120], [x^8 - x^7 - 2*x^6 + 3*x^5 - 3*x^3 + 2*x^2 + x - 1, 1/576], [x^8 - 2*x^7 + 2*x^6 - x^5 + x^3 - 2*x^2 + 2*x - 1, 1/60], [x^8 - 2*x^7 + 2*x^6 - 2*x^5 + 2*x^3 - 2*x^2 + 2*x - 1, 1/1536], [x^8 - 2*x^7 + x^6 - x^2 + 2*x - 1, 1/144], [x^8 - 2*x^7 + 2*x^5 - 2*x^3 + 2*x - 1, 1/384], [x^8 - 2*x^7 - 2*x^6 + 6*x^5 - 6*x^3 + 2*x^2 + 2*x - 1, 1/18432], [x^8 - 3*x^7 + 5*x^6 - 4*x^5 + 4*x^3 - 5*x^2 + 3*x - 1, 1/2592], [x^8 - 3*x^7 + 4*x^6 - 3*x^5 + 3*x^3 - 4*x^2 + 3*x - 1, 1/576], [x^8 - 3*x^7 + 2*x^6 + x^5 - x^3 - 2*x^2 + 3*x - 1, 1/8640], [x^8 - 4*x^7 + 6*x^6 - 4*x^5 + 4*x^3 - 6*x^2 + 4*x - 1, 1/15360], [x^8 - 6*x^7 + 14*x^6 - 14*x^5 + 14*x^3 - 14*x^2 + 6*x - 1, 1/5806080], [x^8 - 8*x^7 + 28*x^6 - 56*x^5 + 70*x^4 - 56*x^3 + 28*x^2 - 8*x + 1, 1/696729600], [x^8 - 5*x^7 + 10*x^6 - 11*x^5 + 10*x^4 - 11*x^3 + 10*x^2 - 5*x + 1, 1/311040], [x^8 - 4*x^7 + 4*x^6 + 4*x^5 - 10*x^4 + 4*x^3 + 4*x^2 - 4*x + 1, 1/184320], [x^8 - 4*x^7 + 8*x^6 - 12*x^5 + 14*x^4 - 12*x^3 + 8*x^2 - 4*x + 1, 1/18432], [x^8 - 4*x^7 + 10*x^6 - 16*x^5 + 19*x^4 - 16*x^3 + 10*x^2 - 4*x + 1, 1/155520], [x^8 - 3*x^7 + 2*x^6 + 3*x^5 - 6*x^4 + 3*x^3 + 2*x^2 - 3*x + 1, 1/6912], [x^8 - 3*x^7 + 3*x^6 - x^5 - x^3 + 3*x^2 - 3*x + 1, 1/1200], [x^8 - 3*x^7 + 5*x^6 - 6*x^5 + 6*x^4 - 6*x^3 + 5*x^2 - 3*x + 1, 1/864], [x^8 - 2*x^7 + 2*x^5 - 2*x^4 + 2*x^3 - 2*x + 1, 1/768], [x^8 - 2*x^7 + x^6 - 2*x^5 + 4*x^4 - 2*x^3 + x^2 - 2*x + 1, 1/2592], [x^8 - 2*x^7 + x^6 + x^5 - 2*x^4 + x^3 + x^2 - 2*x + 1, 1/108], [x^8 - 2*x^7 + x^6 + 2*x^5 - 4*x^4 + 2*x^3 + x^2 - 2*x + 1, 1/288], [x^8 - 2*x^7 + 2*x^6 - 2*x^5 + 2*x^4 - 2*x^3 + 2*x^2 - 2*x + 1, 1/128], [x^8 - 2*x^7 + 2*x^6 - x^4 + 2*x^2 - 2*x + 1, 1/288], [x^8 - 2*x^7 + 3*x^6 - 4*x^5 + 5*x^4 - 4*x^3 + 3*x^2 - 2*x + 1, 1/600], [x^8 - x^7 - 2*x^6 + x^5 + 2*x^4 + x^3 - 2*x^2 - x + 1, 1/1152], [x^8 - x^7 - x^6 + 2*x^4 - x^2 - x + 1, 1/144], [x^8 - x^7 - x^6 + x^5 + x^3 - x^2 - x + 1, 1/80], [x^8 - x^7 - x + 1, 1/28], [x^8 - x^7 + x^5 - 2*x^4 + x^3 - x + 1, 1/96], [x^8 - x^7 + x^5 - x^4 + x^3 - x + 1, 1/30], [x^8 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1, 1/54], [x^8 - x^7 + x^6 + 2*x^5 - 2*x^4 + 2*x^3 + x^2 - x + 1, 1/7776], [x^8 - x^7 + 2*x^6 - 3*x^5 + 2*x^4 - 3*x^3 + 2*x^2 - x + 1, 1/576], [x^8 - 4*x^6 + 6*x^4 - 4*x^2 + 1, 37/221184], [x^8 - 2*x^6 + 2*x^4 - 2*x^2 + 1, 1/64], [x^8 - 2*x^6 + 3*x^4 - 2*x^2 + 1, 1/288], [x^8 - x^6 - x^2 + 1, 1/36], [x^8 - x^6 + x^4 - x^2 + 1, 1/20], [x^8 - x^5 - x^3 + 1, 1/60], [x^8 - 2*x^4 + 1, 9/1024], [x^8 - x^4 + 1, 1/24], [x^8 + 2*x^4 + 1, 1/192], [x^8 + x^5 + x^3 + 1, 1/60], [x^8 + x^6 + x^2 + 1, 1/72], [x^8 + 2*x^6 + 3*x^4 + 2*x^2 + 1, 1/1728], [x^8 + 4*x^6 + 6*x^4 + 4*x^2 + 1, 1/46080], [x^8 + x^7 - 2*x^6 - x^5 + 2*x^4 - x^3 - 2*x^2 + x + 1, 1/1152], [x^8 + x^7 - x^6 - x^5 - x^3 - x^2 + x + 1, 1/80], [x^8 + x^7 - x^6 + 2*x^4 - x^2 + x + 1, 1/144], [x^8 + x^7 - x^5 - 2*x^4 - x^3 + x + 1, 1/96], [x^8 + x^7 - x^5 - x^4 - x^3 + x + 1, 1/30], [x^8 + x^7 + x + 1, 1/28], [x^8 + x^7 + x^6 - 2*x^5 - 2*x^4 - 2*x^3 + x^2 + x + 1, 1/7776], [x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1, 1/54], [x^8 + x^7 + 2*x^6 + 3*x^5 + 2*x^4 + 3*x^3 + 2*x^2 + x + 1, 1/576], [x^8 + 2*x^7 - 2*x^5 - 2*x^4 - 2*x^3 + 2*x + 1, 1/768], [x^8 + 2*x^7 + x^6 - 2*x^5 - 4*x^4 - 2*x^3 + x^2 + 2*x + 1, 1/288], [x^8 + 2*x^7 + x^6 - x^5 - 2*x^4 - x^3 + x^2 + 2*x + 1, 1/108], [x^8 + 2*x^7 + x^6 + 2*x^5 + 4*x^4 + 2*x^3 + x^2 + 2*x + 1, 1/2592], [x^8 + 2*x^7 + 2*x^6 - x^4 + 2*x^2 + 2*x + 1, 1/288], [x^8 + 2*x^7 + 2*x^6 + 2*x^5 + 2*x^4 + 2*x^3 + 2*x^2 + 2*x + 1, 1/128], [x^8 + 2*x^7 + 3*x^6 + 4*x^5 + 5*x^4 + 4*x^3 + 3*x^2 + 2*x + 1, 1/600], [x^8 + 3*x^7 + 2*x^6 - 3*x^5 - 6*x^4 - 3*x^3 + 2*x^2 + 3*x + 1, 1/6912], [x^8 + 3*x^7 + 3*x^6 + x^5 + x^3 + 3*x^2 + 3*x + 1, 1/1200], [x^8 + 3*x^7 + 5*x^6 + 6*x^5 + 6*x^4 + 6*x^3 + 5*x^2 + 3*x + 1, 1/864], [x^8 + 4*x^7 + 4*x^6 - 4*x^5 - 10*x^4 - 4*x^3 + 4*x^2 + 4*x + 1, 1/184320], [x^8 + 4*x^7 + 8*x^6 + 12*x^5 + 14*x^4 + 12*x^3 + 8*x^2 + 4*x + 1, 1/18432], [x^8 + 4*x^7 + 10*x^6 + 16*x^5 + 19*x^4 + 16*x^3 + 10*x^2 + 4*x + 1, 1/155520], [x^8 + 5*x^7 + 10*x^6 + 11*x^5 + 10*x^4 + 11*x^3 + 10*x^2 + 5*x + 1, 1/311040], [x^8 + 8*x^7 + 28*x^6 + 56*x^5 + 70*x^4 + 56*x^3 + 28*x^2 + 8*x + 1, 1/696729600]];

allocatemem(2^30);
B=[];
for(i=1,#A,for(j=1,#A,B=concat([[A[i][1]*A[j][1],A[i][2]*A[j][2]/2]],B)));
for(i=1,#A,B=concat(B,[[subst(A[i][1],x,x^2),A[i][2]/2]]));
C=vecsort(vector(#B,u,[B[u],subst(B[u][1],x,10)]),[2]);
C=vector(#B,u,C[u][1]);
regroupe(l)={FIN=[[],l[1]]; for(i=2,#l,if(l[i][1]==l[i-1][1],FIN[2]=[FIN[2][1],FIN[2][2]+l[i][2]],FIN=[concat(FIN[1],[FIN[2]]),l[i]])); concat(FIN[1],[FIN[2]])};
D=regroupe(C);
E=[];
for(i=1,#D,if(polcoeff(D[i][1],0)==1,E=concat(E,[[D[i][1],2*D[i][2]]])));

%A= liste des polynomes caracteristiques des elements de O(E_8) avec leur proportion (merci a Thomas Megarbane pour cette liste).
%E= liste des polynomes caracteristiques des elements de SO(E_8+E_8) avec leur proportion

indice(p)={v=[[p[1],1]]; for(i=2,#p,if(p[i]==v[#v][1],v[#v][2]=v[#v][2]+1,v=concat(v,[[p[i],1]]))); v};
ord(p)=1/prod(i=1,#p,p[i][2]!*p[i][1]^(p[i][2]));
base2(i,n)=-2*concat(vector(n-#binary(i)),binary(i))+vector(n,u,1);
base2even(i,n)=concat(base2(i,n-1),prod(j=1,n-1,base2(i,n-1)[j]));
carw(p)=vector(2^(#p-1),u,[prod(k=1,#p-1,x^(p[k])-base2(u-1,#p-1)[k])*(x^(p[#p])-prod(k=1,#p-1,base2(u-1,#p-1)[k])),ord(indice(p))*2^(1-#p)]);
r=16;
{B=[]; liste=partitions(r); for(i=1,#liste,B=concat(B,carw(liste[i]))); B};
C=vecsort(vector(#B,u,[B[u],subst(B[u][1],x,10)]),[2]);
C=vector(#B,u,C[u][1]);
regroupe(l)={FIN=[[],l[1]]; for(i=2,#l,if(l[i][1]==l[i-1][1],FIN[2]=[FIN[2][1],FIN[2][2]+l[i][2]],FIN=[concat(FIN[1],[FIN[2]]),l[i]])); concat(FIN[1],[FIN[2]])};
D=regroupe(C);
E=[];
for(i=1,#D,if(polcoeff(D[i][1],0)==1,E=concat(E,[[D[i][1],2*D[i][2]]])));

%E= liste des polynomes caracteristiques des elements de SO(E_16) avec leur proportion