Q1. Calcul des valeurs numériques 

>	listvals:=[seq(Sum(1/16^k/(8*k+i),k=0..infinity),i=1..8)];

 

(1.1.1)

 

>	Lvals:=evalf(listvals,100);

 

(1.1.2)

 

Q2. La découverte 

>	IntegerRelations[LinearDependency]([evalf(Pi,100),op(Lvals)],method=LLL);

 

(1.2.1)

 

Q3. Une procédure de recherche d'identités 

>	test:=proc(cte,listvals,prec) local lnum,rel; lnum:=evalf([cte,op(listvals)],prec);rel:=IntegerRelations[LinearDependency](lnum,method=LLL); if rel[1]<>0 then cte=add(-rel[i+1]/rel[1]*listvals[i],i=1..nops(listvals)) else 0=add(rel[i+1]*listvals[i],i=1..nops(listvals)) fi end;

 

(1.3.1)

 

>	test(Pi,listvals,100);

 

(1.3.2)

 

Q4. Des identités pour ln(2), ln(3),... 

>	for val in [ln(2),ln(3),ln(5),arctan(2),arctan(3),sqrt(2)*arctan(1/sqrt(2)),sqrt(2)*ln(1+sqrt(2))] do test(val,listvals,100) od;

 

 

 

 

 

 

 

	(1.4.1)

 

Q5. Autres constantes, autres identités 

>	L:=[seq(seq(Sum(2^(-j*i)/i^m,i=1..infinity),j=1..6),m=1..5)];

 

(1.5.1)

 

>	test(Pi^2,L,100);

 

(1.5.2)

 

Cette identité pour 0 gêne la suite de la recherche. Il suffit pour s'en débarasser d'enlever un des éléments liés par cette relation. 

>	L:=subsop(6=NULL,L):

 

>	test(Pi^2,L,100);

 

(1.5.3)

 

>	test(ln(7),L,100);

 

(1.5.4)

 

>	test(ln(2)^2,L,100);

 

(1.5.5)

 

Q6. Les sommes  

>	Sum(t^(8*i+j)/(8*i+j),i=0..infinity);

 

(2.1.1)

 

>	diff(%,t);

 

(2.1.2)

 

>	diffS:=value(%);

 

(2.1.3)

 

>	for jj to 8 do S[jj]:=int(subs(j=jj,diffS),t=0..z) assuming 0<z,z<1 od;

 

 

 

 

 

 

 

 

	(2.1.4)

 

Q7. Les sommes et la preuve 

>	for j to 8 do Sigma[j]:=simplify(eval(S[j]/z^j,z=(1/16)^(1/8))) od;

 

 

 

 

 

 

 

 

	(2.2.1)

 

vérification 

>	for j to 8 do evalf(Sigma[j]-listvals[j],100) od;

 

 

 

 

 

 

 

 

	(2.2.2)

 

>	4*Sigma[1]-2*Sigma[4]-Sigma[5]-Sigma[6];

 

(2.2.3)

 

>	convert(%,ln);

 

(2.2.4)

 

>	combine(%,ln);

 

(2.2.5)

 

Q8. Exponentiation binaire 

>	expbin:=proc(m,k) local q,r,rec; if m=0 then 1 else q:=iquo(m,2,r); rec:=expbin(q,k)^2 mod k; if r=0 then rec else 16*rec mod k fi fi end;

 

>

 

(3.1.1)

 

>	seq(expbin(i,37),i=1..10);

 

(3.1.2)

 

>	seq(16^i mod 37,i=1..10);

 

(3.1.3)

 

L'efficacité de expbin est handicapée par le fait que l'execution mène à de nombreux calculs sur des petits entiers, et donc coûte cher à un langage interprêté comme Maple. Il vaut mieux recourir à une version écrite dans le noyau (donc en C), pour la suite : 

>	seq(`&^`(16,i) mod 37,i=1..10);

 

(3.1.4)

 

Q9. Les sommes 

>	S1:=proc(N,q,p,k) local s, i, modulo, res; if q=0 then 0 else s:=1./eval(p,k=N);for i from 0 to N-1 do modulo:=eval(p,k=i); s:=s+(`&^`(16,N-i) mod modulo)/modulo od; res:=frac(q*s); if res<0 then 1+res else res fi fi end;

 

