/* Dans cette section on definit plusieurs functions qui seront utilisees pour le programe */ /* Donne la valeur absolue d'un entier */ abs:=function(x) return AbsoluteValue(x); end function; /* It is the same function as Generators but it really gives matrices and not elements of PSL_2(Z) */ Generators2:=function(G) gen:=Generators(G); n:=#gen; gen2:=[]; for i in [1..n] do gen2[i]:=Matrix(Integers(),2,2,Eltseq(gen[i])); end for; return gen2; end function; /* Soit "table" une table donnee, alors cette fonction construit une nouvelle table de dimension n x dimensions(table) */ expansion_table:=function(n,table) v:=[]; for i in [1..n] do v[i]:=table; end for; return v; end function; /* Soit v=[n1,n2,...,nr] une liste donnee alors cette fonction construit une table de dimension n1 x n2 x ... nr. */ table:=function(v) n:=#v; w:=[]; for i in [1..v[n]] do w:=Append(w,[]); end for; if n eq 1 then return w; else for j in [1..n-1] do w:=expansion_table(v[n-j],w); end for; end if; return w; end function; /* calcul les dimensions d'une table */ dimension:=function(v) n:=#v; w:=[]; w[1]:=n; x:=v; if #x[n] eq 0 then return w; else i:=1; while ( Type(x[n]) eq SeqEnum and #x[n] ne 0 ) do i:=i+1; x:=x[n]; w[i]:=#x; n:=w[i]; end while; end if; return w; end function; /* Soit v=[l_1,l_2,...,l_r] une liste donnee ou les elements sont des liste (de meme dimension) ou des entiers et w=[n_1,..,n_s] une liste d'entiers donnee alors la fonction retourne v x w = [[l_1,n_1],[l_1,n_2],...,[l_r,n_s]] si les l_i sont des entiers ou v x w =[l_1 append n_1,...,l_r append n_s] si les l_i sont des listes. */ cart_prod_sequence:=function(v,w) n:=#v; m:=#w; S:={}; for i in [1..n] do for j in [1..m] do if Type(v[i]) ne SeqEnum then S:={[v[i],w[j]]} join S; else S:={Append(v[i],w[j])} join S; end if; end for; end for; return SetToSequence(S); end function; /* Soit v=[n1,..,ns] une liste d'entiers donnes alors cette fonction applique le programme cart_prod_sequence aux listes suivantes: [1..n1] x [1..n2] x ... x [1..ns] */ sequence_tuple:=function(v) w:=[]; n:=#v; for i in [1..n] do w[i]:=[1..v[i]]; end for; S:=w[1]; if n eq 1 then return S; else for i in [2..n] do S:=cart_prod_sequence(S,w[i]); end for; end if; return S; end function; /* Soit une table donnee "table" de dimemnsion n_1 x n_2 x ... x n_r et une liste donnee [i_1,i_2,...,i_r] avec i_k =< n_k alors evalue la table pour l'entree [i_1,i_2,...,i_r]. */ eval_table:=function(tab,v) x:=tab; n:=#v; for i in [1..n] do x:=x[v[i]]; end for; return x; end function; /* Cette fonction retourne les valeurs admissibles de discriminants */ discriminant:=function(p,f,N_0,N) v:=PrimeDivisors(N_0); n:=#v; for i in [1..N] do epsilon:=0; if ( GCD(i,f) eq 1 ) and ( Sprint(IsFundamentalDiscriminant(i)) eq "true" ) and ( KroneckerSymbol(i,p) eq -1 ) then for j in [1..n] do if KroneckerSymbol(i,v[j]) eq 1 then epsilon:=epsilon+1; end if; end for; if epsilon eq n then print(i); end if; end if; end for; return "termine"; end function; discriminant2:=function(p,f,N_0,N) v:=PrimeDivisors(N_0); n:=#v; list:=[]; for i in [1..N] do epsilon:=0; if ( GCD(i,f) eq 1 ) and ( Sprint(IsFundamentalDiscriminant(i)) eq "true" ) and ( KroneckerSymbol(i,p) eq -1 ) then for j in [1..n] do if KroneckerSymbol(i,v[j]) eq 1 then epsilon:=epsilon+1; end if; end for; if epsilon eq n then list[i]:=i; else list[i]:=0; end if; else list[i]:=0; end if; end for; return list; end function; /*------------------------------------------------------------------------------------------------------------------------------------------*/ /* Ce programme genere une liste de nombres de Bernoulli */ Bernoulli_liste:=function(p,precision) N:=p*(precision+Floor(Log(precision))+1); liste:=table([N]); /* i-eme entree est le i-eme nombre de Bernoulli */ for i in [1..N] do liste[i][1]:=BernoulliNumber(i); printf "le %o eme nombre de Bernoulli sur %o est calcule\n",i,N; end for; PrintFile("Bernoulli_liste.txt","liste:="); PrintFile("Bernoulli_liste.txt",liste); PrintFile("Bernoulli_liste.txt",";"); return liste; end function; /* Ce programme construit un polynome Bernoulli a partir de la liste precedente */ BernoulliPolynomial2:=function(n,liste) P:=PolynomialRing(Rationals()); sum:=x^n; for i in [1..n] do sum:=sum+Binomial(n,n-i)*liste[i][1]*x^(n-i); end for; return sum; end function; /*-----------------------------------------------------------------------------------------------------------------------------------------*/ /* NOUS ALLONS MAINTENANT DEFINIR DIFFERENTES FONCTIONS AFIN DE CALCULER DIVERS QUANTITES ASSOCIEES AUX INTEGRALES QUE NOUS VOULONS CALCULES */ /* Cette fonction prend en entree un ideal I d'un corps quadratique reel K=Q(\sqrt{D}) copremier a f ou f un entier >0 (le conducteur). Cette fonction verifie si il existe un generateur lambda tel que (lambda)=I , lamba>>0 et lambda = 1 mod 3, le corps de nombre K doit etre donne sous la forme K:=QuadraticField(D) */ is_principal:=function(K,w_t,I,f) Real:=RealField(100); O:=MaximalOrder(K); D:=Discriminant(K); epsilon:=FundamentalUnit(QuadraticField(D)); if Type(K) eq FldNum then epsilon:=O!(epsilon[1]+epsilon[2]*(2*w_t-1) ); end if; R,pi:=quo; N:=#R; true_false,gen:=IsPrincipal(I); if Type(K) eq FldNum then gen:=O!( gen[1]+gen[2]*w_t ); end if; if Sprint(true_false) eq "false" then printf "ideal not principal\n"; return true_false; end if; DD:=Real!D; if Type(K) eq FldNum then w_D:=Real!(SquareRoot(DD)+1)/2; else w_D:=Real!SquareRoot(DD); end if; for j in {1,-1} do for i in [1..N] do if pi(j*epsilon^i*O!gen) eq pi(1) and ( j*( epsilon[1]+epsilon[2]*w_D )^i*(gen[1]+gen[2]*w_D) ) gt 0 and Norm(j*epsilon^i*O!gen) gt 0 then return true_false,(j*epsilon^i*gen)[1]+(j*epsilon^i*gen)[2]*w_t; end if; end for; end for; printf "ideal principal but not generated by a totally positive element congruent to 1 modulo f\n"; return "false"; end function; /* Soit tau un element de K:=Q(sqrt{D})\bs Q alors cette fonction calcule End_K(Z*tau+Z) qui est un cerain ordre dans O_K */ endomorphism_order:=function(K,tau) D:=Discriminant(K); L:=QuadraticField(D); if Type(K) eq FldQuad and tau in K then tau2:=tau[1]+tau[2]*t; end if; if Type(K) eq FldNum and tau in K then tau2:=tau[1]+tau[2]*(t+1)/2; end if; poly:=MinimalPolynomial(tau2,Rationals()); /* poly=x^2+a/d1*x+b/d2 ou a,b in Z */ coeff:=Coefficients(poly); d1:=Denominator(coeff[2]); d2:=Denominator(coeff[1]); n:=Lcm(d1,d2); conducteur:=abs(n*tau2[2]); if (D mod 4) eq 1 then conducteur:=2*conducteur; end if; return "Ordre de conducteur",conducteur; end function; /* Ceci calcule une matrice dans SL_2(Z) associee a la cusp a/c */ matrice_cusp:=function(a,c) g,d,b:=ExtendedGreatestCommonDivisor(a,-c); return Matrix(Integers(),2,2,[a,b,c,d]); end function; /* Polynome binomial x(x-1)...(x-(n-1))/n! */ poly_binomial:=function(n) P:=PolynomialRing(Rationals()); t:=1; for i in [0..n-1] do t:=t*(x-i); end for; return t/Factorial(n); end function; /* function indicatrice sur les entiers (Z) */ Indicator:=function(r) if r in Integers() then return 1; else return 0; end if; end function; /* function "partie rationelle" pour Q */ Ratp:=function(r) return r-Floor(r); end function; /* Polynomes de Bernoulli */ Bernoulli:=function(n,x) return Evaluate(BernoulliPolynomial(n),x); end function; /* Polynomes periodiques de Bernoulli */ pBernoulli:=function(n,x) if n eq 1 then return (x-Floor(x))-1/2+1/2*Indicator(x); else return Evaluate(BernoulliPolynomial(n),x-Floor(x)); end if; end function; /* Sommes de Dedekind D_{(s,t)}^{r mod f}(a,c) L'entier c doit etre divisible par f et >=0 */ DedSum:=function(s,t,r,f,a,c) R,pi:=quo; /* R correspond a Z/fZ et pi:Z -> R */ c:=Integers()!c; /* "c" est donne parfois comme un quotient de deux entiers */ S:={h: h in [1..c] | pi(h) eq pi(r) }; sum:=c^(s-1)/(s*t)*(&+[ pBernoulli(s,h/c)*pBernoulli(t,h*a/c) : h in S ]); return sum; end function; /* Sommes de Dedekind D_{(s,t)}^{r mod f}(a,c) calculees dans Qp avec une certaine precision. L'ENTIER c doit etre divisible par f et >=0 APRES QUELQUES ESSEAIS JE ME SUIS RENDU COMPTE QUE CE PROGRAMME N'EST PAS PLUS RAPIDE QUE DedSum. */ DedSum_padic:=function(s,t,r,f,a,c,p,precision) Qp:=pAdicField(p,precision); R,pi:=quo; /* R correspond a Z/fZ et pi:Z -> R */ c:=Integers()!c; /* "c" est donne parfois comme un quotient de deux entiers */ S:={h: h in [1..c] | pi(h) eq pi(r) }; ff:=mapQp|x:-> ( Qp!pBernoulli(s,x/c)*Qp!pBernoulli(t,x*a/c) )>; sum:=c^(s-1)/(s*t)*(&+[ff(h): h in S ]); return sum; end function; /* Sommes de Dedkind ponderees par le vecteur v={n(d_0,r)}_{d_0,r} */ DedSum3:=function(s,t,j,f,a,c,l,S_1,S_2,delta) sum:=0; for r in S_1 do for d_0 in S_2 do if delta() ne 0 then sum:=sum+delta()*(d_0^(-l))*DedSum(s,t,r*j,f,a,c/d_0); end if; end for; end for; return sum; end function; /* Construction d'une fonction "Delta1":Divisors(N_0) x (Z/f)^* --> Z a partir d'une liste "v" pour un r donne dans (Z/fZ)^{\times}. La liste "v" doit etre donnee sous la forme v=[ n_{d_0} ]_{d_0|N_0} ou les coordonnes sont ordonnees en ordre croissant selon les diviseur de N_0. On suppose que sum_{d_0|N_0} d_0*n_d_0=0. L'element x est une liste de longeur 2: x=[d_0,r]\in Divisors(N_0)x(\ZZ/f)^{\times}. Pour connaitre la composante en [d_0,r] du vecteur Delta on doit faire: Delta1(p,f,N_0,r,[d_0,r],v) Le vecteur Delta est supporte sur les paires [d_0,r*p^n] ou n parcourt les entiers. */ Delta1:=function(p,f,N_0,r,x,v) S_2:=SequenceToSet(Divisors(N_0)); S_22:=Divisors(N_0); n:=#S_2; pi:=homIntegers(f)|>; m:=Modorder(p,f); residues:={}; for i in [1..m] do residues:=(residues join {pi(p^i*r)}); end for; w:=[]; for i in [1..n] do w[S_22[i]]:=i; end for; if pi(x[2]) in residues then return v[w[x[1]]]; else return 0; end if; end function; /* Meme chose que le programme precedent mais version symmetrique */ Delta1_sym:=function(p,f,N_0,r,x,v) S_2:=SequenceToSet(Divisors(N_0)); S_22:=Divisors(N_0); n:=#S_2; pi:=homIntegers(f)|>; m:=Modorder(p,f); residues:={}; for sign in {1,-1} do for i in [1..m] do residues:=(residues join {pi(sign*p^i*r)}); end for; end for; w:=[]; for i in [1..n] do w[S_22[i]]:=i; end for; if pi(x[2]) in residues then return v[w[x[1]]]; else return 0; end if; end function; /* Creation d'une fonction allant de: Divisors(N_0)x (Z/f)^* ------> Integers() a partir du vecteur "Delta1" */ Delta2:=function(p,f,N_0,r,v) S_1:={h: h in [1..f]| GCD(h,f) eq 1}; S_2:=SequenceToSet(Divisors(N_0)); S_3:=CartesianProduct(S_2,S_1); delta:=mapIntegers()|x:->Delta1(p,f,N_0,r,x,v)>; return S_2,S_1,delta; end function; /* Meme chose que le programme precedent mais version symmetrique */ Delta2_sym:=function(p,f,N_0,r,v) S_1:={h: h in [1..f]| GCD(h,f) eq 1}; S_2:=SequenceToSet(Divisors(N_0)); S_3:=CartesianProduct(S_2,S_1); delta:=mapIntegers()|x:->Delta1_sym(p,f,N_0,r,x,v)>; return S_2,S_1,delta; end function; /* Sommes de Dedkind ponderees par la liste v=[ n(d_0,r) ]_{d_0,r}, ce programme calcule: \sum_{d_0,r} n(d_0,r)D_{s,t}^{rj mod f}(a,c/d_0), l'entier c doit etre divisible par f*N_0 et >=0. */ DedSum_delta:=function(s,t,j,f,a,c,l,S_1,S_2,delta) sum:=0; for d_0 in S_2 do for r in S_1 do if delta() ne 0 then sum:=sum+delta()*(d_0^(-l))*DedSum(s,t,r*j,f,a,c/d_0); end if; end for; end for; return sum; end function; /* Calcule la valeur de la mesure \mu_j{i\infty\rightarrow a/c} (cette mesure depend de "delta") evaluee sur la boule "rectangulaire": (u+p^sZ)x(v+p^sZ). */ measure_ball:=function(p,u,v,s,delta,f,N_0,a,c,j) c:=Integers()!c; /* "c" might be given as a quotient of two integers */ phi:=homIntegers(f)|>; t:=Valuation(c,p); A:=Integers()!phi(p^(-s)); if Gcd(A,f) ne 1 then A:=A+p^s; else A:=A; end if; S_1:={h: h in [1..f]| GCD(h,f) eq 1}; S_2:=SequenceToSet(Divisors(N_0)); sum:=0; for r in S_1 do for d_0 in S_2 do d:=Integers()!(c/(f*d_0)); for h in [1..d] do B:=f*Ratp(j*r/f)*A*p^s+f*v; sum:=sum + delta()*pBernoulli(1, a/( c/(f*d_0) )*( h+B/(f*p^s) ) - (d_0*f*u)/p^s )*pBernoulli(1, 1/( c/(f*d_0) )*( h+B/(f*p^s) ) ); end for; end for; end for; return -12*sum; end function; /* logarithme p-adic pour des elements de Z,Q,Zp ou Qp */ log_padic:=function(p,n,x) Qp:=pAdicField(p,n); P:=PolynomialRing(Qp); f:=P!(z^(p-1)-1); h:=homIntegers(p)|>; u:=Rationals()!x; s:=Valuation(u,p); y:=Qp!(x/p^s); w:=Rationals()!y; v:=[]; for i in [1..p-1] do v[i]:=HenselLift(f,Qp!i,n); end for; r:=Integers()!h(w); y:=y/v[r]; return(Log(y)); end function; /* logarithme p-adic pour des elements de Kp=Qp(tp) ou D est un non residu modulo p et tp est tel que tp^2=D. */ log_padic2:=function(Kp,tp,a) n:=Precision(Kp); p:=Prime(Kp); P:=PolynomialRing(Kp); b:=Kp!a; f:=y^(p^2-1)-1; v:=ZeroMatrix(Kp,p,p); for i in [1..p] do for j in [1..p] do if i ne p or j ne p then v[i,j]:=HenselLift(f,Kp!(i+j*tp),n); end if; end for; end for; s_1:=Valuation(Rationals()!Coefficients(b)[1],p); s_2:=Valuation(Rationals()!Coefficients(b)[2],p); s:=Minimum(s_1,s_2); b:=b/p^s; r1:=Integers()!Coefficients(b)[1] mod p; r2:=Integers()!Coefficients(b)[2] mod p; if r1 eq 0 then r1:=p; end if; if r2 eq 0 then r2:=p; end if; b:=b/v[r1,r2]; return Log(b); end function; /* Cette fonction retourne une racine de l'unite w primitive d'ordre p^2-1 et un entier m tel que x/w^m \in 1+p w ne depend pas de x. */ log_beta:=function(Kp,tp,x) P:=PolynomialRing(Kp); R:=RingOfIntegers(Kp); F,h:=ResidueClassField(R); p:=Characteristic(F); n:=Precision(Kp); w:=PrimitiveElement(F); w:=Kp!Eltseq(w)[1]+(Kp!Eltseq(w)[2])*tp; w:=HenselLift(z^(p^2-1)-1,w,n); m:=Log(h(x)); return w,m; end function; /* Valuation p-adique Retourne la valuation p-adique d'un element x dans Kp */ valuation_padic:=function(Kp,x) a:=Kp!x; p:=Prime(Kp); r1:=Rationals()!Coefficients(a)[1]; r2:=Rationals()!Coefficients(a)[2]; s:=Minimum(Valuation(r1,p),Valuation(r2,p)); return s; end function; /* Approximation polynomiale de log_p(x) sur Zp^{\times} par un polynome f de degree p*(n+Floor(log(n))+1) avec coefficients dans (Qp,n) ou "n" est la precision, on a donc pour tout x dans Zp^* que |f(x)-log_p(x)|_p<1/p^n */ approx_log_padic:=function(p,n) m:=n+Floor(Log(n))+1; Qp:=pAdicField(p,n); P:=PolynomialRing(Qp); R:=PowerSeriesRing(Qp,m); u:=[]; v:=[]; w:=[]; for i in [1..p-1] do s:=1; for j in [1..p-1] do if j ne i then s:=s*(x-j); else s:=s; end if; end for; v[i]:=s^n; end for; f:=&+[ ((-1)^(j+1))*(y^j)/j: j in [1..m] ]; for i in [1..p-1] do w[i]:=( Integers()!log_padic(p,n,i)+Evaluate(f, y/i ) )/Evaluate(v[i],y+i) ; end for; s:=0; for i in [1..p-1] do u[i]:=Evaluate(w[i],x)*Evaluate(v[i],x+i); s:=s+Evaluate(u[i],x-i); end for; return s; end function; /* Cette fontion calcule certaines distrubution F_{k,r} sur Z=Z/e x Zp faisant intervenir le polynome de Bernoulli B_k et ou l'entier r appartient a (Z/f)^*. L'entier e doit etre copremier a p. Cette fonction retourne la valeur pour l'ouvert a+ep^n*Z. */ distribution:=function(p,N_0,k,r,delta,a,e,n) S_2:=SequenceToSet(Divisors(N_0)); sum:=0; for d_0 in S_2 do sum:=sum+delta()*(e*p^n/d_0)^(k-1)*pBernoulli(k,a*d_0/(e*p^n))/k; end for; return sum; end function; /* Cette fontion calcule certaines distrubution F_{k,r} sur Z=Z/e x Zp faisant intervenir le polynome de Bernoulli B_k et ou l'entier r appartient a (Z/f)^*. L'entier e doit etre copremier a p. Cette fonction retourne la valeur pour l'ouvert a+ep^n*Z. Ce programme prend entree un vecteur de polynomes de Bernoulli */ distribution2:=function(p,N_0,k,r,delta,a,e,n,Poly_Bernoulli) S_2:=SequenceToSet(Divisors(N_0)); sum:=0; for d_0 in S_2 do sum:=sum+delta()*(e*p^n/d_0)^(k-1)*Evaluate(Poly_Bernoulli[k],Ratp(a*d_0/(e*p^n) ) )/k; end for; return sum; end function; /* Calcul donnant la periode de l'unite modulaire beta_{j*\delta} de niveau f*N_0 pour un j dans (Z/f)^* construite a partie du vecteur delta: Donne la valeur de l'integrale suivante: 1/(2\pi i)*\int_{i\infty}^{M i\infty} dlog \beta_{\delta,j}(\tau) Prend comme entree un matrice M dans Gamma_0(fN_0). */ period:=function(p,f,N_0,j,delta,M) pi:=homIntegers(p)|>; S_1:={h: h in [1..f]| GCD(h,f) eq 1}; S_2:=SequenceToSet(Divisors(N_0)); a:=Integers()!M[1,1]; c:=Integers()!M[2,1]; D:=GCD(a,c); a:=Integers()!(a/D); c:=Integers()!(c/D); a:=Sign(c)*a; c:=Sign(c)*c; sum:=0; for d_0 in S_2 do for r in S_1 do if delta() ne 0 then sum:=sum+delta()*( DedSum(1,1,j*r,f,a,c/d_0)); end if; end for; end for; return -12*sum; end function; /* Calcul donnant la periode de l'unite modulaire beta_{j*\delta,p} de niveau p*f*N_0 pour un j dans (Z/f)^* construite a partie du vecteur delta: Donne la valeur de l'integrale suivante: 1/(2\pi i)*\int_{i\infty}^{M i\infty} dlog \beta_{j*\delta,p}(\tau) Prend comme entree un matrice M dans Gamma_0(fN_0). Remarquons que cette entier correspond a \nu_j{i\infty\rightarrow M(i\infty)}(Z_p) ou \nu_j{i\infty\rightarrow M(i\infty)} est la mesure pushforward de \mu_j{i\infty\rightarrow M(i\infty)} par rapport a la projection naturelle pi: (Zp x Zp)' ------> P^1(Qp) */ period_p:=function(p,f,N_0,j,delta,M) pi:=homIntegers(p)|>; S_1:={h: h in [1..f]| GCD(h,f) eq 1}; S_2:=SequenceToSet(Divisors(N_0)); a:=Integers()!M[1,1]; c:=Integers()!M[2,1]; D:=GCD(a,c); a:=Integers()!(a/D); c:=Integers()!(c/D); a:=Sign(c)*a; c:=Sign(c)*c; sum:=0; for d_0 in S_2 do for r in S_1 do sum:=sum+delta()*( DedSum(1,1,j*r,f,a,c/d_0)-DedSum(1,1,j*r,f,p*a,c/d_0) ); end for; end for; return -12*sum; end function; /* On a une mesure \nu_j{i\infty\rightarrow A/C} sur P^1(Qp) correspondant au pushforward de \mu_j{i\infty\rightarrow A/C}. Pour une matrice gamma dans Gamma_0(fN_0 Z[1/p]) cette fonction retourne la valeur: \nu_j{i\infty\rightarrow A/C}(gamma*Zp). */ measure_nu:=function(p,f,N_0,A,C,j,delta,gamma) a:=Integers()!gamma[1,1]; b:=Integers()!gamma[1,2]; c:=Integers()!gamma[2,1]; d:=Integers()!gamma[2,2]; alpha:=d*A-b*C; beta:=-c*A+a*C; q:=alpha/beta; alpha:=Numerator(q); beta:=Denominator(q); M:=matrice_cusp(alpha,beta); sum1:=period_p(p,f,N_0,(a*j) mod f,delta,M); if c eq 0 then sum2:=0; else q:=-d/c; alpha:=Numerator(q); beta:=Denominator(q); M:=matrice_cusp(alpha,beta); sum2:=period_p(p,f,N_0,(a*j) mod f,delta,M); end if; return (sum1-sum2); end function; /* Fonction retournant les matrices associees a la relation de distribution en p. */ matrices_distribution:=function(p) v:=[]; for i in [0..p-1] do v[i+1]:=Matrix(Rationals(),2,2,[p,i,0,1]); end for; return v; end function; /* Fonction retournant l'inverse des matrices associees a la relation de distribution en p avec coefficients entiers */ inv_matrices_distribution:=function(p) v:=[]; for i in [0..p-1] do v[i+1]:=Matrix(Rationals(),2,2,[1,-i,0,p]); end for; return v; end function; /* Cette function calcule l'intgerale suivant: \int_{h+ep^sZ} x^-t dF_{1,j}(x) pour un entier t >= 0 et un element j dans (Z/f)^*. L'entier e doit etre copremier a p et divisible par f*N_0. Pour que cet integrale ait un sens on suppose que s>=1 et que (h,p)=1. Cette algorithme utilise la serie de Taylor de x^-t centree en h. Notons que la serie de Taylor converge car s>=1. On utilise la formule suivante: \int_{U} x_p^(k-1)dF_{1,j}(x)=F_{k,j}(U) ou U est un compact-open dans Z_p x Z/e. */ negative_moment:=function(p,h,e,s,j,t,delta,precision) Qp:=pAdicField(p,precision); P:=PolynomialRing(Qp); N_0:=Maximum(Domain(delta)[1]); poly:=h^(-t)+h^(-t)*&+[Evaluate(poly_binomial(i),-t)*((x-h)/h)^i:i in [1..precision] ]; coeff:=Coefficients(poly); d:=#coeff; sum:=0; for i in [1..d] do sum:=sum+coeff[i]*distribution(p, N_0, i, j, delta, h, e, s); /* le degree du monome correspondant a coeff[i] est i-1 */ end for; return sum; end function; /* Cette function store sous la forme d'une matrice tous les moments negatifs de F_{1,j} sur les boules de la forme h+ep^sZ pour tous les h mod ep^s tels que (h,p)=1. Afin d'avoir convergence on suppose que s>=1. L'entier e doit etre divisible par f*N_0. Le programme cree une table de dimension "N+1" par "ep^s" et pour l'entree (m,h): on a la quantite \int_{h+ep^s} x^-m dF_{1,j}(x), On utilise le fait que: \int_U x_p^(m-1)dF_{1,j}(x)=F_{m,j}(U) pour n'importe quel compact-open dans Z_p x Z/e. */ negative_moments_matrix:=function(p,e,s,j,N,delta,precision) Qp:=pAdicField(p,precision); N_0:=Maximum(Domain(delta)[1]); P:=PolynomialRing(Qp); S:={h: h in [1..e*p^s]|GCD(h,p) eq 1}; M:=table([N+1]); for h in [1..e*p^s] do if h in S then M[1][h]:=distribution(p, N_0, 1, j, delta, h, e, s); else M[1][h]:=0; end if; end for; for i in [1..N] do poly1:=&+[Evaluate(poly_binomial(m),-i)*(x)^m:m in [1..precision] ]; for h in [1..e*p^s] do if h in S then poly2:=h^(-i)+h^(-i)*Evaluate(poly1,(x-h)/h); coeff:=Coefficients(poly2); sum:=0; for k in [1..#coeff] do sum:=sum + coeff[k]*distribution(p, N_0, k, j, delta, h, e, s); /* le degree du monome de coeff[k] est k-1 */ end for; M[i+1][h]:=sum; else M[i+1][h]:=0; end if; end for; end for; return M; end function; /* Cette function store sous la forme d'une matrice tous les moments negatifs de F_{1,j} sur les boules de la formes h+eqp^sZ pour les h copremiers a p. Si s=0 alors on pose q=p dans le cas contraire on pose q=1. L'entier e doit etre divisible par f*N_0 et copremier a p. Le programme cree une table "N+1" par "eqp^s" et pour l'entree (m,h): on a la quantite \int_{h+eqp^s} x^-m dF_{1,j}(x) ON NE SUPPOSE PLUS QUE s>=1. */ negative_moments_matrix2:=function(p,e,s,j,N,delta,precision,P,T) Qp:=pAdicField(p,precision); N_0:=Maximum(Domain(delta)[1]); Poly_Ring:=PolynomialRing(Qp); if s eq 0 then q:=p; epsilon:=1; else q:=1; epsilon:=0; end if; S:={h: h in [1..e*q*p^s]|GCD(h,p) eq 1}; M:=table([N+1]); for h in [1..e*q*p^s] do if h in S then M[1][h]:=distribution2(p, N_0, 1, j, delta, h, e, s+epsilon,P); else M[1][h]:=0; end if; end for; for i in [1..N] do poly1:=T[i]; for h in [1..e*q*p^s] do if h in S then poly2:=h^(-i)+h^(-i)*Evaluate(poly1,(x-h)/h); coeff:=Coefficients(poly2); sum:=0; for k in [1..#coeff] do sum:=sum+coeff[k]*distribution2(p, N_0, k, j, delta, h, e, s+epsilon,P); /* le degree du monome de coeff[k] est k-1 */ end for; M[i+1][h]:=sum; else M[i+1][h]:=0; end if; end for; end for; return M; end function; /* Cette function store sous la forme d'une matrice toutes les expressions de la forme lim_{u-->\infty,g=(p-1)*p^u} F_{g-l+1,j}(h+e*p^s) pour 1=< h =< ep^s et les entiers 0=< l =< N et j dans (Z/f)^* L'entier e doit etre divisible par f*N_0 et copremier a p. Le programme cree une matrice "N+1" par "ep^s" et pour l'entree (m+1,h) on a: lim_{u-->\infty,g=(p-1)*p^u} F_{g-m+1,j}(h+e*p^s). Remarquons que lorsque s>=1 et h divisible par p alors cette limite des egale a 0. De plus, si s>=1 et h est copremier a p alors: lim_{u-->\infty,g=(p-1)*p^u} F_{g-m+1,j}(h+e*p^s)=\int_{h+ep^sZ} x_p^(-m) F_{1,j}(x). */ negative_moments_matrix3:=function(p,e,s,j,N,delta,precision,P,T) Qp:=pAdicField(p,precision); M:=negative_moments_matrix2(p,e,s,j,N,delta,precision,P,T); S1:={h: h in [1..e]}; A:=table([N+1]); if s eq 0 and ( GCD(e,p) eq 1 ) then for i in [1..N+1] do for k in S1 do S2:={h: h in [1..e*p]|GCD(h,p) eq 1 and (h mod e) eq k }; sum:=0; for h in S2 do sum:=sum+M[i][h]; end for; A[i][k]:=sum; end for; end for; return A; else return M; end if; end function; /* Cette function retourne une liste de f elements ou le j ieme element correspond a la matrice obtenue a partir du programme: negative_moments_matrix3(p,e,s,j,N,delta,precision). */ negative_moments_matrix4:=function(p,f,e,s,N,delta,precision,P,T) S:={h:h in [1..f]|GCD(h,f) eq 1}; V:=[]; for j in S do V[j]:=negative_moments_matrix3(p,e,s,j,N,delta,precision,P,T); end for; return V; end function; /* Ce programme donne les tables du programme negative_moments_matrix4 pour les generateurs de Gamma0(12): gen[1]:=[1,1,0,1]; gen[2]:=[5,-1,36,-7]; gen[3]:=[5,-4,24,-19]; gen[4]:=[7,-5,24,-17]; gen[5]:=[5,-3,12,-7]; Les entiers possibles pour e sont: {36,24,12} join {p*36,p*24,p*12}; On a en sorties les tables: negative_moments_matrix4(p,3,e,1,N,delta,precision) pour les valeurs de e mentionnees precedemment. */ negative_moments_gamma0_12:=function(p,delta,N,precision,liste) Qp:=pAdicField(p,precision); Poly_Ring:=PolynomialRing(Qp); P:=[]; T:=[]; G:=Gamma0(12); gen1:=Generators(G); r:=#gen1; gen2:=[]; for i in [1..r] do gen2[i]:=Matrix(Integers(),2,2,Eltseq(gen1[i])); end for; printf "generate the list of periodic Bernoulli polynomials\n"; for i in [1..precision+1] do P[i]:=BernoulliPolynomial2(i,liste); end for; printf "generate the list of Taylor series of (1-x)^(-i)\n"; for i in [1..N] do T[i]:=&+[Evaluate(poly_binomial(m),-i)*(x)^m:m in [1..precision]]; end for; W2_1:=negative_moments_matrix4(p,3,36,1,N,delta,precision,P,T); printf "1/6 des moments calcules\n"; W2_0:=negative_moments_matrix4(p,3,36,0,N,delta,precision,P,T); printf "2/6 des moments calcules\n"; W3_1:=negative_moments_matrix4(p,3,24,1,N,delta,precision,P,T); printf "3/6 des moments calcules\n"; W3_0:=negative_moments_matrix4(p,3,24,0,N,delta,precision,P,T); printf "4/6 des moments calcules\n"; W5_1:=negative_moments_matrix4(p,3,12,1,N,delta,precision,P,T); printf "5/6 des moments calcules\n"; W5_0:=negative_moments_matrix4(p,3,12,0,N,delta,precision,P,T); /* On transforme tout en nombres rationels */ T2_0:=table([2]); V2_0:=table([N+1]); for j in [1..2] do for h in [1..36] do for m in [1..N+1] do V2_0[m][h]:=Rationals()!W2_0[j][m][h] ; end for; end for; T2_0[j][1]:=V2_0; end for; T2_1:=table([2]); V2_1:=table([N+1]); for j in [1..2] do for h in [1..p*36] do for m in [1..N+1] do V2_1[m][h]:=Rationals()!W2_1[j][m][h]; end for; end for; T2_1[j][1]:=V2_1; end for; T3_0:=table([2]); V3_0:=table([N+1]); for j in [1..2] do for h in [1..24] do for m in [1..N+1] do V3_0[m][h]:=Rationals()!W3_0[j][m][h]; end for; end for; T3_0[j][1]:=V3_0; end for; T3_1:=table([2]); V3_1:=table([N+1]); for j in [1..2] do for h in [1..p*24] do for m in [1..N+1] do V3_1[m][h]:=Rationals()!W3_1[j][m][h]; end for; end for; T3_1[j][1]:=V3_1; end for; T5_0:=table([2]); V5_0:=table([N+1]); for j in [1..2] do for h in [1..12] do for m in [1..N+1] do V5_0[m][h]:=Rationals()!W5_0[j][m][h]; end for; end for; T5_0[j][1]:=V5_0; end for; T5_1:=table([2]); V5_1:=table([N+1]); for j in [1..2] do for h in [1..p*12] do for m in [1..N+1] do V5_1[m][h]:=Rationals()!W5_1[j][m][h]; end for; end for; T5_1[j][1]:=V5_1; end for; PrintFile("negative_moments.txt","T2_1:="); PrintFile("negative_moments.txt",T2_1); PrintFile("negative_moments.txt",";"); PrintFile("negative_moments.txt","T2_0:="); PrintFile("negative_moments.txt",T2_0); PrintFile("negative_moments.txt",";"); PrintFile("negative_moments.txt","T3_1:="); PrintFile("negative_moments.txt",T3_1); PrintFile("negative_moments.txt",";"); PrintFile("negative_moments.txt","T3_0:="); PrintFile("negative_moments.txt",T3_0); PrintFile("negative_moments.txt",";"); PrintFile("negative_moments.txt","T5_1:="); PrintFile("negative_moments.txt",T5_1); PrintFile("negative_moments.txt",";"); PrintFile("negative_moments.txt","T5_0:="); PrintFile("negative_moments.txt",T5_0); PrintFile("negative_moments.txt",";"); return T2_0,T2_1,T3_0,T3_1,T5_0,T5_1; end function; /* Ce programme donne les tables du programme negative_moments_matrix4 pour les Les entiers e suivants: {96,84,108,72} avec pour p ne 7: s=1 et pour p=7 on a s=1 pour {96,108,72} et s=2 et on remplace 84 par 12; On a en sorties les tables: negative_moments_matrix4(p,3,e,1 ou 2,N,delta,precision) pour les valeurs de e mentionnees precedemment. Ces moments interviendront lors du programme moments_inverse */ negative_moments_gamma0_12_inverse:=function(p,delta,N,precision,liste) G:=Gamma0(12); gen1:=Generators(G); r:=#gen1; gen2:=[]; Qp:=pAdicField(p,precision); Poly_Ring:=PolynomialRing(Qp); P:=[]; T:=[]; for i in [1..r] do gen2[i]:=Matrix(Integers(),2,2,Eltseq(gen1[i])); end for; /*--------------------------------------------------------------*/ printf "generate the list of periodic Bernoulli polynomials\n"; for i in [1..precision+1] do P[i]:=BernoulliPolynomial2(i,liste); end for; printf "generate the list of Taylor series of (1-x)^(-i)\n"; for i in [1..N] do T[i]:=&+[Evaluate(poly_binomial(m),-i)*(x)^m:m in [1..precision]]; end for; /*--------------------------------------------------------------*/ I2_1:=negative_moments_matrix4(p,3,96,1,N,delta,precision,P,T); printf "1/4 des moments calcules\n"; if p eq 7 then I3_1:=negative_moments_matrix4(p,3,12,2,N,delta,precision,P,T); else I3_1:=negative_moments_matrix4(p,3,84,1,N,delta,precision,P,T); end if; printf "2/4 des moments calcules\n"; I4_1:=negative_moments_matrix4(p,3,108,1,N,delta,precision,P,T); printf "3/4 des moments calcules\n"; I5_1:=negative_moments_matrix4(p,3,72,1,N,delta,precision,P,T); printf "4/4 des moments calcules\n"; /* On transforme tout en nombres rationels */ J2_1:=table([2]); V2_1:=table([N+1]); for j in [1..2] do for h in [1..p*96] do for m in [1..N+1] do V2_1[m][h]:=Rationals()!I2_1[j][m][h] ; end for; end for; J2_1[j][1]:=V2_1; end for; J3_1:=table([2]); V3_1:=table([N+1]); for j in [1..2] do for h in [1..p*84] do for m in [1..N+1] do V3_1[m][h]:=Rationals()!I3_1[j][m][h]; end for; end for; J3_1[j][1]:=V3_1; end for; J4_1:=table([2]); V4_1:=table([N+1]); for j in [1..2] do for h in [1..p*108] do for m in [1..N+1] do V4_1[m][h]:=Rationals()!I4_1[j][m][h]; end for; end for; J4_1[j][1]:=V4_1; end for; J5_1:=table([2]); V5_1:=table([N+1]); for j in [1..2] do for h in [1..p*72] do for m in [1..N+1] do V5_1[m][h]:=Rationals()!I5_1[j][m][h] ; end for; end for; J5_1[j][1]:=V5_1; end for; PrintFile("negative_moments_inverse.txt","J2_1:="); PrintFile("negative_moments_inverse.txt",J2_1); PrintFile("negative_moments_inverse.txt",";"); PrintFile("negative_moments_inverse.txt","J3_1:="); PrintFile("negative_moments_inverse.txt",J3_1); PrintFile("negative_moments_inverse.txt",";"); PrintFile("negative_moments_inverse.txt","J4_1:="); PrintFile("negative_moments_inverse.txt",J4_1); PrintFile("negative_moments_inverse.txt",";"); PrintFile("negative_moments_inverse.txt","J5_1:="); PrintFile("negative_moments_inverse.txt",J5_1); PrintFile("negative_moments_inverse.txt",";"); return J2_1,J3_1,J4_1,J5_1; end function; /* Cette fonction calcule tous les moments d'ordre plus petit ou egal a N de la mesure \nu_{j}{i\infty\rightarrow b/(ep^s)} sur Zp i.e. \int_{Zp} x^i d\nu_j{i\infty\rightarrow b/ep^s}(x) pour 0=< i=< N La fonction retourne un vecteur w ou w[i+1] correspond a: \int_{Zp} x^i d\nu_j{i\infty\rightarrow b/ep^s}(x) Ce programme prend en entree une matrice V obtenue a partir du programme " negative_moments_matrix4(p,f,e,s,N*,delta,precision) = V ". pour N* >= N. On a 6 choix pour V: T2_0,T2_1,T3_0,T3_1,T5_0,T5_1 correspondant aux denominateurs {36,24,12} join {36*p,24*p,12*p} on peut aussi avoir e=12 avec s=2 "Les matrices doivent etre donnees en nombre rationels" */ moments_nu:=function(p,f,b,e,s,j,delta,N,V,precision,Poly_Bernoulli) Qp:=pAdicField(p,precision); T:={h:h in [1..f]|GCD(h,f) eq 1}; W:=[]; for i in [0..N] do sum:=0; for l in [0..i] do for t in T do if s ne 0 then S:={h: h in [1..e*p^s]|GCD(h,p) eq 1 and (h mod f) eq ((t*j) mod f) }; else S:={h: h in [1..e]|(h mod f) eq ( (t*j) mod f) }; end if; for h in S do sum:=sum + Binomial(i,l)*(b/(e*p^s))^(i-l)*(-1)^l*Evaluate( Poly_Bernoulli[l+1], Ratp(h*b/(e*p^s)) )/(l+1)*V[t][l+1][h]; end for; end for; end for; W[i+1]:=-12*Qp!(Rationals()!sum); /* Contient le ieme moment */ end for; return W; end function; /* Coerce les vecteurs coordonnees de V1 de Qp a Z. */ change1:=function(V,f,Ring) v:=dimension(V); n:=#v; w:=Remove(v,n); U:=table(w); for i in {h: h in [1..f]| GCD(h,f) eq 1} do for j in [1..v[2]] do U[i][j]:=Ring!V[i][j]; end for; end for; return U; end function; /* Coerce les vecteurs coordonnees de V2 de Qp a Z. */ change2:=function(V,f,Ring) v:=dimension(V); n:=#v; w:=Remove(v,n); U:=table(w); for i in {h: h in [1..f]| GCD(h,f) eq 1} do for j in [1..v[2]] do U[i][j]:=Ring!V[i][j]; end for; end for; return U; end function; /* Coerce les vecteurs coordonnees de V3 de Qp a Z. */ change3:=function(V,f,Ring) v:=dimension(V); n:=#v; w:=Remove(v,n); U:=table(w); for i in [1..v[1]] do for j in {h: h in [1..f]| GCD(h,f) eq 1} do for k in [1..v[3]] do for l in [1..v[4]] do if i ne 1 then U[i][j][k][l]:=Ring!V[i][j][k][l]; else U[i][j][k][l]:=Ring!0; end if; end for; end for; end for; end for; return U; end function; /* Coerce les vecteurs coordonnees de V4 de Qp a Z. */ change4:=function(V,f,Ring) v:=dimension(V); print(v[3]); n:=#v; w:=Remove(v,n); U:=table([5,2]); for i in [1..v[1]] do for j in {h: h in [1..f]| GCD(h,f) eq 1} do for k in [1..v[3]] do if i ne 1 then U[i][j][k]:=Ring!V[i][j][k]; else U[i][j][k]:=Ring!0; end if; end for; end for; end for; return U; end function; change4_new:=function(V,f,Ring,v) print(v[3]); n:=#v; w:=Remove(v,n); U:=table(w); for i in [1..v[1]] do for j in {h: h in [1..f]| GCD(h,f) eq 1} do for k in [1..v[3]] do if (i ne 1) and (i ne 3) then U[i][j][k]:=Ring!V[i][j][k]; else U[i][j][k]:=Ring!0; end if; end for; end for; end for; return U; end function; /* Ce programme construit un vecteur de longueur 5 ou seules les entrees 2=< k =<5 sont definies. Une telle entree k contient la valeur moments_shifted2(p,f,N_0,delta,N,precision). On retourne un vecteur WW de longueur 5. */ moments_shifted3:=function(p,f,N_0,delta,N,precision,T2_0,T2_1,T3_0,T3_1,T5_0,T5_1,liste) Poly_Ring:=PolynomialRing(Rationals()); Poly:=[]; /*--------------------------------------------------------------*/ printf "generate the list of periodic Bernoulli polynomials\n"; for i in [1..precision+1] do Poly[i]:=Poly_Ring!BernoulliPolynomial2(i,liste); end for; /*--------------------------------------------------------------*/ /* Premiere Etape */ G:=Gamma0(f*N_0); gen1:=Generators(G); r:=#gen1; gen2:=[]; for i in [1..r] do gen2[i]:=Matrix(Integers(),2,2,Eltseq(gen1[i])); end for; Set:={}; for i in [2..r] do Set:=Set join { abs( gen2[i][2,1] ) }; end for; WW:=[]; /* On va tout storer dans une liste WW */ /* deuxieme etape */ for c in Set do /* 1er do */ for j in [2..r] do /* 2e do */ if (gen2[j][2,1]) eq c then /* 1er if */ aa:=Sign(gen2[j][2,1])*gen2[j][1,1]; cc:=Sign(gen2[j][2,1])*gen2[j][2,1]; A:=matrice_cusp(aa,cc); P:=inv_matrices_distribution(p); w:=[]; v2:=[]; /* Troisieme etape */ V:=[]; S:={h: h in [1..f]| GCD(h,f) eq 1}; for h in S do /* 3e do */ for i in [0..p-1] do /* 4e do */ B:=P[i+1]*A; q:=B[1,1]/B[2,1]; aa:=Numerator(q); cc:=Denominator(q); epsilon:=Sign(c); aa:=epsilon*aa; cc:=epsilon*cc; b:=aa; s:=Valuation(cc,p); e:=Integers()!(cc/p^s); if s eq 1 and e eq 36 then v1:=moments_nu(p,f,b,e,s,h,delta,N,T2_1,precision,Poly); end if; if s eq 0 and e eq 36 then v1:=moments_nu(p,f,b,e,s,h,delta,N,T2_0,precision,Poly); end if; if s eq 1 and e eq 24 then v1:=moments_nu(p,f,b,e,s,h,delta,N,T3_1,precision,Poly); end if; if s eq 0 and e eq 24 then v1:=moments_nu(p,f,b,e,s,h,delta,N,T3_0,precision,Poly); end if; if s eq 1 and e eq 12 then v1:=moments_nu(p,f,b,e,s,h,delta,N,T5_1,precision,Poly); end if; if s eq 0 and e eq 12 then v1:=moments_nu(p,f,b,e,s,h,delta,N,T5_0,precision,Poly); end if; for k in [1..N+1] do v2[k]:=v1[k]*p^(k-1); end for; w[i+1]:=v2; end for; /* 4e do */ V[h]:=w; end for; /* 3e do */ WW[j]:=V; end if; /* 1er if */ end for; /* 2e do */ end for; /* 1er do */ WWW:=change3(WW,f,Integers()); PrintFile("moments_3.txt","M3:="); PrintFile("moments_3.txt",WWW); PrintFile("moments_3.txt",";"); return WW; end function; /* Cette fonction retourne le developpement de Taylor centre en 1=Integers(p)|>; P:=PolynomialRing(Rationals()); R:=PowerSeriesRing(Rationals(),precision); if special_case eq 0 then j:=Integers()!h( 1/(f*N_0) ); else j:=Integers()!h( -1/(f*N_0) ); end if; if special_case eq 0 then poly:=( (y+j)/( -(f*N_0)*(y+j)+1 ) )^(-n); /* Serie de Taylor centree en 0 */ else poly:=( (y+j)/( (f*N_0)*(y+j)+1 ) )^(-n); /* Serie de Taylor centree en 0 */ end if; return Evaluate(poly,x); end function; /* Calcule le rayclass group rapidement */ ray_group:=function(D,f) K:=QuadraticField(D); O:=MaximalOrder(K); I:=ideal; return RayClassGroup(I,[1,2]); end function; ray_group2:=function(D,f) K:=QuadraticField(D); O:=MaximalOrder(K); I:=ideal; return RayClassGroup(I,[1,2]),RayClassGroup(I); end function; ray_field:=function(D,f) K:=QuadraticField(D); O:=MaximalOrder(K); I:=ideal; G,pi:=RayClassGroup(I,[1,2]); L:=RayClassField(pi); return L; end function; ray_group3:=function(p,f,N_0,N) list:=discriminant2(p,f,N_0,N); M:=#list; for i in [1..M] do if list[i] ne 0 then D:=list[i]; K:=QuadraticField(D); O:=MaximalOrder(K); I:=ideal; h:=#RayClassGroup(I,[1,2]); printf "D = %o, h = %o\n\n",D,h; end if; end for; return "fin"; end function; ray_field2:=function(D,f) K:=QuadraticField(D); O:=MaximalOrder(K); I:=ideal; G1,pi1:=RayClassGroup(I,[1,2]); G2,pi2:=RayClassGroup(I); L1:=RayClassField(pi1); L2:=RayClassField(pi2); return L1,L2; end function; /* ------------------------------------------------------------------------------------------------------------ */ /* Soit i = 1/(f*N_0) (mod p) avec 0=< i =< p-1 et c dans {96,84,108,72} et a tel que (a-i*c) n'est pas divisible par p. ATTENTION ATTENTION, lorsque p=7 on a que p|84!!! alors cette fonction calcule les integrales suivantes: \int_{i+pZp} (t-i)^n d\mu_j{i\infty\rightarrow a/c}(t) pour n=0..N et pour les j dans (Z/f)^* donnes. Elle retourne un vecteur de longueur N+1 contenant les moments calcules precedemment. On utilise l'identite: \int_{i+pZp} (t-i)^n d\nu_j{i\infty\rightarrow a/c}(t)= p^n \int_{Zp} u^n d\nu_j{i\infty\rightarrow (a-i*c)/p*c}(u) On peut aussi avoir e=12 avec s=0 */ /* ATTENTION ATTENTION ATTENTION ATTENTION ATTENTION */ /* Lorsque p = 7 on a que p|84 et ceci provoque un probleme de convergence. Dans ce cas au lieu de prendre i = 1/f*N_0 mod p on prend i = -1/f*N_0 mod p et alors 5/24 --> -5/36 avec p ne divisant pas 36 */ moments_shifted2_inverse:=function(p,f,a,c,i,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly) /* Premiere etape */ a:=Sign(c)*a; c:=Sign(c)*c; if c eq 96 then W1:=J2_1; end if; if c eq 84 then /* ATTENTION si p=7 alors p|84 */ W1:=J3_1; end if; if c eq 108 then W1:=J4_1; end if; if c eq 72 then W1:=J5_1; end if; if c eq 36 then /* In the case where p=7 avec le 3e generateur M{5}{24}{-4}{-19} pour la matrice M{1}{0}{-12}{1} */ W1:=T2_1; end if; if c eq 12 then W1:=T5_0; end if; /* Deuxieme etape */ A:=matrice_cusp(a,c); P:=inv_matrices_distribution(p); w:=[]; v2:=[]; /* Troisieme etape */ V:=[]; S:={h: h in [1..f]| GCD(h,f) eq 1}; for h in S do B:=P[i+1]*A; q:=B[1,1]/B[2,1]; a:=Numerator(q); c:=Denominator(q); epsilon:=Sign(c); a:=epsilon*a; c:=epsilon*c; b:=a; s:=Valuation(c,p); e:=Integers()!(c/p^s); v1:=moments_nu(p,f,b,e,s,h,delta,N,W1,precision,Poly); for k in [1..N+1] do v2[k]:=v1[k]*p^(k-1); end for; w[i+1]:=v2; V[h]:=w; end for; return V; end function; /* Cette fonction calcule les N+1 moments negatifs suivants: \int_{ P^1(Qp)-Zp } t^(-n)d\nu_j{i\infty\rightarrow a/c} (t) pour un certain j dans (Z/f)^* et 0=< n =< N pour c dans {36,24,12} On utilise le fait que: \int_{ P^1(Qp)-Zp } t^(-n)d\nu_j{i\infty\rightarrow a/c} (t) = \int_{ j+pZp } ( u/(-f*N_0*u+1) )^(-n) d\nu_j{i\infty\rightarrow a/c} ( u/(-f*N_0*u+1) ) en utilisant l'invariance de \nu on obtient: \int_{ j+pZp } ( u/(-f*N_0*u+1) )^(-n) d\nu_j{1/(f*N_0) \rightarrow a/( f*N_0*a+c ) }(u) Nous faisons correspondre la matrice J2_1 --> 96 J3_1 --> 84 J4_1 --> 108 J5_1 --> 72 Il faut faire tres attention car lorsque p=7 on a que p|84 (ATTENTION!! ATTENTION!! ATTENTION!! ATTENTION!!) On utilise dans ce cas la matrice M{1}{0}{-12}{1} */ inverse_moments:=function(p,f,N_0,a,c,j,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly) Qp:=pAdicField(p,precision); pi:=homIntegers(p)|>; alpha:=Integers()!pi( 1/(f*N_0) ); alpha2:=Integers()!pi( -1/(f*N_0) ); if (p eq 7) and (a eq 5) and (c eq 24) then a:=a; c:=-f*N_0*a+c; else a:=a; c:=(f*N_0*a+c); end if; a:=Sign(c)*a; c:=Sign(c)*c; V:=[]; /* On va storen les N+1 moments inverses dans V */ if (p eq 7) and (a eq -5) and (c eq 36) then M1:=moments_shifted2_inverse(p,f,a,c,alpha2,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly); M2:=moments_shifted2_inverse(p,f,-1,f*N_0,alpha2,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly); else M1:=moments_shifted2_inverse(p,f,a,c,alpha,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly); M2:=moments_shifted2_inverse(p,f,1,f*N_0,alpha,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly); end if; if (p eq 7) and (a eq -5) and (c eq 36) then special_case:=1; alpha:=alpha2; else special_case:=0; end if; for n in [0..N] do poly:=approx_moebius(p,f,N_0,n,precision,special_case); coeff:=Coefficients(poly); d:=#coeff; sum:=0; for k in [1..d] do sum:=sum+coeff[k]*Rationals()!( M1[j][alpha+1][k]-M2[j][alpha+1][k] ); end for; V[n+1]:=Qp!sum; end for; return V; end function; /* Ce programme applique le programme inverse moments pour chaque entier j dans (Z/f)^* et retourne une vecteur de longueur f ou la j ieme entree est inverse_moments:=function(p,f,N_0,a,c,j,delta,N,precision) */ inverse_moments2:=function(p,f,N_0,a,c,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly) H:={h: h in [1..f]| GCD(h,f) eq 1}; pi:=homIntegers(p)|>; alpha:=Integers()!pi( 1/(f*N_0) ); alpha2:=Integers()!pi( -1/(f*N_0) ); if (p eq 7) and (a eq 5) and (c eq 24) then a:=5; c:=-f*N_0*a+c; else a:=a; c:=(f*N_0*a+c); end if; a:=Sign(c)*a; c:=Sign(c)*c; V:=table([f]); /* On va storen les N+1 moments inverses dans V */ print(a); print(c); if (p eq 7) and (a eq -5) and (c eq 36) then M1:=moments_shifted2_inverse(p,f,a,c,alpha2,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly); M2:=moments_shifted2_inverse(p,f,-1,f*N_0,alpha2,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly); else M1:=moments_shifted2_inverse(p,f,a,c,alpha,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly); M2:=moments_shifted2_inverse(p,f,1,f*N_0,alpha,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly); end if; if (p eq 7) and (a eq -5) and (c eq 36) then special_case:=1; alpha:=alpha2; else special_case:=0; end if; for h in [1..f] do for n in [0..N] do if h in H then poly:=approx_moebius(p,f,N_0,n,precision,special_case); coeff:=Coefficients(poly); d:=#coeff; sum:=0; for k in [1..d] do sum:=sum+coeff[k]*( M1[h][alpha+1][k] - M2[h][alpha+1][k] ); end for; V[h][n+1]:=sum; else V[h][n+1]:=0; end if; end for; end for; return V; end function; /* Ce programme construit un vecteur de longueur 5 ou seules les entrees 2=< k =<5 sont definies. Une telle entree k contient la valeur inverse_moments2(p,f,N_0,delta,delta,N,precision). On retourne un vecteur WW de longueur 5. */ inverse_moments3:=function(p,f,N_0,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,liste) Poly:=[]; /*--------------------------------------------------------------*/ printf "generate the list of periodic Bernoulli polynomials\n"; for i in [1..precision+1] do Poly[i]:=BernoulliPolynomial2(i,liste); end for; /*--------------------------------------------------------------*/ G:=Gamma0(f*N_0); gen1:=Generators(G); r:=#gen1; gen2:=[]; for i in [1..r] do gen2[i]:=Matrix(Integers(),2,2,Eltseq(gen1[i])); end for; V:=[]; for i in [2..r] do a:=Sign(gen2[i][2,1])*gen2[i][1,1]; c:=Sign(gen2[i][2,1])*gen2[i][2,1]; print(a/c); V[i]:=inverse_moments2(p,f,N_0,a,c,delta,N,precision,J2_1,J3_1,J4_1,J5_1,T2_1,T5_0,Poly); end for; WWW:=change4(V,f,Integers()); PrintFile("moments_4.txt","M4:="); PrintFile("moments_4.txt",WWW); PrintFile("moments_4.txt",";"); return V; end function; /* table T de hauteur 4 de dimensions v:=[2,1,N+1,e*p] -> on retourne S de dimension [2,N+1,e*p] */ swap_Qp:=function(Qp,T) v:=dimension(T); S:=table([v[1],v[3]]); for i in [1..v[1]] do for j in [1..v[3]] do for k in [1..v[4]] do S[i][j][k]:=Qp!T[i][1][j][k]; end for; end for; end for; return S; end function; /* */ swap_all_Qp:=function(Qp,T2_0,T2_1,T3_0,T3_1,T5_0,T5_1,J2_1,J3_1,J4_1,J5_1) S2_0:=swap_Qp(Qp,T2_0); S2_1:=swap_Qp(Qp,T2_1); S3_0:=swap_Qp(Qp,T3_0); S3_1:=swap_Qp(Qp,T3_1); S5_0:=swap_Qp(Qp,T5_0); S5_1:=swap_Qp(Qp,T5_1); L2_1:=swap_Qp(Qp,J2_1); L3_1:=swap_Qp(Qp,J3_1); L4_1:=swap_Qp(Qp,J4_1); L5_1:=swap_Qp(Qp,J5_1); return S2_0,S2_1,S3_0,S3_1,S5_0,S5_1,L2_1,L3_1,L4_1,L5_1; end function; /* Ce programme compare les entrees de deux tables T1 et T2. On doit avoir dim(T1)=n1 x n2 x x n3 x n4 et dim(T2)=m1 x m2 x m3 x m4 avec ni =< mi. En particulier les tables T1 et T2 ont la meme hauteur i.e. 4 */ compare3:=function(T1,T2,p,precision) Qp:=pAdicField(p,precision); v1:=dimension(T1); v2:=dimension(T2); n1:=#v1; n2:=#v2; if n1 ne n2 then return "hauteur des tables incompatibles"; end if; for i in [1..n1] do if v1[i] gt v2[i] then return "dimensions incompatibles"; end if; end for; T3:=table([ v1[1],v1[2],v1[3] ]); for i in [ 2..v1[1] ] do for j in [1.. v1[2] ] do for k in [1..v1[3] ] do for l in [1.. v1[4]] do T3[i][j][k][l]:=Integers(p^precision)!( Rationals()!T1[i][j][k][l] - Rationals()!T2[i][j][k][l] ); end for; end for; end for; end for; return T3; end function; /* Ce programme compare les entrees de deux tables T1 et T2. On doit avoir dim(T1)=n1 x n2 x x n3 et dim(T2)=m1 x m2 x m3 avec ni =< mi. En particulier les tables T1 et T2 ont la meme hauteur i.e. 4 */ compare4:=function(T1,T2,p,precision) Qp:=pAdicField(p,precision); v1:=dimension(T1); v2:=dimension(T2); n1:=#v1; n2:=#v2; if n1 ne n2 then return "hauteur des tables incompatibles"; end if; for i in [1..n1] do if v1[i] gt v2[i] then return "dimensions incompatibles"; end if; end for; T3:=table([ v1[1],v1[2]]); for i in [ 2..v1[1] ] do for j in [1.. v1[2] ] do for k in [1..v1[3] ] do T3[i][j][k]:=Integers(p^precision)!( Rationals()!T1[i][j][k] - Rationals()!T2[i][j][k] ); end for; end for; end for; return T3; end function; /* -------------------------------------------------------------------------------------------------------------- */ /* Formules pour les moments sur l'espace Zp x Zp^(\times), le programme calcule l'integrale suivant pour m et n positifs: \int_X x^ny^m d\wt{\mu}_j(i\infty \rightarrow a/c)(x,y) */ moments_1:=function(p,f,N_0,j,a,c,n,m,delta) S_1:={h: h in [1..f]| GCD(h,f) eq 1}; S_2:=SequenceToSet(Divisors(N_0)); R,pi:=quo; /* R correspond a Z/fZ et pi:Z -> R */ sum:=0; for i in [0..n] do for d_0 in S_2 do for r in S_1 do if delta() ne 0 then S:={h: h in [1..Integers()!(c/d_0)] | pi(h) eq pi(r) }; for h in S do sum:=sum+Binomial(n,i)*(a/c)^(n-i)*( (-1)^i )*( 1/((n+m-i+1)*(i+1)) )*(c/d_0)^(n+m-i)*delta()*d_0^(-i) *pBernoulli(n+m-i+1,h/(c/d_0))*( pBernoulli(i+1,a*h/(c/d_0))-p^(n+m-i)*pBernoulli(i+1, p*a*h/(c/d_0)) ); end for; end if; end for; end for; end for; return -12*(1/f^(m+n))*sum; end function; /* Formules pour les moments sur l'espace Zp x Zp^(\times), le programme calcule l'integrale suivant pour m et n positifs: \int_X x^ny^m d\wt{\mu}_j(i\infty \rightarrow a/c)(x,y) */ moments_1_modify:=function(p,f,N_0,j,a,c,n,m,delta,P) S_1:={h: h in [1..f]| GCD(h,f) eq 1}; S_2:=SequenceToSet(Divisors(N_0)); R,pi:=quo; /* R correspond a Z/fZ et pi:Z -> R */ sum:=0; for i in [0..n] do for d_0 in S_2 do for r in S_1 do if delta() ne 0 then S:={h: h in [1..Integers()!(c/d_0)] | pi(h) eq pi(r*j) }; for h in S do sum:=sum+Binomial(n,i)*(a/c)^(n-i)*( (-1)^i )*( 1/((n+m-i+1)*(i+1)) )*(c/d_0)^(n+m-i)*delta()*d_0^(-i) *Evaluate(P[n+m-i+1], Ratp(h/(c/d_0)) )*( Evaluate( P[i+1], Ratp( a*h/(c/d_0) ) )-p^(n+m-i)*Evaluate(P[i+1],Ratp(p*a*h/(c/d_0)) ) ); end for; end if; end for; end for; end for; return -12*(1/f^(m+n))*sum; end function; /* L'integrale: \int_{Zp x Zp^*} log_p y dmu{i\infty\rightarrow a/c}(x,y) */ integral_1:=function(p,f,N_0,j,a,c,delta,precision,P) poly:=approx_log_padic(p,precision); coeff:=Coefficients(poly); d:=#coeff; sum:=0; printf "Le polynome est de degree %o\n",d; for i in [1..d] do sum:=sum+coeff[i]*moments_1_modify(p,f,N_0,j,a,c,0,i-1,delta,P); print(i); end for; return sum; end function; /* Formules pour les moments sur l'espace Zp^* x pZp: \int_X x^ny^m d\wt{\mu}_j(i\infty rightarrow a/c)(x,y) calcules dans Qp avec une certaine precision. */ moments_2:=function(p,f,N_0,j,a,c,n,m,delta,precision) S_1:={h: h in [1..f]| GCD(h,f) eq 1}; S_2:=SequenceToSet(Divisors(N_0)); R,pi:=quo; /* R correspond a Z/fZ et pi:Z -> R */ sum:=0; for i in [0..n] do for d_0 in S_2 do for r in S_1 do if delta() ne 0 then if (n+m-i) le precision then if n+m le precision then S:={h: h in [1..Integers()!(c/d_0)] | pi(h) eq pi(j*r) }; for h in S do sum:=sum+Binomial(n,i)*(a/c)^(n-i)*( (-1)^i )*( 1/((n+m-i+1)*(i+1)) )*(c/d_0)^(n+m-i)*delta()*d_0^(-i) *pBernoulli(n+m-i+1,h/(c/d_0))*( p^(n+m-i)*pBernoulli(i+1,p*a*h/(c/d_0))-p^(n+m)*pBernoulli(i+1,a*h/(c/d_0)) ); end for; else S:={h: h in [1..Integers()!(c/d_0)] | pi(h) eq pi(j*r) }; for h in S do sum:=sum+Binomial(n,i)*(a/c)^(n-i)*( (-1)^i )*( 1/( (n+m-i+1)*(i+1)) )*(c/d_0)^(n+m-i)*delta()*d_0^(-i) *pBernoulli(n+m-i+1,h/(c/d_0))*( p^(n+m-i)*pBernoulli(i+1,p*a*h/(c/d_0)) ); end for; end if; else sum:=sum; end if; end if; end for; end for; end for; return -12*(1/f^(n+m))*sum; end function; /* Formules pour les moments sur l'espace Zp^* x pZp: \int_X x^ny^m d\wt{\mu}_j(i\infty rightarrow a/c)(x,y) calcules dans Qp avec une certaine precision. */ moments_2_modify:=function(p,f,N_0,j,a,c,n,m,delta,precision,P) S_1:={h: h in [1..f]| GCD(h,f) eq 1}; S_2:=SequenceToSet(Divisors(N_0)); R,pi:=quo; /* R correspond a Z/fZ et pi:Z -> R */ sum:=0; for i in [0..n] do for d_0 in S_2 do for r in S_1 do if delta() ne 0 then if (n+m-i) le precision then if n+m le precision then S:={h: h in [1..Integers()!(c/d_0)] | pi(h) eq pi(j*r) }; for h in S do sum:=sum+Binomial(n,i)*(a/c)^(n-i)*( (-1)^i )*( 1/((n+m-i+1)*(i+1)) )*(c/d_0)^(n+m-i)*delta()*d_0^(-i) *Evaluate(P[n+m-i+1],Ratp( h/(c/d_0) ) )*( p^(n+m-i)*Evaluate(P[i+1], Ratp( p*a*h/(c/d_0) ) )-p^(n+m)*Evaluate(P[i+1], Ratp(a*h/(c/d_0)) ) ); end for; else S:={h: h in [1..Integers()!(c/d_0)] | pi(h) eq pi(j*r) }; for h in S do sum:=sum+Binomial(n,i)*(a/c)^(n-i)*( (-1)^i )*( 1/( (n+m-i+1)*(i+1)) )*(c/d_0)^(n+m-i)*delta()*d_0^(-i) *Evaluate(P[n+m-i+1], Ratp( h/(c/d_0) ) )*( p^(n+m-i)*Evaluate(P[i+1], Ratp( p*a*h/(c/d_0)) ) ); end for; end if; else sum:=sum; end if; end if; end for; end for; end for; return -12*(1/f^(n+m))*sum; end function; /* L'integrale: \int_{Zp^* x pZp} log_p x dmu{i\infty\rightarrow a/c}(x,y) */ integral_2:=function(p,f,N_0,j,a,c,delta,precision,P) poly:=approx_log_padic(p,precision); coeff:=Coefficients(poly); d:=#coeff; printf "Le polynome est de degree %o\n",d; sum:=0; for i in [1..d] do sum:=sum+coeff[i]*moments_2_modify(p,f,N_0,j,a,c,i-1,0,delta,precision,P); print(i); end for; return sum; end function; /* Calcul l'integrale: \int_{Zp x Zp^*} log_p y dmu_j{i\infty\rightarrow a_i/c_i}(x,y) pour une ensemble generateurs M_i=\M{a_i}{b_i}{c_i}{d_i} du groupe de congruence Gamma_0(fN_0) pour tous les j in (Z/f)^{\times}. La fonction retourne un vecteur de dimensions f. */ integral_1_all:=function(p,f,N_0,delta,precision,liste) Qp:=pAdicField(p,precision); printf "generate the list of peridic Bernoulli polynomials\n"; d:=(p-1)*( precision+Floor(Log(precision))+1 ); P:=[]; for i in [1..d] do P[i]:=BernoulliPolynomial2(i,liste); end for; S_1:={h:h in [1..f]| GCD(h,f) eq 1 and h le f/2}; G:=Gamma0(f*N_0); gen1:=Generators(G); gen2:=[]; w:=[]; s:=#gen1; for i in [1..s] do gen2[i]:=Matrix(2,2,Eltseq(gen1[i])); end for; for h in S_1 do v:=[]; for i in [1..s] do a:=gen2[i][1,1]; c:=gen2[i][2,1]; epsilon:=Sign(c); a:=epsilon*a; c:=epsilon*c; if c eq 0 then v[i]:=Qp!0; else v[i]:=integral_1(p,f,N_0,h,a,c,delta,precision,P); end if; end for; w[h]:=v; end for; W:=[]; for h in S_1 do W[h]:=w[h]; W[f-h]:=w[h]; end for; M1:=change1(W,f,Rationals()); PrintFile("moments_1.txt","M1:="); PrintFile("moments_1.txt",M1); PrintFile("moments_1.txt",";"); return W,gen2; end function; /* Calcul l'integrale: \int_{Zp^* x pZp} log_p x dmu_j{i\infty\rightarrow a_i/c_i}(x,y) pour une ensemble generateurs M_i=\M{a_i}{b_i}{c_i}{d_i} du groupe de congruence Gamma_0(fN_0) pour tous les j in (Z/f)^{\times}. La fonction retourne un vecteur de dimensions f. */ integral_2_all:=function(p,f,N_0,delta,precision,liste) Qp:=pAdicField(p,precision); printf "generate the list of peridic Bernoulli polynomials\n"; d:=(p-1)*( precision+Floor(Log(precision))+1 ); P:=[]; for i in [1..d] do P[i]:=BernoulliPolynomial2(i,liste); end for; S_1:={h:h in [1..f]| GCD(h,f) eq 1 and h le f/2}; G:=Gamma0(f*N_0); gen1:=Generators(G); gen2:=[]; w:=[]; s:=#gen1; for i in [1..s] do gen2[i]:=Matrix(2,2,Eltseq(gen1[i])); end for; for h in S_1 do v:=[]; for i in [1..s] do a:=gen2[i][1,1]; c:=gen2[i][2,1]; epsilon:=Sign(c); a:=epsilon*a; c:=epsilon*c; if c eq 0 then v[i]:=Qp!0; else v[i]:=integral_2(p,f,N_0,h,a,c,delta,precision,P); end if; end for; w[h]:=v; printf "L'indice %o est calcule\n",h; end for; W:=[]; for h in S_1 do W[h]:=w[h]; W[f-h]:=w[h]; end for; M2:=change2(W,f,Rationals()); PrintFile("moments_2.txt","M2:="); PrintFile("moments_2.txt",M2); PrintFile("moments_2.txt",";"); return W,gen2; end function; /* Determine si l'ideal de norme 2 de Q(\sqrt{D}) pour D = 1 mod 8 est principal ou non. */ is_2_principal:=function(D) K:=QuadraticField(D); O:=MaximalOrder(K); I:=ideal; true_false:=Sprint(IsPrincipal(I)); if true_false eq "true" then string,alpha:=IsPrincipal(I); return string,alpha; else return "false"; end if; end function; /* ----------------------------------------------------------------------------------------------------------- */ /* Soit une matrice \M{a}{b}{c}{d} et un vecteur [\tau,1] alors cette matrice retourne: (a\tau+b)/(c\tau+d) */ multiply:=function(M,v) return (M[1,1]*v[1]+M[1,2]*v[2])/(M[2,1]*v[1]+M[2,2]*v[2]); end function; /* Soit une matrice \M{a}{b}{c}{d} et un vecteur [\tau,1] alors cette matrice retourne: (a\tau+b)/(c\tau+d) */ multiply2:=function(M,v) return [(M[1,1]*v[1]+M[1,2]*v[2]),(M[2,1]*v[1]+M[2,2]*v[2])]; end function; /* Ce programme prend les entrees de tables Zi et les coerce lorsque c'est possible dans Qp */ swap_Qp_foot:=function(f,Qp,Z1,Z2,Z3,Z4) Z1:=change1(Z1,f,Rationals()); Z2:=change2(Z2,f,Rationals()); Z3:=change3(Z3,f,Rationals()); Z4:=change4(Z4,f,Rationals()); Z1:=change1(Z1,f,Qp); Z2:=change2(Z2,f,Qp); Z3:=change3(Z3,f,Qp); Z4:=change4(Z4,f,Qp); return Z1,Z2,Z3,Z4; end function; /* Cette fonction retourne la plus petite unite epsilon de K=Q(\sqrt{D}) pour laquelle epsilon = 1 mod 3. D correspond au discriminant du corps quadratique reel K. Le programme retourne une unite fondamental epsilon de K et un entier n. L'unite desiree est epsilon^n. */ smallest_unit:=function(D,f,N_0) K:=QuadraticField(D); R:=MaximalOrder(K); epsilon:=R!FundamentalUnit(K); I:=ideal; G:=UnitGroup(quo); n:=#G; Real:=RealField(100); DD:=Real!D; S:=SequenceToSet(Divisors(n)); for h in S do residue:=epsilon^h mod I; if (residue mod I) eq (R!1 mod I) and Norm(epsilon^h) eq 1 then a:=Eltseq(epsilon)[1]; /* l'unite est ecrite sous la forme a + b(sqrt{D}+1)/2 */ b:=Eltseq(epsilon)[2]; u:=Integers()!(a+b/2); v:=Integers()!(b/2); if u+v*SquareRoot(DD) gt 0 then return (u+v*t),h; else return -(u+v*t),h; end if; end if; end for; end function; /* Ce programme prend en entree un ideal entier I d'un corps quadratique K de discriminant D congruent a 1 modulo 8. Considerons l'ideal J=( 2,(1+sqrt{D})/2 )^2 alors la paire (I,I*J) a la propriete que I/(J*I) simeq Z/4. Il existe donc une Z-base de I et J*I telle que la matrice representant I est donnee par \M{a/4}{b/4}{0}{d} et la matrice representant I*J par \M{a}{b}{0}{d}. L'ideal I doit avoir ete construit a partir du corps K donne par K:=NUmberField(x^2-x+(1-D)/4) ou w:=( 1+sqrt{D} )/2 */ special_basis:=function(f,K,O,w,I,range) D:=Discriminant(K); L:=QuadraticField(D); n:=Floor(abs(f*(-SquareRoot(D)/2)))+1; J:=ideal; JJ:=ideal; I1:=I; I2:=I*J^2; I3:=I*JJ^2; B1:=Basis(I1,K); B2:=Basis(I2,K); B3:=Basis(I3,K); a1:=Eltseq(B1[1])[2]; b1:=Eltseq(B1[1])[1]; c1:=Eltseq(B1[2])[2]; d1:=Eltseq(B1[2])[1]; a2:=Eltseq(B2[1])[2]; b2:=Eltseq(B2[1])[1]; c2:=Eltseq(B2[2])[2]; d2:=Eltseq(B2[2])[1]; a3:=Eltseq(B3[1])[2]; b3:=Eltseq(B3[1])[1]; c3:=Eltseq(B3[2])[2]; d3:=Eltseq(B3[2])[1]; M1:=Matrix(Integers(),2,2,[a1,b1,c1,d1]); M2:=Matrix(Integers(),2,2,[a2,b2,c2,d2]); M3:=Matrix(Integers(),2,2,[a3,b3,c3,d3]); A1,T1:=HermiteForm(M1); A2,T2:=HermiteForm(M2); A3,T2:=HermiteForm(M3); print(A1); print(A2); print(A3); if A2[1,2] eq 0 then I2:=I3; end if; a:=Integers()!( A2[2,2]/A1[2,2] ); if a eq 2 then for k2 in [1..range] do for k1 in [k2*n..k2*n+range] do if ( (1+f*k1)+( f*k2*(-SquareRoot(D) + 1 )/2 ) ) gt 0 then while a eq 2 do B1:=Basis(I1*(1+f*k1+f*k2*w),K); B2:=Basis(I2*(1+f*k1+f*k2*w),K); a1:=Eltseq(B1[1])[2]; b1:=Eltseq(B1[1])[1]; c1:=Eltseq(B1[2])[2]; d1:=Eltseq(B1[2])[1]; a2:=Eltseq(B2[1])[2]; b2:=Eltseq(B2[1])[1]; c2:=Eltseq(B2[2])[2]; d2:=Eltseq(B2[2])[1]; M1:=Matrix(Integers(),2,2,[a1,b1,c1,d1]); M2:=Matrix(Integers(),2,2,[a2,b2,c2,d2]); A1,T1:=HermiteForm(M1); A2,T2:=HermiteForm(M2); print(Determinant(T1)); print(Determinant(T2)); a:=Integers()!( A2[2,2]/A1[2,2] ); end while; end if; end for; end for; end if; if a eq 1 then /* if fort */ for k in [0..range] do if A2[1,2]+k*A2[2,2] eq 4*A1[1,2] then B1:=A1; B2:=Matrix(2,2,[A2[1,1],4*A1[1,2],A2[2,1],A2[2,2] ]); return B1,B2; end if; if A2[1,2]-k*A2[2,2] eq 4*A1[1,2] then B1:=A1; B2:=Matrix(2,2,[A2[1,1],4*A1[1,2],A2[2,1],A2[2,2] ]); return B1,B2; end if; end for; end if; /* if fort */ if A1[1,2] eq 0 then for k in [0..range] do if k*A1[2,2] eq A2[1,2] then B1:=Matrix(2,2,[ A1[1,1],A2[1,2],A1[2,1],A1[2,2] ]); B2:=A2; return B1,B2; end if; end for; end if; if a eq 1 then /* if bravo */ for k in [1..range] do /* for all */ if A2[1,2] eq k*A2[2,2]+A2[1,2] then B1:=A1; B2:=Matrix(2,2,[A2[1,1],A2[1,2],A2[2,1],A2[2,2] ]); return B1,B2; end if; if A2[1,2] eq -k*A2[2,2]+A2[1,2] then B1:=A1; B2:=Matrix(2,2,[A2[1,1],A2[1,2],A2[2,1],A2[2,2] ]); return B1,B2; end if; end for; /* for allo */ end if; /* if bravo */ end function; /* Pour un element tau donne dans Q(\sqrt{D}), cette matrice calcule l'action de epsilon in O(f\infty)^{\tau} sur le reseau Z+tauZ par rapport a la base {1,tau} << Le reseau Z+tauZ doit etre un ideal >> L'element tau doit etre donne sous la forme a+b*sqrt{D} pour a et b dans Q. */ stabilisateur:=function(f,N_0,D,tau) K:=QuadraticField(D); a:=Eltseq(tau)[2]; /* On a ici que tau=a*sqrt{D}+b */ b:=Eltseq(tau)[1]; d:=Maximum(Denominator(a),Denominator(b)); a:=a*d; b:=b*d; epsilon,n:=smallest_unit(D,f,N_0); unit:=epsilon^n; u:=Eltseq(unit)[1]; /* L'unite epsilon peut etre ecrite comme u+v*sqrt{D} */ v:=Eltseq(unit)[2]; M:=Matrix(2,2,[u,D*v,v,u]); B:=Matrix(Rationals(),2,2,[a,b,0,d]); A:=B*M*B^(-1); n:=1; B:=A; while ( Integers()!B[2,1] mod (f*N_0) ) ne ( 0 mod (f*N_0) ) do B:=B*A; n:=n+1; end while; B:=Matrix(Integers(),2,2,Eltseq(B)); return B; end function; /* ------------------------------------------------------------------------------------------------------------------------ */ /* Cette fonction calcule les produits consecutifs: Soit w:=[g1,g2,...,gr] alors cette fonction retourne un vecteur de longueur r appele gamma ou gamma_i=g1*g2...*gi. */ gamma:=function(G,gen,w) n:=#w; A:=G!Matrix(2,2,[1,0,0,1]); gamma:=[]; for i in [1..n] do B:=A; for k in [1..i] do B:=B*gen[abs(w[k])]^Sign(w[k]); end for; gamma[i]:=B; end for; return gamma; end function; /* Ceci est une somme: \sum_{i=0}^{p-1} log_p(tau-i)\mu_{j}{i\infty\rightarrow A(i\infty)}(i+pZp) ou \mu est la mesure sur P^1(\QQ_p) */ sum1_integral_3:=function(p,f,N_0,j,delta,A,K,t,tau) M:=inv_matrices_distribution(p); sum:=0; for i in [0..p-1] do sum:=sum+log_padic2(K,t,tau-i)*period_p(p,f,N_0,j,delta,M[i+1]*A); end for; return sum; end function; /* Cette fonction calcule l'integrale suivant: \int_{(Z_p x Z_p)'} (x-\tau*y)d\mu_j{i\infty\rightarrow g(i\infty)}(x,y) Ici g est un des generateurs de Gamma_0(f*N_0). ou gen[r]:=g, K est egal Qp(\sqrt{D}) avec une certaine precision */ integrale1_tau:=function(Kp,tp,tau,gen,r,j,f,N_0,delta,V1,V2,V3,V4) p:=Prime(Kp); sum:=0; sum1:=0; sum2:=0; sum3:=0; sum:=sum+V1[j][r]; sum:=sum+V2[j][r]; A:=Matrix(2,2,Eltseq(gen[r])); sum1:=sum1+sum1_integral_3(p,f,N_0,j,delta,A,Kp,tp,tau); for i in [0..p-1] do for k in [1..( dimension(V3)[4]-1 )] do sum2:=sum2+1/(k*(tau-i)^k)*V3[r][j][i+1][k+1]; end for; end for; for k in [1..dimension(V4)[3]-1] do sum3:=sum3+tau^k/k*V4[r][j][k+1]; end for; return sum+sum1-sum2-sum3; end function; /* Cette fonction calcule l'integrale suivante: \int_{(Zp x Zp)'}log_p(x-\tau y)d\mu_j{i\infty\rightarrow g\infty}(x,y) ou w:=FindWord(G,g); */ integrale2_tau:=function(Kp,tp,tau,G,gen,w,j,f,N_0,delta,V1,V2,V3,V4) M:=gamma(G,gen,w); r:=#w; sum:=0; for i in [1..r] do if i ne r then A:=Matrix(2,2, Eltseq(M[r-i]) ); B:=A^(-1); else A:=Matrix(2,2,[1,0,0,1]); B:=A^(-1); end if; if abs(w[r-i+1]) ne 1 then theta:=multiply(B,[tau,1]); if Sign(w[r-i+1]) eq 1 then T:=integrale1_tau(Kp,tp,theta,gen,abs(w[r-i+1]), ( Integers()!(j*A[1,1]) ) mod f,f,N_0,delta,V1,V2,V3,V4); sum:=sum+T; else C:=Matrix(2,2, Eltseq(gen[abs(w[r-i+1])]) ); D:=C*B; theta:=multiply(D,[tau,1]); T:=integrale1_tau(Kp,tp,theta,gen,abs(w[r-i+1]), ( Integers()!(j*D[2,2]) ) mod f ,f,N_0,delta,V1,V2,V3,V4); sum:=sum-T; end if; else sum:=sum; end if; end for; return sum; end function; /* Cette fonction calcul \int_{i\infty}^{\gamma i\infty} \wt{F}_{2,\delta}(-r,z)dz pour \gamma un matrice dand G. Le vecteur w est donne ici par FindWord(G,gamma)=w. */ valuation:=function(p,f,N_0,r,delta,G,gen,w) M:=[]; valeurs:=table([f]); gam:=gamma(G,gen,w); for i in [1..#gen] do M[i]:=Matrix(Integers(),2,2,Eltseq(gen[i])); end for; for j in {h:h in [1..f]|GCD(h,f) eq 1} do for i in [1..#gen] do if i eq 1 then valeurs[j][1]:=0; else valeurs[j][i]:=period(p,f,N_0,j,delta,M[i]); end if; end for; end for; sum:=0; for i in [1..#w] do if i eq 1 then sum:=sum+Sign(w[i])*valeurs[r][abs(w[i])]; else sum:=sum+Sign(w[i])*valeurs[ ( Integers()!r*Eltseq(gam[i-1])[1] ) mod f][abs(w[i])]; end if; end for; return sum; end function; /* Calcule un certain entier dans Z/(p^2-1)Z pour un genrateur g de G et un element tau\in Kp ce programme calcule \sum_{(u,v) mod p\neq (0,0)} log_beta(u-v\tau)\mu_j{i\infty\rightarrow g(i\infty)}(u,v) */ int_log_beta:=function(Kp,tp,tau,g,p,f,N_0,j,delta) M:=Matrix(Integers(),2,2,Eltseq(g) ); a:=M[1,1]; c:=M[2,1]; a:=a*Sign(c); c:=c*Sign(c); sum:=0; for u in [0..p-1] do for v in [0..p-1] do if u ne 0 or v ne 0 then beta,k:=log_beta(Kp,tp,u-v*tau); sum:=sum+k*measure_ball(p,u,v,1,delta,f,N_0,a,c,j); end if; end for; end for; sum:=Integers()!sum; return sum,beta; end function; /* Calcule un certain entier dans Z/(p^2-1)Z pour un element quelconque g de G et un element tau\in Kp ce programme calcule \sum_{(u,v) mod p\neq (0,0)} log_beta(u-v\tau)\mu_j{i\infty\rightarrow g(i\i ici tau est un element de K_p de la forme tau=a+b*tp */ int_log_beta2:=function(Kp,tp,tau,G,g,p,f,N_0,j,delta) gen:=Generators(G); w:=FindWord(G,G!g); M:=gamma(G,gen,w); r:=#w; sum:=0; for i in [1..r] do if i ne r then A:=Matrix(2,2, Eltseq(M[r-i]) ); B:=A^(-1); else A:=Matrix(2,2,[1,0,0,1]); B:=A^(-1); end if; if abs(w[r-i+1]) ne 1 then theta:=multiply(B,[tau,1]); if Sign(w[r-i+1]) eq 1 then sum:=sum+int_log_beta(Kp,tp,theta,gen[abs(w[r-i+1])],p,f,N_0, (j*A[1,1]) mod f ,delta); else C:=Matrix(2,2, Eltseq(gen[abs(w[r-i+1])]) ); D:=C*B; theta:=multiply(D,[tau,1]); sum:=sum-int_log_beta(Kp,tp,theta,gen[abs(w[r-i+1])],p,f,N_0, (j*D[2,2]) mod f ,delta); end if; else sum:=sum; end if; end for; sum:=Integers()!sum; return sum,log_beta(Kp,tp,tau); end function; /* Soit g un element du groupe Gamma0(fN_0) et (u,v) une paire d'entiers donnes telle que 0=< u,v =< p-1 et (u,v) \neq (0,0) alors ce programme calcule \mu_j{i\infty\rightarrow g(i\infty)}( (u+pZp) x (v+pZp) ) */ measure_ball_fast:=function(g,u,v,p,f,N_0,j,delta) G:=Gamma0(f*N_0); gen:=Generators(G); gg:=G!g; w:=FindWord(G,G!gg); M:=gamma(G,gen,w); r:=#w; uu:=u; vv:=v; sum:=0; for i in [1..r] do if i ne r then A:=Matrix(2,2, Eltseq(M[r-i]) ); B:=A^(-1); else A:=Matrix(2,2,[1,0,0,1]); B:=A^(-1); end if; if abs(w[r-i+1]) ne 1 then if Sign(w[r-i+1]) eq 1 then uu:=( B[1,1]*u+B[1,2]*v ) mod p; vv:=( B[2,1]*u+B[2,2]*v ) mod p; C:=Matrix(Integers(),2,2,Eltseq( gen[abs(w[r-i+1])] ) ); sum:=sum+measure_ball(p, uu, vv, 1, delta, f, N_0, C[1,1], C[2,1], (j*B[2,2]) mod f ); else C:=Matrix(2,2, Eltseq(gen[abs(w[r-i+1])]) ); D:=C*B; uu:=( D[1,1]*u+D[1,2]*v ) mod p; vv:=( D[2,1]*u+D[2,2]*v ) mod p; sum:=sum - measure_ball(p, uu, vv, 1, delta, f, N_0, C[1,1], C[2,1], (j*D[2,2]) mod f ); end if; else sum:=sum; end if; end for; return sum; end function; /* Soit g un element de Gamma_0(fN_0) alors cette fonction calcule le GCD des entiers suivants: \mu_j{i\infty\rightarrow g(i\infty)}((u+pZp)x(v+pZp)) pour u,v parcourant les entiers entre 0 et p-1 tels que u,v ne sont pas tous deux congrus a 0 */ gcd_measure_ball:=function(g,p,delta,f,N_0,j) d:=Integers()!measure_ball_fast(g,1,0,p,f,N_0,j,delta); for u in [0..p-1] do for v in [0..p-1] do if (u ne 0) or (v ne 0) then d:=GCD(d,Integers()!measure_ball_fast(g,u,v,p,f,N_0,j,delta)); end if; end for; end for; return d; end function; /* Ce programme calcule le polynome minimal associe aux unites p-adiques construites pour un corps reel quadratique Q(\sqrt{D}) pour lequel p est inert et D=1 mod 8 et copremier a 3. Nos unites corresponde a l'unite modulaire associee au vecteur [2,-3,1]; */ calculs_polynome:=function(p,f,N_0,D,delta,j,V1,V2,V3,V4,range,precision) Qp:=pAdicField(p,precision); P:=PolynomialRing(Rationals()); K:=QuadraticField(D); KK:=NumberField(x^2-x+(1-D)/4); O:=MaximalOrder(KK); Kp:=UnramifiedExtension(Qp,x^2-D); R:=PolynomialRing(Kp); Real:=RealField(100); DD:=Real!D; G:=Gamma0(f*N_0); gen:=Generators(G); V1,V2,V3,V4:=swap_Qp_foot(f,Qp,V1,V2,V3,V4); I:=ideal; H,pi:=RayClassGroup(I,[1,2]); n:=#H; vecteur_ideaux:=[]; i:=1; for h in Set(H) do vecteur_ideaux[i]:=pi(h); i:=i+1; end for; vecteur_log_unites:=[]; vecteur_valuation:=[]; vecteur_racines_unite:=[]; vecteur_pgcd:=[]; vecteur_unites:=[]; for i in [1..n] do W1,W2:=special_basis(f,KK,O,w,vecteur_ideaux[i],range); a:=W2[2,2]/W1[2,2]; print(W1); print(W2); tau:=1/a*multiply(W1,[(t+1)/2,1]); print(tau); tau_p:=1/a*multiply(W1,[(tp+1)/2,1]); r:=Eltseq(tau)[1]; s:=Eltseq(tau)[2]; M:=stabilisateur(f,N_0,D,tau); print(multiply(M,[tau,1])); print( M[2,1]*tau+M[2,2] ); print(M); if Norm(M[2,1]*tau+M[2,2]) ne 1 then M:=M^(2); end if; if Norm(M[2,1]*tau+M[2,2]) ne 1 then return "Il y a erreur"; end if; if M[2,1]*(r*Real!SquareRoot(DD)+s)+M[2,2] gt 1 then M:=M; else M:=M^(-1); end if; print(M); word:=FindWord(G,G!M); print(word); vecteur_log_unites[i]:=integrale2_tau(Kp, tp, tau_p, G, gen, word, (Integers()!(j*W1[2,2])) mod f, f, N_0, delta, V1, V2, V3, V4); sum:=vecteur_log_unites[i]; printf "vecteur log unites est calcule\n"; vecteur_valuation[i]:=valuation(p, f, N_0, (Integers()!(j*W1[2,2])) mod f, delta, G, gen, word); val:=vecteur_valuation[i]; printf "vecteur valuation est calcule\n"; m,root:=int_log_beta2( Kp, tp, tau_p, G, M, p, f, N_0,(Integers()!(j*W1[2,2])) mod f , delta); vecteur_racines_unite[i]:=m; vecteur_pgcd[i]:=gcd_measure_ball(M, p, delta, f, N_0,(Integers()!(j*W1[2,2])) mod f); d:=vecteur_pgcd[i]; printf "vecteur pgcd est calcule\n"; print(sum); vecteur_unites[i]:=p^( Integers()!(val/d) )*root^(Integers()!(m/d))*Exp(sum/d); end for; poly:=1; for i in [1..n] do poly:=poly*(y-vecteur_unites[i]); end for; return poly,vecteur_valuation,vecteur_pgcd,vecteur_racines_unite,vecteur_log_unites; end function; /* Cette fonction retourne un polynome apres avoir fait une reconstruction rationel */ rational_reconstruction:=function(poly,p,D,precision) Ring:=Integers(p^precision); K:=QuadraticField(D); P:=PolynomialRing(K); coeff:=Coefficients(poly); d:=#coeff; yes_no_vecteur1:=[]; yes_no_vecteur2:=[]; poly2:=0; for i in [1..d] do w:=coeff[i]; v:=Append(Eltseq(w),[0]); a:=v[1]; b:=v[2]; a:=Rationals()!a; b:=Rationals()!b; n_a:=Valuation(a,p); n_b:=Valuation(b,p); if Sprint(n_a) ne "Infinity" or Sprint(n_b) ne "Infinity" then min:=Minimum(n_a,n_b); yes_no1,aa:=RationalReconstruction(Ring!Integers()!(a*p^(-min))); yes_no2,bb:=RationalReconstruction(Ring!Integers()!(b*p^(-min))); yes_no_vecteur1[i]:=yes_no1; yes_no_vecteur2[i]:=yes_no2; if Sprint(yes_no1) eq "true" and Sprint(yes_no2) eq "true" then poly2:=poly2+p^(min)*(aa+bb*t)*x^(i-1); end if; end if; end for; return poly2, yes_no_vecteur1,yes_no_vecteur2; end function; /* Ce programme calcule le polynome minimal associe aux unites p-adiques construites pour un corps reel quadratique Q(\sqrt{D}) pour lequel p est inert et D=1 mod 8 et copremier a 3. Nos unites corresponde a l'unite modulaire associee au vecteur [2,-3,1]; Pour un GCD=d */ calculs_polynome2:=function(p,f,N_0,D,delta,j,V1,V2,V3,V4,range,precision,d1,d2) Qp:=pAdicField(p,precision); P:=PolynomialRing(Rationals()); K:=QuadraticField(D); KK:=NumberField(x^2-x+(1-D)/4); O:=MaximalOrder(KK); Kp:=UnramifiedExtension(Qp,x^2-D); R:=PolynomialRing(Kp); Real:=RealField(100); G:=Gamma0(f*N_0); gen:=Generators(G); V1,V2,V3,V4:=swap_Qp_foot(f,Qp,V1,V2,V3,V4); I:=ideal; H,pi:=RayClassGroup(I,[1,2]); n:=#H; vecteur_ideaux:=[]; i:=1; for h in Set(H) do vecteur_ideaux[i]:=pi(h); i:=i+1; end for; vecteur_log_unites:=[]; vecteur_valuation:=[]; vecteur_racines_unite:=[]; vecteur_pgcd:=[]; vecteur_unites:=[]; for i in [1..n] do W1,W2:=special_basis(f,KK,O,w,vecteur_ideaux[i],range); a:=W2[2,2]/W1[2,2]; print(W1); print(W2); tau:=1/a*multiply(W1,[(t+1)/2,1]); print(tau); tau_p:=1/a*multiply(W1,[(tp+1)/2,1]); r:=Eltseq(tau)[1]; s:=Eltseq(tau)[2]; M:=stabilisateur(f,N_0,D,tau); print(multiply(M,[tau,1])); print( M[2,1]*tau+M[2,2] ); print(M); if Norm(M[2,1]*tau+M[2,2]) ne 1 then M:=M^(2); end if; if Norm(M[2,1]*tau+M[2,2]) ne 1 then return "Il y a erreur"; end if; if M[2,1]*(r*Real!SquareRoot(D)+s)+M[2,2] gt 1 then M:=M; else M:=M^(-1); end if; print(M); word:=FindWord(G,G!M); print(word); vecteur_log_unites[i]:=integrale2_tau(Kp, tp, tau_p, G, gen, word, (Integers()!(j*W1[2,2])) mod f, f, N_0, delta, V1, V2, V3, V4); sum:=vecteur_log_unites[i]; vecteur_valuation[i]:=valuation(p, f, N_0, (Integers()!(j*W1[2,2])) mod f, delta, G, gen, word); val:=vecteur_valuation[i]; m,root:=int_log_beta2( Kp, tp, tau_p, G, M, p, f, N_0,(Integers()!(j*W1[2,2])) mod f , delta); vecteur_racines_unite[i]:=m; vecteur_pgcd[i]:=gcd_measure_ball(M, p, delta, f, N_0,(Integers()!(j*W1[2,2])) mod f); print(sum); vecteur_unites[i]:=p^( Integers()!(val/d1) )*root^(Integers()!(m/d2))*Exp(sum/d1); end for; poly:=1; for i in [1..n] do poly:=poly*(y-vecteur_unites[i]); end for; return poly,vecteur_valuation,vecteur_pgcd,vecteur_racines_unite,vecteur_log_unites; end function; /* Ce programme prend en entree un ideal entier I d'un corps quadratique K de discriminant D congruent a 1 modulo 8 ainsi qu'un ideal I. Considerons l'ideal J=( 2,(1+sqrt{D})/2 )^2 alors la paire (I,I*J) a la propriete que I/(J*I) simeq Z/4. Il existe donc une Z-base de I et J*I telle que la matrice representant I est donnee par \M{a/4}{b/4}{0}{d} et la matrice representant I*J par \M{a}{b}{0}{d}. L'ideal I doit avoir ete construit a partir du corps K donne par K:=NUmberField(x^2-x+(1-D)/4) ou w:=( 1+sqrt{D} )/2 */ special_basis2:=function(f,K,O,w,I,range) printf "L'ideal est engendre par\n"; print(Basis(I,K)); D:=Discriminant(K); R,map:=quo; L:=QuadraticField(D); O_L:=MaximalOrder(L); A:=ideal; A:=A^2; n:=Floor(abs(3*(-SquareRoot(D)/2)))+1; J:=ideal; I1:=I; I2:=I*J^2; B1:=Basis(I1,K); B2:=Basis(I2,K); a1:=Eltseq(B1[1])[2]; b1:=Eltseq(B1[1])[1]; c1:=Eltseq(B1[2])[2]; d1:=Eltseq(B1[2])[1]; a2:=Eltseq(B2[1])[2]; b2:=Eltseq(B2[1])[1]; c2:=Eltseq(B2[2])[2]; d2:=Eltseq(B2[2])[1]; M1:=Matrix(Integers(),2,2,[a1,b1,c1,d1]); M2:=Matrix(Integers(),2,2,[a2,b2,c2,d2]); M1:=M1; M2:=M2; A1,T1:=HermiteForm(M1); A2,T2:=HermiteForm(M2); a:=Integers()!( A2[2,2]/A1[2,2] ); if a eq 1 and Sprint(is_principal(L,t,ideal,f)) eq "true" then sign:=1; e1:=w; e2:=1; else sign:=-1; end if; if a eq 2 then for k2 in [1..range] do for k1 in [k2*n..k2*n+range] do if ( (1+f*k1)+( f*k2*(-SquareRoot(D) + 1 )/2 ) ) gt 0 then while a eq 2 do B1:=Basis(I1*(1+f*k1+f*k2*w),K); B2:=Basis(I2*(1+f*k1+f*k2*w),K); a1:=Eltseq(B1[1])[2]; b1:=Eltseq(B1[1])[1]; c1:=Eltseq(B1[2])[2]; d1:=Eltseq(B1[2])[1]; a2:=Eltseq(B2[1])[2]; b2:=Eltseq(B2[1])[1]; c2:=Eltseq(B2[2])[2]; d2:=Eltseq(B2[2])[1]; M1:=Matrix(Integers(),2,2,[a1,b1,c1,d1]); M2:=Matrix(Integers(),2,2,[a2,b2,c2,d2]); A1,T1:=HermiteForm(M1); A2,T2:=HermiteForm(M2); a:=Integers()!( A2[2,2]/A1[2,2] ); end while; end if; end for; end for; if a eq 1 and Sprint(is_principal(L,t,ideal,f)) eq "true" then e1:=w; e2:=1; sign:=1; end if; if a eq 4 and Sprint( is_principal(L,t,ideal,f) ) eq "true" then e1:=w; e2:=1; B1:=A1; B2:=A2; sign:=1; end if; end if; /* Ici les ideaux B1 et B2 sont fixes et on varie la base */ /* On fait une premiere vague */ if a eq 2 or sign eq -1 then k:=0; while sign eq -1 and k le range do S:={h:h in [n+3*k,n+1+3*k,n+2+3*k]| (h mod 3) eq 1 }; v:=Setseq(S)[1]; P:=matrice_cusp(3,v); T:=Transpose(P); e1:=T[1,1]*w+T[1,2]; e2:=T[2,1]*w+T[2,2]; a1:=Eltseq(B1[1])[2]; b1:=Eltseq(B1[1])[1]; c1:=Eltseq(B1[2])[2]; d1:=Eltseq(B1[2])[1]; a2:=Eltseq(B2[1])[2]; b2:=Eltseq(B2[1])[1]; c2:=Eltseq(B2[2])[2]; d2:=Eltseq(B2[2])[1]; M1:=Matrix(Integers(),2,2,[a1,b1,c1,d1]); M2:=Matrix(Integers(),2,2,[a2,b2,c2,d2]); M1:=T*M1*T^(-1); M2:=T*M2*T^(-1); A1,T1:=HermiteForm(M1); A2,T2:=HermiteForm(M2); a:=Integers()!( A2[2,2]/A1[2,2] ); if a eq 4 and A2[1,2] eq 0 and Sprint( is_principal(K,w,ideal,f) ) eq "true" then sign:=1; B1:=A1; B2:=A2; end if; k:=k+1; end while; end if; /*------------------------------------------------------------------------------------------ */ if a eq 1 and sign eq 1 then /* if fort */ B1:=A1; B2:=Matrix(2,2,[A2[1,1],4*A1[1,2],A2[2,1],A2[2,2] ]); return B1,B2,e1,e2,1; end if; /* if fort */ /* ---------------------------------------------------------------------------------- */ if a eq 4 and sign eq 1 and B1[1,2] eq 0 then return B1,B2,e1,e2,1; end if; /* ----------------------------------------------------------------------------------- */ if a eq 4 and sign eq 1 then B1:=Matrix(2,2,[A2[1,1],A2[1,2],A2[2,1],A1[2,2]]); B2:=A2; return B1,B2,e1,e2,1; end if; /* ---------------------------------------------------------------------------------- */ R,map:=quo; epsilon:=FundamentalUnit(L); epsilon:=epsilon[1]+epsilon[2]*(2*w-1); if sign eq -1 and map(epsilon) ne map(1) and (Integers()!epsilon[2] mod 3) ne 0 then if Norm(epsilon) eq -1 then a:=Eltseq( map(epsilon^(-1)) )[1]; if Eltseq( map(epsilon^(-1)) )[2] eq 1 then b:=1; else b:=-1; end if; e1:=1; e2:=a+b*w; end if; if Norm(epsilon) eq 1 then a:=Eltseq( map(epsilon^(-1)) )[1]; if Eltseq( map(epsilon^(-1)) )[2] eq 1 then b:=1; else b:=-1; end if; e1:=1; e2:=a+3*( Floor(b/2*SquareRoot(D)) + 1 )+b*w; end if; a:=0; b:=1; c:=Integers()!e2[2]; d:=Integers()!e2[1]; T:=Matrix(2,2,[a,b,c,d]); for k2 in [1..range] do for k1 in [k2*n..k2*n+range] do tau_1:=multiply(M1,[w,1]); tau_2:=multiply(M2,[w,2]); /* ---------------------------------------------------- */ B1:=Basis(I1,K); B2:=Basis(I2,K); a1:=Eltseq(B1[1])[2]; b1:=Eltseq(B1[1])[1]; c1:=Eltseq(B1[2])[2]; d1:=Eltseq(B1[2])[1]; a2:=Eltseq(B2[1])[2]; b2:=Eltseq(B2[1])[1]; c2:=Eltseq(B2[2])[2]; d2:=Eltseq(B2[2])[1]; M1:=Matrix(Integers(),2,2,[a1,b1,c1,d1]); M2:=Matrix(Integers(),2,2,[a2,b2,c2,d2]); /* ------------------------------------------------- */ B1:=Basis(I1*(1+f*k1+f*k2*w),K); B2:=Basis(I2*(1+f*k1+f*k2*w),K); a1:=Eltseq(B1[1])[2]; b1:=Eltseq(B1[1])[1]; c1:=Eltseq(B1[2])[2]; d1:=Eltseq(B1[2])[1]; a2:=Eltseq(B2[1])[2]; b2:=Eltseq(B2[1])[1]; c2:=Eltseq(B2[2])[2]; d2:=Eltseq(B2[2])[1]; M1:=Matrix(Integers(),2,2,[a1,b1,c1,d1]); M2:=Matrix(Integers(),2,2,[a2,b2,c2,d2]); M1:=T*M1*T^(-1); /* Matrix M1 with respect to the basis {e1,e2} */ M2:=T*M2*T^(-1); /* Matrix M2 with respect to the basis {e1,e2} */ /* ------------------------------------------------------- */ A1,T1:=HermiteForm(M1); A2,T2:=HermiteForm(M2); a:=Integers()!( A2[2,2]/A1[2,2] ); rap2:=(e2*epsilon)[1]+(e2*epsilon)[2]*(SquareRoot(D)+1)/2; orientation:=Sign(rap2); /* ------------------------------------------------------ */ if a eq 1 and map( A1[2,2]*e2*epsilon*orientation ) eq map(1) then B1:=A1; B2:=Matrix(2,2,[A2[1,1],4*A1[1,2],A2[2,1],A2[2,2] ]); return B1,B2,(e1*epsilon),(e2*epsilon),orientation; end if; end for; end for; end if; /*------------------------------------On fait une derniere vague-----------------------*/ B1:=Basis(I1,K); B2:=Basis(I2,K); if a eq 2 or sign eq -1 then k:=1; while sign eq -1 and k le range do i:=1; while (i le 20) and sign eq -1 do if GCD(i,k) eq 1 then T:=matrice_cusp(k,i); e1:=T[1,1]*w+T[1,2]; e2:=T[2,1]*w+T[2,2]; a1:=Eltseq(B1[1])[2]; b1:=Eltseq(B1[1])[1]; c1:=Eltseq(B1[2])[2]; d1:=Eltseq(B1[2])[1]; a2:=Eltseq(B2[1])[2]; b2:=Eltseq(B2[1])[1]; c2:=Eltseq(B2[2])[2]; d2:=Eltseq(B2[2])[1]; M1:=Matrix(Integers(),2,2,[a1,b1,c1,d1]); M2:=Matrix(Integers(),2,2,[a2,b2,c2,d2]); M1:=T*M1*T^(-1); M2:=T*M2*T^(-1); A1,T1:=HermiteForm(M1); A2,T2:=HermiteForm(M2); a:=Integers()!( A2[2,2]/A1[2,2] ); if a eq 4 and Sprint( is_principal(K,w,ideal,f) ) eq "true" then sign:=1; B1:=A1; B2:=A2; end if; if a eq 1 and Sprint( is_principal(K,w,ideal,f) ) eq "true" then sign:=1; B1:=A1; B2:=A2; end if; end if; i:=i+1; end while; k:=k+1; end while; end if; /*------------------------------------------------------------------------------------------ */ if a eq 1 and sign eq 1 then /* if fort */ B1:=A1; B2:=Matrix(2,2,[A2[1,1],4*A1[1,2],A2[2,1],A2[2,2] ]); return B1,B2,e1,e2,1; end if; /* if fort */ /* -------------------------------------------------------------------------------------- */ if a eq 4 and sign eq 1 then B1:=Matrix(2,2,[A2[1,1],A2[1,2],A2[2,1],A1[2,2]]); B2:=A2; return B1,B2,e1,e2,1; end if; /* ---------------------------------------------------------------------------------- */ return "aucunne bonne base"; end function; /* Pour un tau dans K reduit par rapport a tau, ce programme renvoit une matrice dans Gamma0(f*N_0) inter Gamma1(f) qui stabilise le point tau L'element tau doit etre donne sous la form a+b*sqrt{D} pour a et b dans Q. */ stabilisateur2:=function(f,N_0,D,tau) K:=QuadraticField(D); a:=Eltseq(tau)[2]; /* On a ici que tau=a*sqrt{D}+b */ b:=Eltseq(tau)[1]; d:=Maximum(abs(Denominator(a)),abs(Denominator(b))); a:=a*d; b:=b*d; epsilon,n:=smallest_unit(D,f,N_0); unit:=epsilon^n; u:=Eltseq(unit)[1]; /* L'unite epsilon peut etre ecrite comme u+v*sqrt{D} */ v:=Eltseq(unit)[2]; M:=Matrix(2,2,[u,D*v,v,u]); B:=Matrix(Rationals(),2,2,[a,b,0,d]); A:=B*M*B^(-1); A:=Matrix(Integers(),2,2,Eltseq(A)); return A; end function; /* Ce programme calcule le polynome minimal associe aux unites p-adiques construites pour un corps reel quadratique Q(\sqrt{D}) pour lequel p est inert et D=1 mod 8 et copremier a 3. Nos unites corresponde a l'unite modulaire associee au vecteur [2,-3,1]; Si on ne veut pas considerer les racines de l'unites or on set "no_roots_unity=1" sinon on veut les considerer alors on set "no_roots_unity=0" */ calculs_polynome3:=function(p,f,N_0,D,delta,j,V1,V2,V3,V4,range,precision,no_roots_unity) Qp:=pAdicField(p,precision); P:=PolynomialRing(Rationals()); K:=QuadraticField(D); KK:=NumberField(x^2-x+(1-D)/4); O:=MaximalOrder(KK); Kp:=UnramifiedExtension(Qp,x^2-D); R:=PolynomialRing(Kp); Real:=RealField(200); wp:=(tp+1)/2; DD:=Real!D; G:=Gamma0(f*N_0); gen:=Generators(G); V1,V2,V3,V4:=swap_Qp_foot(f,Qp,V1,V2,V3,V4); I:=ideal; H,pi:=RayClassGroup(I,[1,2]); n:=#H; residue_ring,pi2:=quo; vecteur_ideaux:=[]; i:=1; for h in Set(H) do vecteur_ideaux[i]:=pi(h); i:=i+1; end for; vecteur_log_unites:=[]; vecteur_valuation:=[]; vecteur_racines_unite:=[]; vecteur_pgcd:=[]; vecteur_racines_unite_pgcd:=[]; vecteur_unites:=[]; for i in [1..n] do W1,W2,e1,e2,orientation:=special_basis2(f,KK,O,w,vecteur_ideaux[i],range); a:=W2[2,2]/W1[2,2]; printf "a est egal a %o\n\n",a; printf "W1 est egal a \n"; print(W1); printf "W2 est egal a \n"; print(W2); printf "e1 est egal a %o\n",e1; printf "e2 est egal a %o\n\n",e2; printf "Le signe associe a e2 est %o\n",orientation; /*------------------------------------------------------------- */ if a eq 1 then tau_w:=( W1[1,1]*e1+W1[1,2]*e2 )/(W1[2,2]*e2*orientation); tau:=tau_w[1]+tau_w[2]*(t+1)/2; epsilon:=Sign(tau[2]); tau:=epsilon*tau; tau_p:=tau_w[1]+tau_w[2]*(tp+1)/2; tau_p:=epsilon*tau_p; end if; if a eq 4 then tau_w:=W1[2,2]*e2/(W1[1,1]*e1+W1[1,2]*e2); tau:=tau_w[1]+tau_w[2]*(t+1)/2; epsilon:=Sign(tau[2]); tau:=epsilon*tau; tau_p:=tau_w[1]+tau_w[2]*(tp+1)/2; tau_p:=epsilon*tau_p; end if; /* --------------------------------------------------- */ r:=Eltseq(tau)[2]; s:=Eltseq(tau)[1]; printf "tau est egal a %o\n",tau; M:=stabilisateur2(f,N_0,D,tau); printf "Le stabilisateur de tau est \n"; print(M); printf "Verifions %o\n\n",multiply(M,[tau,1]); printf "L'unite est egale a %o\n",(M[2,1]*tau+M[2,2]); if Norm(M[2,1]*tau+M[2,2]) ne 1 then M:=M^(2); end if; if Norm(M[2,1]*tau+M[2,2]) ne 1 then return "Il y a erreur"; end if; if (M[2,1]*(r*Real!SquareRoot(DD)+s)+M[2,2]) gt 1 then M:=M; else M:=M^(-1); end if; printf "Apres orientation\n"; print(M); printf "Unite apres orientation\n"; print( M[2,1]*tau+M[2,2] ); word:=FindWord(G,G!M); print(word); vecteur_log_unites[i]:=integrale2_tau(Kp, tp, tau_p, G, gen, word, (Integers()!(j*W1[2,2])) mod f, f, N_0, delta, V1, V2, V3, V4); sum:=vecteur_log_unites[i]; printf "etape no 1\n"; vecteur_valuation[i]:=valuation(p, f, N_0, (Integers()!(j*W1[2,2])) mod f, delta, G, gen, word); val:=vecteur_valuation[i]; printf "etape no 2\n"; m,root:=int_log_beta2( Kp, tp, tau_p, G, M, p, f, N_0,(Integers()!(j*W1[2,2])) mod f , delta); printf "etape no 3\n"; vecteur_racines_unite[i]:=m; vecteur_pgcd[i]:=gcd_measure_ball(M, p, delta, f, N_0,(Integers()!(j*W1[2,2])) mod f); d:=vecteur_pgcd[i]; printf "etape no 4\n"; printf "d est egal a %o\n\n",d; vecteur_racines_unite_pgcd[i]:=GCD(m,d*48); print(sum); if no_roots_unity eq 0 then vecteur_unites[i]:=p^( Integers()!(val/d) )*root^(Integers()!(m/d))*Exp(sum/d); else vecteur_unites[i]:=p^( Integers()!(val/d) )*Exp(sum/d); end if; end for; poly:=1; for i in [1..n] do poly:=poly*(y-vecteur_unites[i]); end for; return poly,vecteur_valuation,vecteur_pgcd,vecteur_racines_unite,vecteur_racines_unite_pgcd,vecteur_log_unites; end function; /* Ce programme calcule le polynome minimal associe aux unites p-adiques construites pour un corps reel quadratique Q(\sqrt{D}) pour lequel p est inert et D=1 mod 8 et copremier a 3. Nos unites corresponde a l'unite modulaire associee au vecteur [2,-3,1]; Si on ne veut pas considerer les racines de l'unites or on set "no_roots_unity=1" sinon on veut les considerer alors on set "no_roots_unity=0" */ calculs_polynome4:=function(p,f,N_0,D,delta,j,V1,V2,V3,V4,range,precision,no_roots_unity,d) Qp:=pAdicField(p,precision); P:=PolynomialRing(Rationals()); K:=QuadraticField(D); KK:=NumberField(x^2-x+(1-D)/4); O:=MaximalOrder(KK); Kp:=UnramifiedExtension(Qp,x^2-D); R:=PolynomialRing(Kp); Real:=RealField(200); wp:=(tp+1)/2; DD:=Real!D; G:=Gamma0(f*N_0); gen:=Generators(G); V1,V2,V3,V4:=swap_Qp_foot(f,Qp,V1,V2,V3,V4); I:=ideal; H,pi:=RayClassGroup(I,[1,2]); n:=#H; residue_ring,pi2:=quo; vecteur_ideaux:=[]; i:=1; for h in Set(H) do vecteur_ideaux[i]:=pi(h); i:=i+1; end for; vecteur_log_unites:=[]; vecteur_valuation:=[]; vecteur_racines_unite:=[]; vecteur_pgcd:=[]; vecteur_racines_unite_pgcd:=[]; vecteur_unites:=[]; for i in [1..n] do W1,W2,e1,e2,orientation:=special_basis2(f,KK,O,w,vecteur_ideaux[i],range); a:=W2[2,2]/W1[2,2]; printf "a est egal a %o\n\n",a; printf "W1 est egal a \n"; print(W1); printf "W2 est egal a \n"; print(W2); printf "e1 est egal a %o\n",e1; printf "e2 est egal a %o\n\n",e2; printf "Le signe associe a e2 est %o\n",orientation; /*------------------------------------------------------------- */ if a eq 1 then tau_w:=( W1[1,1]*e1+W1[1,2]*e2 )/(W1[2,2]*e2*orientation); tau:=tau_w[1]+tau_w[2]*(t+1)/2; epsilon:=Sign(tau[2]); tau:=epsilon*tau; tau_p:=tau_w[1]+tau_w[2]*(tp+1)/2; tau_p:=epsilon*tau_p; end if; if a eq 4 then tau_w:=W1[2,2]*e2/(W1[1,1]*e1+W1[1,2]*e2); tau:=tau_w[1]+tau_w[2]*(t+1)/2; epsilon:=Sign(tau[2]); tau:=epsilon*tau; tau_p:=tau_w[1]+tau_w[2]*(tp+1)/2; tau_p:=epsilon*tau_p; end if; /* --------------------------------------------------- */ r:=Eltseq(tau)[2]; s:=Eltseq(tau)[1]; printf "tau est egal a %o\n",tau; M:=stabilisateur2(f,N_0,D,tau); printf "Le stabilisateur de tau est \n"; print(M); printf "Verifions %o\n\n",multiply(M,[tau,1]); printf "L'unite est egale a %o\n",(M[2,1]*tau+M[2,2]); if Norm(M[2,1]*tau+M[2,2]) ne 1 then M:=M^(2); end if; if Norm(M[2,1]*tau+M[2,2]) ne 1 then return "Il y a erreur"; end if; if (M[2,1]*(r*Real!SquareRoot(DD)+s)+M[2,2]) gt 1 then M:=M; else M:=M^(-1); end if; printf "Apres orientation\n"; print(M); printf "Unite apres orientation\n"; print( M[2,1]*tau+M[2,2] ); word:=FindWord(G,G!M); print(word); vecteur_log_unites[i]:=integrale2_tau(Kp, tp, tau_p, G, gen, word, (Integers()!(j*W1[2,2])) mod f, f, N_0, delta, V1, V2, V3, V4); sum:=vecteur_log_unites[i]; vecteur_valuation[i]:=valuation(p, f, N_0, (Integers()!(j*W1[2,2])) mod f, delta, G, gen, word); val:=vecteur_valuation[i]; m,root:=int_log_beta2( Kp, tp, tau_p, G, M, p, f, N_0,(Integers()!(j*W1[2,2])) mod f , delta); vecteur_racines_unite[i]:=m; vecteur_pgcd[i]:=gcd_measure_ball(M, p, delta, f, N_0,(Integers()!(j*W1[2,2])) mod f); printf "d est egal a %o\n\n",d; vecteur_racines_unite_pgcd[i]:=GCD(m,d*48); print(sum); if no_roots_unity eq 0 then vecteur_unites[i]:=p^( Integers()!(val/d) )*root^(Integers()!(m/d))*Exp(sum/d); else vecteur_unites[i]:=p^( Integers()!(val/d) )*Exp(sum/d); end if; end for; poly:=1; for i in [1..n] do poly:=poly*(y-vecteur_unites[i]); end for; return poly,vecteur_valuation,vecteur_pgcd,vecteur_racines_unite,vecteur_racines_unite_pgcd,vecteur_log_unites; end function; /* Ce programme permet de determiner l'entier n_D={1,2,3,4,6,12} */ rational_reconstruction2:=function(p,D,vecteur_valuation,vecteur_racines_unite,vecteur_log_unites,precision,no_roots_unity,d) Qp:=pAdicField(p,precision); P:=PolynomialRing(Rationals()); Kp:=UnramifiedExtension(Qp,x^2-D); R:=PolynomialRing(Kp); n:=#vecteur_valuation; root:=log_beta(Kp,tp,tp); vecteur_unites:=[]; for i in [1..n] do sum:=Kp!( Integers()!Coefficients(vecteur_log_unites[i])[1]+Integers()!Coefficients(vecteur_log_unites[i])[2]*tp ); val:=vecteur_valuation[i]; m:=vecteur_racines_unite[i]; if no_roots_unity eq 0 then vecteur_unites[i]:=p^( Integers()!(val/d) )*root^(Integers()!(m/d))*Exp(sum/d); else vecteur_unites[i]:=p^( Integers()!(val/d) )*Exp(sum/d); end if; end for; poly:=1; for i in [1..n] do poly:=poly*(y-vecteur_unites[i]); end for; return poly; end function; /* ---------------- Fonction test --------------------------- */ test:=function(f,K,O,w,I,range,a,b) D:=Discriminant(K); L:=QuadraticField(D); n:=Floor(abs(3*(-SquareRoot(D)/2)))+1; S:={h:h in [3*b+n,3*b+n+1,3*b+n+2]| (h mod 3) eq 1 }; v:=Setseq(S)[1]; P:=matrice_cusp(a,b); P:=Transpose(P); T:=P; e1:=T[1,1]*w+T[1,2]; e2:=T[2,1]*w+T[2,2]; J:=ideal; I1:=I; I2:=I*J^2; B1:=Basis(I1,K); B2:=Basis(I2,K); a1:=Eltseq(B1[1])[2]; b1:=Eltseq(B1[1])[1]; c1:=Eltseq(B1[2])[2]; d1:=Eltseq(B1[2])[1]; a2:=Eltseq(B2[1])[2]; b2:=Eltseq(B2[1])[1]; c2:=Eltseq(B2[2])[2]; d2:=Eltseq(B2[2])[1]; M1:=Matrix(Integers(),2,2,[a1,b1,c1,d1]); M2:=Matrix(Integers(),2,2,[a2,b2,c2,d2]); M1:=T*M1*T^(-1); M2:=T*M2*T^(-1); A1,T1:=HermiteForm(M1); A2,T2:=HermiteForm(M2); return A1,A2,e1,e2; end function; /* ----------------------------------------------------------------------------------- */ swap_coeff_poly:=function(poly,p,D,precision) Qp:=pAdicField(p,precision); P:=PolynomialRing(Qp); Kp:=UnramifiedExtension(Qp,x^2-D); PP:=PolynomialRing(Kp); coeff:=Coefficients(poly); d:=#coeff; poly2:=0; for i in [1..d] do w:=coeff[i]; v:=Append(Eltseq(w),[0]); a:=v[1]; b:=v[2]; a:=Qp!a; b:=Qp!b; poly2:=poly2+(a+b*tp)*z^(i-1); end for; return poly2; end function; /* Ce programme suppose que f=3 et cacul les divers moments qu'il met dans un fichier texte. */ calculs_moments:=function(p,r,v,precision) liste:=Bernoulli_liste(p,precision); f:=3; N_0:=4; S_2,S_1,delta:=Delta2(p,f,N_0,r,v); M1:=integral_1_all(p, f, N_0, delta, precision,liste); M2:=integral_2_all(p, f, N_0, delta, precision,liste); T2_0,T2_1,T3_0,T3_1,T5_0,T5_1:=negative_moments_gamma0_12(p, delta, precision, precision,liste); J2_1,J3_1,J4_1,J5_1:=negative_moments_gamma0_12_inverse(p, delta, precision, precision,liste); Qp:=pAdicField(p,precision); T2_0,T2_1,T3_0,T3_1,T5_0,T5_1,J2_1,J3_1,J4_1,J5_1:=swap_all_Qp(Qp, T2_0, T2_1, T3_0, T3_1, T5_0, T5_1, J2_1, J3_1, J4_1, J5_1); M3:=moments_shifted3(p, f, N_0, delta, precision, precision, T2_0, T2_1, T3_0, T3_1, T5_0, T5_1,liste); V3:=change3(M3,f,Qp); M4:=inverse_moments3(p, f, N_0, delta, precision, precision, J2_1, J3_1, J4_1, J5_1, T2_1, T5_0,liste); V4:=change4(M4,f,Qp); V1:=change1(M1,f,Qp); V2:=change2(M2,f,Qp); printf "poly,vecteur_valuation,vecteur_pgcd,vecteur_racines_unite,vecteur_log_unites:=calculs_polynome3(p,f,N_0,D,delta,j,V1,V2,V3,V4,range,precision,no_roots_unity)\n"; printf "poly:=rational_reconstruction2(p, D, vecteur_valuation, vecteur_racines_unite, vecteur_log_unites, precision2, no_roots_unity, d)"; return V1,V2,V3,V4; end function; /* Meme chose que le programme precedent sauf qu'il faut loader prealablement les fichiers de moments calcule a partir du programme precedent. Ceci permet de sauver de l'espace de memoire */ calculs_moments_2:=function(p,r,v,precision,M1,M2,M3,M4) f:=3; N_0:=4; S_2,S_1,delta:=Delta2(p,f,N_0,r,v); Qp:=pAdicField(p,precision); V3:=change3(M3,f,Qp); V4:=change4(M4,f,Qp); V1:=change1(M1,f,Qp); V2:=change2(M2,f,Qp); printf "poly,vecteur_valuation,vecteur_pgcd,vecteur_racines_unite,vecteur_log_unites:=calculs_polynome(p, f, N_0, D, delta, j, V1, V2, V3, V4, range, precision)\n"; printf "poly:=rational_reconstruction2(p, D, vecteur_valuation, vecteur_racines_unite, vecteur_log_unites, precision, no_roots_unity, d)"; return V1,V2,V3,V4; end function; /* Ce programme prend en entree un polynome avec des coefficients dans Qp et retourne en sortie un polynome avec des coefficients dans Z */ poly_convert:=function(poly) v:=Coefficients(poly); n:=#v; w:=[]; P:=PolynomialRing(Rationals()); poly_m:=P!0; for i in [1..n] do poly_m:=poly_m+x^(i-1)*(Rationals()!v[i]); end for; return poly_m; end function; /* Ce programme prend en entree un polynome avec des coefficients dans Qp(\sqrt{D) et retourne en sortie un polynome avec des coefficients dans Q(\sqrt{D}) */ poly_convert2:=function(poly,D) K:=QuadraticField(D); P:=PolynomialRing(K); v:=Coefficients(poly); n:=#v; w:=[]; poly_m:=P!0; for i in [1..n] do seq:=Eltseq(v[i]); poly_m:=poly_m+x^(i-1)*(Rationals()!seq[1]+t*Rationals()!seq[2]); end for; return poly_m; end function; /* Ce programme prend en entree un polynome avec des coefficients dans Qp(\sqrt{D) et retourne en sortie un polynome avec des coefficients dans Q(\sqrt{D}) sous une forme tenant compte de la factorization */ poly_convert3:=function(poly,D) poly2:=poly_convert2(2*poly,D); poly2:=1/2*poly2; coeffs:=Coefficients(poly2); n:=#coeffs; for i in [1..n] do den:=Denominator(coeffs[i]); num:=Numerator(coeffs[i]); printf "%o",Factorization(den); printf "(%o)x^{%o}\n",num,n-i; end for; return "fin"; end function; /* Ce programme calcul les moments suivants: \mu_r\{\infty\rightarrow g_j\infty\}(Z_p) and \mu_r\{\infty\rightarrow M_ig_j\infty\}(Z_p) */ moments_5_6:=function(p,f,N_0,delta) G:=Gamma0(f*N_0); gen:=Generators(G); M:=inv_matrices_distribution(p); n:=#gen; gen2:=[]; for i in [1..n] do gen2[i]:=Matrix(Integers(),2,2,Eltseq(gen[i])); end for; Tab1:=table([f,n]); for j in [1..f] do for k in [1..n] do if Gcd(j,f) eq 1 and gen2[k][2,1] ne 0 then Tab1[j][k][1]:=period_p(p,f,N_0,j,delta,gen2[k]); else Tab1[j][k][1]:=0; end if; end for; end for; Tab2:=table([f,n,p]); for j in [1..f] do for k in [1..n] do for l in [1..p] do if Gcd(j,f) eq 1 and gen2[k][2,1] ne 0 then Tab2[j][k][l][1]:=period_p(p,f,N_0,j,delta,M[l]*gen2[k]); else Tab2[j][k][l][1]:=0; end if; end for; end for; end for; PrintFile("moments_5.txt","M5:="); PrintFile("moments_5.txt",Tab1); PrintFile("moments_5.txt",";"); PrintFile("moments_6.txt","M6:="); PrintFile("moments_6.txt",Tab2); PrintFile("moments_6.txt",";"); return Tab1,Tab2; end function; /* Ce programme calcule \mu_j{i\infty\rightarrow g(i\infty)}(u,v) pour tous les generateurs g de G et pour toutes les paires (u,v) mod p */ moments_7:=function(p,f,N_0,delta) G:=Gamma0(f*N_0); gen:=Generators2(G); n:=#gen; liste:=table([n,f,p,p]); for i in [1..n] do for j in [1..f] do for u in [0..p-1] do for v in [0..p-1] do a:=gen[i][1,1]; c:=gen[i][2,1]; a:=Integers()!a*Sign(c); c:=Integers()!c*Sign(c); if (u ne 0 or v ne 0) and (Gcd(f,j) eq 1) and c ne 0 then liste[i][j][u+1][v+1][1]:=measure_ball(p,u,v,1,delta,f,N_0,a,c,j); else liste[i][j][u+1][v+1][1]:=0; end if; end for; end for; end for; end for; PrintFile("moments_7.txt","M7:="); PrintFile("moments_7.txt",liste); PrintFile("moments_7.txt",";"); return liste; end function; /* Nouveaux changements Nouveaux changements Nouveaux changements Nouveaux changements Nouveaux changements Nouveaux changements Nouveaux changements Nouveaux changements Nouveaux changements Nouveaux changements Nouveaux changements Nouveaux changements Nouveaux changements Nouveaux changements */ /* Calcule un certain entier dans Z/(p^2-1)Z pour un genrateur g de G et un element tau\in Kp ce programme calcule \sum_{(u,v) mod p\neq (0,0)} log_beta(u-v\tau)\mu_j{i\infty\rightarrow g(i\infty)}(u,v) r est l'indice correspondant au generateur g du group Gamma_0(fN_0). */ int_log_beta_speed:=function(Kp,tp,tau_p,r,p,f,N_0,j,delta,V7) R:=RingOfIntegers(Kp); F,h:=ResidueClassField(R); p:=Characteristic(F); w:=PrimitiveElement(F); a:=Integers()!Eltseq(tau_p)[1] mod p; b:=Integers()!Eltseq(tau_p)[2] mod p; a:=Integers()!a; b:=Integers()!b; pi:=homIntegers(p)|>; sum:=0; for k in [0..p^2-2] do u:=Integers()!Eltseq(w^k)[1]; /* les coefficients u,v sont par rapport a tau_p*/ v:=Integers()!Eltseq(w^k)[2]; if u ne 0 or v ne 0 then uu:=Integers()! ( Integers()!(u-pi(b^(-1))*a*v) mod p ); vv:=Integers()! ( ( Integers()!(-v*pi(b^(-1))) ) mod p ); sum:=sum+k*V7[r][j][uu+1][vv+1][1]; /* subtle but ok */ end if; end for; sum:=Integers()!sum; return sum; end function; /* Calcule un certain entier dans Z/(p^2-1)Z pour un element quelconque g de G et un element tau\in Kp ce programme calcule \sum_{(u,v) mod p\neq (0,0)} log_beta(u-v\tau)\mu_j{i\infty\rightarrow g(i\i */ int_log_beta2_speed:=function(Kp,tp,tau,G,g,p,f,N_0,j,delta,V7) gen:=Generators(G); w:=FindWord(G,G!g); M:=gamma(G,gen,w); r:=#w; sum:=0; for i in [1..r] do if i ne r then A:=Matrix(2,2, Eltseq(M[r-i]) ); B:=A^(-1); else A:=Matrix(2,2,[1,0,0,1]); B:=A^(-1); end if; if abs(w[r-i+1]) ne 1 then theta:=multiply(B,[tau,1]); if Sign(w[r-i+1]) eq 1 then sum:=sum+int_log_beta_speed(Kp,tp,theta,abs(w[r-i+1]),p,f,N_0, (j*A[1,1]) mod f ,delta,V7); else C:=Matrix(2,2, Eltseq(gen[abs(w[r-i+1])]) ); D:=C*B; theta:=multiply(D,[tau,1]); sum:=sum-int_log_beta_speed(Kp,tp,theta,abs(w[r-i+1]),p,f,N_0, (j*D[2,2]) mod f ,delta,V7); end if; else sum:=sum; end if; end for; /* epsilon:=Sign(sum); printf "La sum est\n"; print(sum); sum:=Integers()!(sum mod (p^2-1)); b:=Integers()!((p^2-1)/2); if sum le b then sum:=sum; else sum:=sum-(p^2-1); end if; if sum eq b then sum:=epsilon*sum; end if; */ return sum,log_beta(Kp,tp,tau); end function; /* logarithme p-adic pour des elements de Kp=Qp(tp) ou D est un non residu modulo p et tp est tel que tp^2=D. */ log_padic2_speed:=function(Kp,tp,a) n:=Precision(Kp); p:=Prime(Kp); P:=PolynomialRing(Kp); b:=Kp!a; f:=y^(p^2-1)-1; s_1:=Valuation(Rationals()!Coefficients(b)[1],p); s_2:=Valuation(Rationals()!Coefficients(b)[2],p); s:=Minimum(s_1,s_2); b:=b/p^s; r1:=Integers()!Coefficients(b)[1] mod p; r2:=Integers()!Coefficients(b)[2] mod p; if r1 eq 0 then r1:=p; end if; if r2 eq 0 then r2:=p; end if; b:=b/(HenselLift(f,Kp!(r1+r2*tp),n)); return Log(b); end function; /* Ceci est une somme: \sum_{i=0}^{p-1} log_p(tau-i)\mu_{j}{i\infty\rightarrow A(i\infty)}(i+pZp) ou \mu est la mesure sur P^1(\QQ_p) */ sum1_integral_3_speed:=function(p,f,N_0,j,delta,r,K,t,tau,V5,V6) M:=inv_matrices_distribution(p); sum:=0; for i in [0..p-1] do sum:=sum+log_padic2_speed(K,t,tau-i)*V6[j][r][i+1][1]; end for; return sum; end function; /* Cette fonction calcule l'integrale suivant: \int_{(Z_p x Z_p)'} (x-\tau*y)d\mu_j{i\infty\rightarrow g(i\infty)}(x,y) Ici g est un des generateurs de Gamma_0(f*N_0). ou gen[r]:=g, K est egal Qp(\sqrt{D}) avec une certaine precision */ integrale1_tau_speed:=function(Kp,tp,tau,gen,r,j,f,N_0,delta,V1,V2,V3,V4,V5,V6) p:=Prime(Kp); sum:=0; sum1:=0; sum2:=0; sum3:=0; sum:=sum+V1[j][r]; sum:=sum+V2[j][r]; sum1:=sum1+sum1_integral_3_speed(p,f,N_0,j,delta,r,Kp,tp,tau,V5,V6); for i in [0..p-1] do for k in [1..( dimension(V3)[4]-1 )] do sum2:=sum2+1/(k*(tau-i)^k)*V3[r][j][i+1][k+1]; end for; end for; for k in [1..dimension(V4)[3]-1] do sum3:=sum3+tau^k/k*V4[r][j][k+1]; end for; return sum+sum1-sum2-sum3; end function; /* Cette fonction calcule l'integrale suivante: \int_{(Zp x Zp)'}log_p(x-\tau y)d\mu_j{i\infty\rightarrow g\infty}(x,y) ou w:=FindWord(G,g); */ integrale2_tau_speed:=function(Kp,tp,tau,G,gen,w,j,f,N_0,delta,V1,V2,V3,V4,V5,V6) M:=gamma(G,gen,w); r:=#w; sum:=0; for i in [1..r] do if i ne r then A:=Matrix(2,2, Eltseq(M[r-i]) ); B:=A^(-1); else A:=Matrix(2,2,[1,0,0,1]); B:=A^(-1); end if; if abs(w[r-i+1]) ne 1 then theta:=multiply(B,[tau,1]); if Sign(w[r-i+1]) eq 1 then T:=integrale1_tau_speed(Kp,tp,theta,gen,abs(w[r-i+1]), ( Integers()!(j*A[1,1]) ) mod f,f,N_0,delta,V1,V2,V3,V4,V5,V6); sum:=sum+T; else C:=Matrix(2,2, Eltseq(gen[abs(w[r-i+1])]) ); D:=C*B; theta:=multiply(D,[tau,1]); T:=integrale1_tau_speed(Kp,tp,theta,gen,abs(w[r-i+1]), ( Integers()!(j*D[2,2]) ) mod f ,f,N_0,delta,V1,V2,V3,V4,V5,V6); sum:=sum-T; end if; else sum:=sum; end if; end for; return sum; end function; /* Ce programme calcule le polynome minimal associe aux unites p-adiques construites pour un corps reel quadratique Q(\sqrt{D}) pour lequel p est inert et D=1 mod 8 et copremier a 3. Nos unites corresponde a l'unite modulaire associee au vecteur [2,-3,1]; */ calculs_polynome_speed:=function(p,f,N_0,D,delta,j,V1,V2,V3,V4,V5,V6,V7,range,precision) Qp:=pAdicField(p,precision); P:=PolynomialRing(Rationals()); K:=QuadraticField(D); KK:=NumberField(x^2-x+(1-D)/4); O:=MaximalOrder(KK); Kp:=UnramifiedExtension(Qp,x^2-D); R:=PolynomialRing(Kp); Real:=RealField(100); DD:=Real!D; G:=Gamma0(f*N_0); gen:=Generators(G); V1,V2,V3,V4:=swap_Qp_foot(f,Qp,V1,V2,V3,V4); I:=ideal; H,pi:=RayClassGroup(I,[1,2]); n:=#H; vecteur_ideaux:=[]; i:=1; for h in Set(H) do vecteur_ideaux[i]:=pi(h); i:=i+1; end for; vecteur_log_unites:=[]; vecteur_valuation:=[]; vecteur_racines_unite:=[]; vecteur_pgcd:=[]; vecteur_unites:=[]; for i in [1..n] do W1,W2:=special_basis(f,KK,O,w,vecteur_ideaux[i],range); a:=W2[2,2]/W1[2,2]; print(W1); print(W2); tau:=1/a*multiply(W1,[(t+1)/2,1]); print(tau); tau_p:=1/a*multiply(W1,[(tp+1)/2,1]); r:=Eltseq(tau)[1]; s:=Eltseq(tau)[2]; M:=stabilisateur(f,N_0,D,tau); print(multiply(M,[tau,1])); print( M[2,1]*tau+M[2,2] ); print(M); if Norm(M[2,1]*tau+M[2,2]) ne 1 then M:=M^(2); end if; if Norm(M[2,1]*tau+M[2,2]) ne 1 then return "Il y a erreur"; end if; if M[2,1]*(r*Real!SquareRoot(DD)+s)+M[2,2] gt 1 then M:=M; else M:=M^(-1); end if; print(M); word:=FindWord(G,G!M); print(word); vecteur_log_unites[i]:=integrale2_tau_speed(Kp, tp, tau_p, G, gen, word, (Integers()!(j*W1[2,2])) mod f, f, N_0, delta, V1, V2, V3, V4,V5,V6); sum:=vecteur_log_unites[i]; printf "vecteur log unites est calcule\n"; vecteur_valuation[i]:=valuation(p, f, N_0, (Integers()!(j*W1[2,2])) mod f, delta, G, gen, word); val:=vecteur_valuation[i]; printf "vecteur valuation est calcule\n"; m,root:=int_log_beta2_speed( Kp, tp, tau_p, G, M, p, f, N_0,(Integers()!(j*W1[2,2])) mod f , delta,V7); vecteur_racines_unite[i]:=m; printf "vecteur racine de l'unite est calcule\n"; d:=12; vecteur_unites[i]:=p^( Integers()!(val/d) )*root^(Integers()!(m/d))*Exp(sum/d); end for; poly:=1; for i in [1..n] do poly:=poly*(y-vecteur_unites[i]); end for; return poly,vecteur_valuation,vecteur_racines_unite,vecteur_log_unites; end function; calculs_moments_3:=function(p,r,v,precision,M1,M2,M3,M4,M5,M6,M7) f:=3; N_0:=4; S_2,S_1,delta:=Delta2(p,f,N_0,r,v); Qp:=pAdicField(p,precision); V3:=change3(M3,f,Qp); V4:=change4(M4,f,Qp); V1:=change1(M1,f,Qp); V2:=change2(M2,f,Qp); V5:=M5; V6:=M6; V7:=M7; printf "poly,vecteur_valuation,vecteur_racines_unite,vecteur_log_unites:=calculs_polynome_speed_3(p, f, N_0, D, delta, j, V1, V2, V3, V4,V5,V6,V7, precision)\n"; printf "poly:=rational_reconstruction2(p, D, vecteur_valuation, vecteur_racines_unite, vecteur_log_unites, precision, no_roots_unity, d)"; return V1,V2,V3,V4,V5,V6,V7; end function; /* LE CALCULS MOMENTS SPEED */ calculs_moments_speed:=function(p,N_0,r,v,precision) liste:=Bernoulli_liste(p,precision); f:=3; S_2,S_1,delta:=Delta2(p,f,N_0,r,v); M1:=integral_1_all(p, f, N_0, delta, precision,liste); M2:=integral_2_all(p, f, N_0, delta, precision,liste); T2_0,T2_1,T3_0,T3_1,T5_0,T5_1:=negative_moments_gamma0_12(p, delta, precision, precision,liste); J2_1,J3_1,J4_1,J5_1:=negative_moments_gamma0_12_inverse(p, delta, precision, precision,liste); Qp:=pAdicField(p,precision); T2_0,T2_1,T3_0,T3_1,T5_0,T5_1,J2_1,J3_1,J4_1,J5_1:=swap_all_Qp(Qp, T2_0, T2_1, T3_0, T3_1, T5_0, T5_1, J2_1, J3_1, J4_1, J5_1); M3:=moments_shifted3(p, f, N_0, delta, precision, precision, T2_0, T2_1, T3_0, T3_1, T5_0, T5_1,liste); V3:=change3(M3,f,Qp); M4:=inverse_moments3(p, f, N_0, delta, precision, precision, J2_1, J3_1, J4_1, J5_1, T2_1, T5_0,liste); V4:=change4(M4,f,Qp); V1:=change1(M1,f,Qp); V2:=change2(M2,f,Qp); V5,V6:=moments_5_6(p,f,N_0,delta); V7:=moments_7(p,f,N_0,delta); printf "poly,vecteur_valuation,vecteur_racines_unite,vecteur_log_unites:=calculs_polynome_speed_3(p, f, N_0, D, delta, j, V1, V2, V3, V4,V5,V6,V7, precision)\n"; printf "poly:=rational_reconstruction2(p, D, vecteur_valuation, vecteur_racines_unite, vecteur_log_unites, precision, no_roots_unity, d)"; return V1,V2,V3,V4,V5,V6,V7; end function; /* Soit I un O_K-ideal de K=Q(w) ou w=(1+\sqrt{D})/2 ou D est un discriminant. Alors ce programme retourne un entier 1=< r<= f copremier a f et un element tau\in Q(w) tels que r\Lambda_{\tau} est equivalent a I. */ special_pair:=function(f,K,O,w,I) D:=Discriminant(K); B:=Basis(I,K); a:=Eltseq(B[1])[2]; b:=Eltseq(B[1])[1]; c:=Eltseq(B[2])[2]; d:=Eltseq(B[2])[1]; M1:=Matrix(2,2,[a,b,c,d]); if c eq 0 then r:=Integers()!(Integers()!M1[2,2] mod f); tau:=(M1[1,1]*w+M1[1,2])/M1[2,2]; return r,tau; end if; if a eq 0 then M1:=Matrix(2,2,[c,d,a,b]); r:=Integers()!(Integers()!M1[2,2] mod f); tau:=(M1[1,1]*w+M1[1,2])/M1[2,2]; return r,tau; end if; g,u,v:=ExtendedGreatestCommonDivisor(Integers()!-c,Integers()!-a); /* u*(-c)+v*(-a)=1 */ M2:=Matrix(2,2,[v,u,c,-a]); M3:=M2*M1; if Sign(M3[2,2]) eq -1 then M3:=-M3; end if; r:=Integers()!(Integers()!M3[2,2] mod f); tau:=(M3[1,1]*w+M3[1,2])/M3[2,2]; return r,tau; end function; /* Soit J:=ideal et I un ideal quelconque alors ce programme donne une base triangulee de (I,I*J^2) sous la forme: [w1,w2], [w1,4*w2] Remarquons que l'anneau des endomorphismes de w2/w1 et 4w2/w1 est l'ordre maximal O_K. */ special_basis_2:=function(f,K,O,w,I) D:=Discriminant(K); L:=QuadraticField(D); J:=ideal; I1:=I; I2:=I*J^2; B1:=Basis(I1,K); B2:=Basis(I2,K); a1:=Eltseq(B1[1])[2]; b1:=Eltseq(B1[1])[1]; c1:=Eltseq(B1[2])[2]; d1:=Eltseq(B1[2])[1]; a2:=Eltseq(B2[1])[2]; b2:=Eltseq(B2[1])[1]; c2:=Eltseq(B2[2])[2]; d2:=Eltseq(B2[2])[1]; M1:=Matrix(Integers(),2,2,[a1,b1,c1,d1]); M2:=Matrix(Integers(),2,2,[a2,b2,c2,d2]); A1,T1:=HermiteForm(M1); /* T1*M1=A1 */ A2,T2:=HermiteForm(M2); /* T2*M2=A2 */ A1:=Matrix(Rationals(),2,2,Eltseq(A1)); A2:=Matrix(Rationals(),2,2,Eltseq(A2)); A3:=A2*A1^(-1); tau1:=(A1[1,1]*w+A1[1,2]); tau2:=A1[2,2]*1; if A3[1,1] eq 1 then w1:=A3[1,1]*tau1+A3[1,2]*tau2; w2:=tau2; return [w1,w2],[w1,4*w2]; end if; if A3[1,1] eq 4 then k:=Integers()!(Integers()!A3[1,2] mod 4); w1:=tau2; w2:=tau1+(A3[1,2]-k)/4*tau2; return [w1,w2],[w1,4*w2]; end if; b:=Integers()!A3[1,2]; /* b est necessairement impair */ g,u,v:=ExtendedGreatestCommonDivisor(-2,Integers()!-b); /* (-2)*u+(-b)*v=1 */ M:=Matrix(2,2,[-b,u,2,v]); v:=multiply2(M^(-1),[tau1,tau2]); v1:=v[1]; v2:=v[2]; MM:=A3*M; if MM[1,2] eq -1 then MM:=Matrix(2,2,[MM[1,1],1,4,MM[2,2]]); end if; /* On swap */ MM:=Matrix(2,2,[MM[2,1],MM[2,2],MM[1,1],MM[1,2]]); if MM[1,1] eq 4 then k:=Integers()!(Integers()!MM[1,2] mod 4); w1:=v2; w2:=v1+(MM[1,2]-k)/4*v2; return [w1,w2],[w1,4*w2]; end if; return "Allo"; end function; /* Soit J:=ideal ou O=Z+wZ et I un ideal quelconque copremier a 3. Alors ce programme donne une base triangulee de (I,I*J^2) sous la forme: [w1,w2], [w1,4*w2] ou w1=a+bw avec b divisible par 3. */ special_basis_copremier_3:=function(f,K,O,w,I) D:=Discriminant(K); L:=QuadraticField(D); J:=ideal; B1,B2:=special_basis_2(f,K,O,w,I); w1:=B1[1]; w2:=B1[2]; a:=Integers()!Eltseq(w1)[1]; /* On ecrit w1 sous la forme w1=(a+bw) */ b:=Integers()!Eltseq(w1)[2]; c:=Integers()!Eltseq(w2)[1]; /* On ecrit w2 sous la forme w2=(c+dw) */ d:=Integers()!Eltseq(w2)[2]; if d eq 0 then w2:=w1+w2; c:=Integers()!Eltseq(w2)[1]; /* On ecrit w2 sous la forme w2=(c+dw) */ d:=Integers()!Eltseq(w2)[2]; end if; if (b mod 3) eq (0 mod 3) then return w1,w2; end if; if (d mod 3) eq (0 mod 3) then w2:=w1+w2; c:=Integers()!Eltseq(w2)[1]; /* On ecrit w2 sous la forme w2=(c+dw) */ d:=Integers()!Eltseq(w2)[2]; end if; pi:=homIntegers(f)|>; k:=Integers()!pi(b*d^(-1)); w1:=w1-4*k*w2; w2:=w2; return w1,w2; return "allo"; end function; /* Soit I un O_K-ideal de K=Q(w) ou w=(1+\sqrt{D})/2 ou D est un discriminant. Alors ce programme retourne un entier 1=< r<= f copremier a f et un element tau\in Q(w) tels que r\Lambda_{\tau} est equivalent a I et tel que \Lambda_{4\tau} est aussi un ideal Le programme retourne aussi un entier epsilon qui est egal a {+1,-1}. Ce programme retourne une base orientee i.e. tau-tau^{\sigma}>0 */ special_pair_2:=function(f,K,O,w,I) w1,w2:=special_basis_copremier_3(f,K,O,w,I); /* La base est donnee par rapport a {1,w} */ a:=Integers()!Eltseq(w1)[1]; /* On ecrit w1=a+b*w ou b=0 (mod 3) */ b:=Integers()!Eltseq(w1)[2]; epsilon:=Sign(Norm(w1)); tau:=w2/w1; r:=Integers()!(a mod f); if Sign(Eltseq(tau)[2]) eq -1 then tau:=-tau; end if; return r,tau,epsilon; end function; /* Meme chose que calculs_polynome_speed sauf que l'on essaie de corriger le probleme "can't find good representatives" Ce programme calcule le polynome minimal associe aux unites p-adiques construites pour un corps reel quadratique Q(\sqrt{D}) pour lequel p est inert et D=1 mod 8 et copremier a 3. Nos unites corresponde a l'unite modulaire associee au vecteur [2,-3,1]; */ calculs_polynome_speed_2:=function(p,f,N_0,D,delta,j,V1,V2,V3,V4,V5,V6,V7,precision) Qp:=pAdicField(p,precision); P:=PolynomialRing(Rationals()); K:=QuadraticField(D); KK:=NumberField(x^2-x+(1-D)/4); O:=MaximalOrder(KK); Kp:=UnramifiedExtension(Qp,x^2-D); R:=PolynomialRing(Kp); Real:=RealField(100); DD:=Real!D; G:=Gamma0(f*N_0); gen:=Generators(G); V1,V2,V3,V4:=swap_Qp_foot(f,Qp,V1,V2,V3,V4); I:=ideal; H,pi:=RayClassGroup(I,[1,2]); n:=#H; vecteur_ideaux:=[]; i:=1; for h in Set(H) do vecteur_ideaux[i]:=pi(h); i:=i+1; end for; vecteur_log_unites:=[]; vecteur_valuation:=[]; vecteur_racines_unite:=[]; vecteur_pgcd:=[]; vecteur_unites:=[]; for i in [1..n] do r,tau,epsilon:=special_pair_2(f,KK,O,w,vecteur_ideaux[i]); a:=Eltseq(tau)[1]; b:=Eltseq(tau)[2]; tau2:=a+b*(1+t)/2; /* On convertit w en (t+1)/2 */ M:=stabilisateur(f,N_0,D,tau2); /* Le stabilisateur */ tau2_p:=a+b*(1+tp)/2; /* tau_p est ecrit par rapport a la base tp^2=D */ if Norm(M[2,1]*tau2+M[2,2]) ne 1 then M:=M^(2); end if; if Norm(M[2,1]*tau2+M[2,2]) ne 1 then return "Il y a erreur"; end if; if M[2,1]*((b/2)*Real!SquareRoot(DD)+(a+b/2))+M[2,2] gt 1 then M:=M; else M:=M^(-1); end if; print([r,tau2]); print(endomorphism_order(K,tau2)); print(endomorphism_order(K,4*tau2)); print(multiply(M,[tau2,1])); print(M); print(M[2,1]*tau2+M[2,2]); printf "epsilon est egal a %o\n\n",epsilon; word:=FindWord(G,G!M); print(word); vecteur_log_unites[i]:=epsilon*integrale2_tau_speed(Kp, tp, tau2_p, G, gen, word, (Integers()!(j*r)) mod f, f, N_0, delta, V1, V2, V3, V4,V5,V6); sum:=vecteur_log_unites[i]; printf "vecteur log unites est calcule\n"; vecteur_valuation[i]:=epsilon*valuation(p, f, N_0, (Integers()!(j*r)) mod f, delta, G, gen, word); val:=vecteur_valuation[i]; printf "vecteur valuation est calcule\n"; m,root:=int_log_beta2_speed(Kp,tp,tau2_p,G,M,p,f,N_0,Integers()!( (j*r) mod f ),delta,V7); m:=epsilon*m; vecteur_racines_unite[i]:=m; printf "vecteur racine de l'unite est calcule\n"; d:=1; vecteur_unites[i]:=p^( Integers()!(val) )*root^( Integers()!( m mod (p^2-1) ) )*Exp(sum); end for; poly:=1; for i in [1..n] do poly:=poly*(y-vecteur_unites[i]); end for; return poly,vecteur_valuation,vecteur_racines_unite,vecteur_log_unites; end function; /* Ce programme construit un ensemble complet de representants de I_K(f)/P_{K,1}(f\infty) modulo le sous-groupe d'ordre 2 engendre par la conjugaison complexe. */ set_representatives:=function(D,f) P:=PolynomialRing(Rationals()); K:=QuadraticField(D); hom:=homIntegers(f)|>; G,pi:=ray_group(D,f); seq:=Setseq(Set(G)); n:=#G; for i in [1..n] do char,generator:=IsPrincipal(pi(seq[i])); if Sprint(char) eq "true" then aa:=generator[1]; /* On ecrit generator comme aa+bb(1+t)/2 */ bb:=generator[2]; a:=aa+bb/2; /* On ecrit generator comme a+bt ou t^2=D */ b:=bb/2; w:=a+b*t; unit:=FundamentalUnit(K); for j in [1..8] do ww:=unit^j*w; if ( hom(ww[1]) eq (1 mod f) or hom(ww[1]) eq (2 mod f) ) and hom(ww[2]) eq (0 mod f) and Sign(Norm(ww)) eq -1 then index:=i; end if; end for; end if; end for; S:=Set(G); rep:=[]; for i in [1..Integers()!(n/2)] do set_list:=Setseq(S); T:={set_list[1],seq[index]*set_list[1]}; S:=(S diff T); rep[i]:=set_list[1]; end for; return G,pi,rep; end function; /* Meme chose que calculs_polynome_speed sauf que l'on essaie de corriger le probleme "can't find good representatives" De plus, on essaie de diminuer le probleme d'ambiguite avec avec les racines p^2-1 de l'unite. Ce programme calcule le polynome minimal associe aux unites p-adiques construites pour un corps reel quadratique Q(\sqrt{D}) pour lequel p est inert et D=1 mod 8 et copremier a 3. Nos unites corresponde a l'unite modulaire associee au vecteur [2,-3,1]; */ calculs_polynome_speed_3:=function(p,f,N_0,D,delta,j,V1,V2,V3,V4,V5,V6,V7,precision) Qp:=pAdicField(p,precision); P:=PolynomialRing(Rationals()); K:=QuadraticField(D); KK:=NumberField(x^2-x+(1-D)/4); O:=MaximalOrder(KK); Kp:=UnramifiedExtension(Qp,x^2-D); R:=PolynomialRing(Kp); hom:=homIntegers(f)|>; Real:=RealField(100); DD:=Real!D; G:=Gamma0(f*N_0); gen:=Generators(G); V1,V2,V3,V4:=swap_Qp_foot(f,Qp,V1,V2,V3,V4); I:=ideal; H,pi:=RayClassGroup(I,[1,2]); n:=#H; seq:=Setseq(Set(H)); for i in [1..n] do char,generator:=IsPrincipal(pi(seq[i])); if Sprint(char) eq "true" then aa:=generator[1]; /* On ecrit generator comme aa+bb(1+t)/2 */ bb:=generator[2]; a:=aa+bb/2; /* On ecrit generator comme a+bt ou t^2=D */ b:=bb/2; Z:=a+b*t; unit:=FundamentalUnit(K); for j in [1..8] do ZZ:=unit^j*Z; if ( hom(ZZ[1]) eq (1 mod f) or hom(ZZ[1]) eq (2 mod f) ) and hom(ZZ[2]) eq (0 mod f) and Sign(Norm(ZZ)) eq -1 then index:=i; end if; end for; end if; end for; S:=Set(H); rep:=[]; for i in [1..Integers()!(n/2)] do set_list:=Setseq(S); T:={set_list[1],seq[index]*set_list[1]}; S:=(S diff T); rep[i]:=set_list[1]; end for; nn:=Integers()!(n/2); vecteur_ideaux:=[]; for i in [1..nn] do vecteur_ideaux[i]:=pi(rep[i]); end for; vecteur_log_unites:=[]; vecteur_valuation:=[]; vecteur_racines_unite:=[]; vecteur_pgcd:=[]; vecteur_unites:=[]; for i in [1..nn] do r,tau,epsilon:=special_pair_2(f,KK,O,w,vecteur_ideaux[i]); a:=Eltseq(tau)[1]; b:=Eltseq(tau)[2]; tau2:=a+b*(1+t)/2; /* On convertit w en (t+1)/2 */ M:=stabilisateur(f,N_0,D,tau2); /* Le stabilisateur */ tau2_p:=a+b*(1+tp)/2; /* tau_p est ecrit par rapport a la base tp^2=D */ if Norm(M[2,1]*tau2+M[2,2]) ne 1 then M:=M^(2); end if; if Norm(M[2,1]*tau2+M[2,2]) ne 1 then return "Il y a erreur"; end if; if M[2,1]*((b/2)*Real!SquareRoot(DD)+(a+b/2))+M[2,2] gt 1 then M:=M; else M:=M^(-1); end if; print([r,tau2]); print(endomorphism_order(K,tau2)); print(endomorphism_order(K,4*tau2)); print(multiply(M,[tau2,1])); print(M); print(M[2,1]*tau2+M[2,2]); printf "epsilon est egal a %o\n\n",epsilon; word:=FindWord(G,G!M); print(word); vecteur_log_unites[i]:=epsilon*integrale2_tau_speed(Kp, tp, tau2_p, G, gen, word, (Integers()!(j*r)) mod f, f, N_0, delta, V1, V2, V3, V4,V5,V6); sum:=vecteur_log_unites[i]; printf "vecteur log unites est calcule\n"; vecteur_valuation[i]:=epsilon*valuation(p, f, N_0, (Integers()!(j*r)) mod f, delta, G, gen, word); val:=vecteur_valuation[i]; printf "vecteur valuation est calcule\n"; m,root:=int_log_beta2_speed(Kp,tp,tau2_p,G,M,p,f,N_0,Integers()!( (j*r) mod f ),delta,V7); m:=epsilon*m; vecteur_racines_unite[i]:=m; printf "vecteur racine de l'unite est calcule\n"; d:=1; vecteur_unites[i]:=p^( Integers()!(val) )*root^( Integers()!( m mod (p^2-1) ) )*Exp(sum); end for; for i in [1..nn] do vecteur_log_unites[nn+i]:=-vecteur_log_unites[i]; vecteur_valuation[nn+i]:=-vecteur_valuation[i]; vecteur_racines_unite[nn+i]:=-vecteur_racines_unite[i]; vecteur_unites[nn+i]:=vecteur_unites[i]^(-1); end for; poly:=1; for i in [1..n] do poly:=poly*(y-vecteur_unites[i]); end for; return poly,vecteur_valuation,vecteur_racines_unite,vecteur_log_unites; end function; /* Ce programme determine si K(f\infty), ou K=\QQ(\sqrt{D}), est un CM field. */ CM_field:=function(D,f) P:=PolynomialRing(Rationals()); K:=QuadraticField(D); hom:=homIntegers(f)|>; unit:=FundamentalUnit(K); if Norm(unit) eq -1 then unit:=unit^2; /* Il faut choisir la plus petite unite totallement positive */ end if; for j in [1..8] do w:=Eltseq(unit^j); a:=w[1]; /* On ecrit unit^j=a+b\sqrt{D} */ b:=w[2]; if hom(a) eq (-1 mod f) and hom(b) eq (0 mod f) then return "yes",j,unit; end if; end for; return "no"; end function;