L’algorithme d’Euclide étendu consiste à calculer la suite de relations
où
dont on tire la relation de Bezout
en prenant \(n\) tel que \(r_n>0\) et \(r_{n+1}=0\).
iquo et irem.Sur \(\k[x]\), les fonctions de quotient et reste deviennent quo et rem.
En outre, le pgcd n’est défini qu’à un inversible près, et il est préférable de
choisir à chaque étape \(r_k\) unitaire.
L’algorithme s’écrit alors
euclideetendu(a,b):={
local la,lb,uk1,uk,vk1,vk,rk1,rk,qk,r,l;
la:=lcoeff(a); lb:=lcoeff(b);
uk1:=1/la; vk1:=0; rk1:=a/la;
uk:=0; vk:=1/lb; rk:=b/lb;
while(rk!=0) {
qk := normal(quo(rk1,rk));
r := normal(rk1 - qk * rk);
l := lcoeff(r);
if(r==0) {l:=1;}
uk1, uk := uk, (uk1 - qk * uk)/l;
vk1, vk := vk, (vk1 - qk * vk)/l;
rk1, rk := rk, r/l;
};
return(normal([uk1,vk1,rk1]));
}:;
Exécuter cet algorithme sur des polynômes entiers à coefficients aléatoires.
Que se passe-t-il ? Comparer les résultats avec la fonction xcas gcdex.
On remarque que les coefficients rationnels explosent. De son côté xcas renvoie un pgcd dans \(\Z[X]\), mais les résultats sont bien associés
a:=randPoly(x,5):; b:=randPoly(x,5):;
A := euclideetendu(a,b);
B := gcdex(a,b);
normal(A*B[2])==B;
Modifier l’algorithme en ajoutant un troisième argument \(t\) qui fasse retourner le triplet \((u_k,v_k,r_k)\) au dernier rang \(k\) tel que \(\deg(v_k)\leq t\) (avec un déroulement normal si \(t>=deg(a)\)).
euclideetendu(a,b,t):={
local la,lb,uk1,uk,vk1,vk,rk1,rk,qk,r,l;
la:=lcoeff(a); lb:=lcoeff(b);
uk1:=1/la; vk1:=0; rk1:=a/la;
uk:=0; vk:=1/lb; rk:=b/lb;
while(rk!=0 and degree(vk)<=t) {
qk := normal(quo(rk1,rk));
r := normal(rk1 - qk * rk);
l := lcoeff(r);
if(r==0) {l:=1;}
uk1, uk := uk, (uk1 - qk * uk)/l;
vk1, vk := vk, (vk1 - qk * vk)/l;
rk1, rk := rk, r/l;
//print( normal([uk,vk,rk]), normal(rk/vk) );
};
return(normal([uk1,vk1,rk1]));
}:;
Soit \(b(x)=\sum_{k\geq 0} b_k x^k\) une série entière. L’approximant de Padé d’ordre \((m,n)\) de \(b\) est la fraction rationnelle \(r/v\) vérifiant
avec \(r,v\) de degrés \(m\) et \(n\).
Calculer le développement en série entière \(b(x)\) de \(\tan(x)\) au voisinage de \(0\) à l’ordre 9, et le convertir en polynôme (de degré 9).
b := series(tan(x),x,0,10);
b := convert(b,polynom);
Représenter graphiquement \(b(x)\) et \(\tan(x)\) sur \(]-3/2,2[\).
plot([b,tan(x)],x=-3/2..2);
Écrire l’algorithme d’Euclide étendu entre \(x^{10}\) et \(b(x)\), et récupérer tous les couples \(v_k,r_k\) issus de l’algorithme.
v := euclideetendu(x^10,b,1):; normal(v[2]/v[1]);
v := euclideetendu(x^10,b,2):; normal(v[2]/v[1]);
v := euclideetendu(x^10,b,3):; normal(v[2]/v[1]);
v := euclideetendu(x^10,b,4):; normal(v[2]/v[1]);
v := euclideetendu(x^10,b,5):; normal(v[2]/v[1]);
v := euclideetendu(x^10,b,6):; normal(v[2]/v[1]);
v := euclideetendu(x^10,b,7):; normal(v[2]/v[1]);
v := euclideetendu(x^10,b,8):; normal(v[2]/v[1]);
Pour chaque fraction \(r_k/v_k\), tracer la courbe représentative. Que dire de la qualité d’approximation de \(\tan(x)\), comparé à la formule de Taylor ?
Écrire une fonction approx_pade(b,n,m) qui calcule le \((n,m)\)
approximant de Padé de la série \(b\). Comparer avec la fonction
correspondante xcas.
pade(tan(x),x,9,3);
pade(tan(x),x,9,5);
pade(tan(x),x,9,6);
Soit \((b_k)_{k\in\N}\) une suite récurrente linéaire d’ordre \(n\), c’est-à-dire qu’il existe un polynôme \(c(x)=x^n-\sum_{i=0}^{n-1}c_ix^i\) tel que
Écrire une fonction linrec(b,c,k) qui à partir des \(n\) premiers
coefficients de la suite et des coefficients \(c_0,\dots c_{n-1}\) du polynôme
annulateur calcule la liste des \(k\) premiers termes.
linrec(b,c,k):= {
local n,l,j;
l := length(c);
for(n:=length(b);n<k;n++) {
b := append(b,sum(c[j]*b[n-l+j],j,0,l-1));
}
return b;
}:;
Les suites récurrentes à coefficients dans des corps finis sont périodiques, pourquoi ? Trouver la période des suites récurrentes sur \(\F_2\) données par les valeurs initiales \(1,1,\dots 1\) et les polynômes
linrec([seq(1%2,5)],[0,1,0,0,1],40)%0
linrec([seq(1%2,5)],[1,0,0,1,0],40)%0
On considère la suite de valeurs initiales \(1,2,3,4,5\) et de polynôme générateur \(x^5-x^4+x^3-2x^2+2x-1\). Former la série \(b(x)=\sum_k b_kx^x\) tronquée à l’ordre \(10\).
b := linrec([1,2,3,4,5],[1,-2,2,-1,1],10)
P := sum(b[k]*x^k,k,0,10)
À l’aide de la fonction pade, retrouver le polynôme générateur à partir
des premiers termes de la suite.
Cet algorithme de résolution des systèmes linéaires \(Ax=b\) est particulièrement intéressant dans le cas de matrices sur les corps finis, de grandes dimension et «creuses», c’est-à-dire avec peu de coefficients non-nuls.
C’est le cas par exemple pour l’algorithme page-rank de google, ou pour les algorithmes de factorisation par crible (matrices de dimension plusieurs millions, mais avec moins d’une centaine de coefficients non-nuls par ligne).
Dans ce genre de cas, on ne stocke pas toute la matrice mais seulement l’emplacement de ses coefficients non-nuls. Inverser \(A\) par une technique de pivot de Gauss est hors de question, puisque le pivot a tendance à remplir la matrice.
La technique de Wiedemann consiste à calculer le polynôme minimal de l’endomorphisme \(A\) restreint au sous-espace de Krylov engendré par les itérés de \(b\) par \(A\). Ce polynôme permet alors d’exprimer \(A^{-1}b\) comme un polynôme en les \(A^kb\).
En réalité, on ne trouve pas le polynôme minimal de \(A\) mais un annulateur de la trace (relativement à la forme linéaire \(u^T\)) du minimal de \(A\) sur le sous-espace engendré par les itérés de \(v\) (appelé sous-espace de Krylov).
Une approche probabiliste consiste à prendre le ppcm des annulateurs trouvés via divers couples \((u,v)\) jusqu’à vérifier que l’on annule \(A\). Pour montrer que c’est efficace, on peut montrer des résultats de densité sur les vecteurs \(u\) et \(v\) qui fournissent le résultat.
L’algorithme de Wiedemann ne s’intéresse en réalité pas au polynôme minimal de \(A\), mais bien au minimal sur l’espace de Krylov des \(A^ku\), le but étant d’inverser \(A\)