(3.2.1)

 

>	S2:=proc(N,q,p,k) local s,i; evalf(Sum(q/16^(i-N)/subs(k=i,p),i=N+1..infinity)) end;

 

(3.2.2)

 

vérification 

>	S1(0,1,8*n+1,n);

 

(3.2.3)

 

>	S2(0,1,8*n+1,n);

 

(3.2.4)

 

>	listvals[1]-%-%%;

 

(3.2.5)

 

>

evalf(%);

 

(3.2.6)

 

>	S1(0,1,8*n+2,n);

 

(3.2.7)

 

>	S2(0,1,8*n+2,n);

 

(3.2.8)

 

>	evalf(listvals[2]-%-%%);

 

(3.2.9)

 

Q10. La procédure finale 

>	somme:=proc(N,B) local j, k; frac(add(S1(N,B[j],8*k+j,k)+S2(N,B[j],8*k+j,k),j=1..8)) end;

 

(3.3.1)

 

>	listpi:=[4,0,0,-2,-1,-1,0,0];

 

(3.3.2)

 

>	somme(0,listpi);

 

(3.3.3)

 

>	evalf(Pi-%);

 

(3.3.4)

 

>	Digits:=100:

 

>	evalf(Pi-somme(0,listpi));

 

(3.3.5)

 

Donc ça fonctionne bien quand on part du tout début.
Voici la valeur du reste de π à partir du 10ième chiffre en base 16 : 

>	evalf(frac(16^10*Pi));

 

(3.3.6)

 

>	somme(10,listpi);

 

(3.3.7)

 

>

%-%%;

 

(3.3.8)

 

>	sommefinale:=proc(N,B,k) local bk; Digits:=k*(ilog10(16)+1)+max(ilog2(N),1);bk:=somme(N,B); convert(trunc(16^k*bk),base,16) end;

 

(3.3.9)

 

Les 20 premières à partir de 0 (il faut lire en partant de la fin) : 

>	sommefinale(0,listpi,20);

 

(3.3.10)

 

les 10 premières : 

>	sommefinale(0,listpi,10);

 

(3.3.11)

 

les 10 suivantes : 

>	sommefinale(10,listpi,10);

 

(3.3.12)

 

>	%%%-[op(%),op(%%)];

 

(3.3.13)

 

Quelques tests de temps : les 5 premières décimales à partir de la 1e, 10e, 100e,... 

>	for i to 6 do st:=time(); print(i,sommefinale(10^i,listpi,5),time()-st) od:

 

 

 

 

 

 

	(3.3.14)

 

À la fin, en un peu plus de 30 secondes de calcul, des décimales au-delà de la  millionième ont été obtenues. 

Pour le calcul de π, où les dénominateurs sont des , il est possible d'améliorer un peu les temps de calculs en utilisant add au lieu d'une boucle : 

>	S1:=proc(N,q,p,k) local s, i, modulo, ini, res; if q = 0 then 0 else ini:=eval(p,k=0);s:=add(1.*(`&^`(16,N-i) mod (8*i+ini))/(8*i+ini),i=0..N);res:=frac(q*s); if res<0 then res+1 else res fi fi end;

 

(3.3.15)

 

>	st:=time():sommefinale(10^6,listpi,5),time()-st;

 

(3.3.16)

 

>

 

Q11. Complexité 

La complexité est approximativement linéaire en , et elle est reflétée par les temps de calcul qui doublent à peu près en doublant . Plus précisément, chaque exponentiation binaire prend étapes, chaque étape étant une opération modulo un entier de taille . Au total, le temps de calcul est donc dominé par , où est le temps utilisé pour multiplier deux entiers de taille C'est donc du avec une multiplication naïve, ou du avec de la FFT.  

Ainsi, cette méthode n'est pas théoriquement meilleure que les meilleurs algorithmes pour à base de moyenne arithmético-géométrique et de FFT, qui calculent toutes les premières décimales. 

Par contre, la consommation mémoire de la méthode présentée ici est négligeable, et une multiplication rapide n'est pas nécessaire pour obtenir des résultats.