Q1. Une équation fonctionnelle 

>	a[0]:=1+x;b[0]:=1-x;a[1]:=(a[0]+b[0])/2;b[1]:=sqrt(a[0]*b[0]);

 

>	eq:=a[0]*M(1,b[0]/a[0])=M(a[1],b[1]);

 

Q2. Les premiers coefficients de  

>	eval(eq,M=proc(toto,u) 1/A(1-u^2) end);

 

>	map(x->1/x,%);

 

>	eqA:=op(1,%)-op(2,%);

 

>	series(eval(eqA,A=proc(t) 1+add(u[i]*t^i,i=1..10) end),x,10);

 

>	solve({seq(coeff(%,x,i),i=0..9)},{seq(u[i],i=1..10)});

 

>	sol:=subs(%,series(1+add(u[i]*x^i,i=1..9),x,10));

 

Au passage, ce calcul montre que l'équation a au plus une solution formelle avec A(0)=1 car le système est triangulaire. Comme l'énoncé annonce que M est analytique, ce développement est celui de A. 

Q3. Une équation différentielle satisfaite par cette série tronquée 

>	with(gfun):

 

>	seriestodiffeq(series(sol,x,11),y(x));

 

>	deq:=op(1,%):

 

>	diffeqtorec(deq,y(x),u(n));

 

Le terme de tête de la récurrence ne s'annule pas sur ℕ, l'équation différentielle admet donc une unique solution série avec ces conditions initiales. 

Q4. La preuve que (et pas seulement sa troncature) est solution de l'équation différentielle 

La preuve consiste à montrer que cette solution de l'équation différentielle vérifie aussi l'équation fonctionnelle.
Pour cela, on construit une équation différentielle satisfaite par  

>	F(x)=map(normal,subs(A=y,eqA));

 

lorsque est solution de l'équation différentielle. On peut commencer, en raisonnant juste sur les ordres, par constater que F satisfait une équation différentielle d'ordre au plus 4. Il serait alors tentant de conclure que la vérification des 4 premiers coefficients conclut la preuve, mais il faut d'abord s'assurer que le coefficient de tête de l'équation différentielle ne s'annule pas en 0. 

>	deq1:=subs(y=y1,algebraicsubs(deq,y-x^2,y(x)));

 

>	deq2:=subs(y=y2,algebraicsubs(deq,numer(y-(4*x/(1+x)^2)),y(x)));

 

>	poltodiffeq(y2(x)/(1+x)-y1(x),[deq1,deq2],[y1(x),y2(x)],F(x));

 

Cette équation différentielle n'a que 0 comme solution, ce qui est la preuve cherchée. 

Q5. Le lien avec la série hypergéométrique 

>	diffeqtorec(deq,y(x),u(n));

 

>	collect(%,u,factor);

 

Cette récurrence est exactement celle que satisfont les coefficients de la série , qui est donc égale à . 

Q6. Une première récurrence 

On peut commencer par la deviner : 

>	recd:=listtorec([seq(int(t^(2*n)/sqrt(1-t^2),t=0..1),n=0..10)],d(n))[1];

 

Ensuite, on la prouve par intégration par parties : 

>	B:=Int(t^(2*n)/sqrt(1-t^2),t=0..1);

 

>	student[intparts](B,t^(2*n-1)) assuming n>1;

 

on lit ainsi d(n)=(2n-1)(d(n-1)-d(n)). 

Q7. Les représentations hypergéométriques des intégrales elliptiques complètes 

On peut faire d'une pierre deux coups en introduisant 

>	ff:=(1-k^2)^alpha;

 

>	deq:=diff(y(k),k)=normal(diff(ff,k)/ff)*y(k);

 

On développe ensuite en k et on intègre terme à terme. Les coefficients du développement sont donnés par 

>	recbis:=diffeqtorec({deq,y(0)=1},y(k),d(n));

 

et ceux de l'intégrale par : 

>	recKE:=`rec*rec`(recd,recbis,d(n));

 

>	recKE:=collect(recKE,d,factor);

 

On récupère ainsi les deux récurrences cherchées : 

>	recK:=subs(alpha=-1/2,recKE);recE:=subs(alpha=1/2,recKE);

 

Q8. Une relation entre et  

est définie par une équation différentielle que nous avons déjà calculé : 

>	deqF:=holexprtodiffeq(hypergeom([1/2,1/2],[1],x),y(x));

 

>	poltodiffeq((1-x)*(2*x*diff(y(x),x)+y(x)),[deqF],[y(x)],y(x));

 

>	diffeqtorec(%,y(x),u(n));

 

On retrouve bien l'équation satisfaite par les coefficients de la série , avec la même condition initiale : 

>	diffeqtorec(holexprtodiffeq(hypergeom([-1/2,1/2],[1],x),y(x)),y(x),u(n));

 

Q9. Le membre gauche de la relation de Legendre 

>	G:=x-> (1-x)*(2*x*D(F)(x)+F(x));

 

>	G(x)*F(1-x)+G(1-x)*F(x)-F(x)*(F(1-x));

 

>	factor(%);

 

>

G:='G':

 

Q10. et satisfont la même équation 

>	algebraicsubs(deqF,y-(1-x),y(x));

 

Il ne reste qu'à calculer ce Wronskien : 

>	dsolve(diff(w(x),x)=-(-8*x+4)/(-4*x^2+4*x)*w(x),w(x));

 

Le reste du calcul se termine sans l'aide du calcul formel. 

Q13. Les 60 premières décimales de π 

>	Digits:=60:

 

>	a[0]:=1;b[0]:=1/sqrt(2.);c[0]:=sqrt(a[0]^2-b[0]^2);

 

>	for i to 7 do a[i]:=(a[i-1]+b[i-1])/2;b[i]:=sqrt(a[i-1]*b[i-1]); c[i]:=c[i-1]^2/4/a[i] od;

 

>	for j to 6 do j,2*a[j]^2/(1-add(2^k*c[k]^2,k=0..j-1)) od;

 

Q14. Le calcul en entiers 

>	decPi:=proc(D) local a, b, c, i; a[0]:=10^D:b[0]:=isqrt(iquo(10^(2*D),2)):c[0]:=isqrt(a[0]^2-b[0]^2):for i while c[i-1]<>0 do a[i]:=iquo(a[i-1]+b[i-1],2);b[i]:=isqrt(a[i-1]*b[i-1]); c[i]:=iquo(c[i-1]^2,4*a[i]) od;iquo(2*a[i-1]^2*a[0],(a[0]^2-add(2^k*(a[k]^2-b[k]^2),k=0..i-2))) end;

 

>	decPi(10^4);

 

>	evalf(%-trunc(10^(10^4)*Pi),10^4);

 

>	time(decPi(10^4)),time(decPi(2*10^4)),time(decPi(4*10^4));

 

>	time(decPi(8*10^4)),time(decPi(16*10^4));

 

Le temps est multiplié par un peu plus que 2, mais moins que 3 à chaque étape, indiquant que la multiplication est effectuée par FFT, et que la complexité est bonne. On obtient ainsi un million de décimales en moins d'une minute : 

>	time(decPi(10^6));