Ok ok...
Gros spoiler, donc. La fin n'est pas très amusante et franchement pas nécessaire pour comprendre le cœur de l'algo.
Bonne lecture (déconseillée aux moins de 18 ans) !
On remarque qu'une fois que l'on a choisi les applications effectuées sur la première ligne, on peut déduire celles qui
sont effectuées ensuite. En effet, un 1 sur une ligne où l'on ne peut plus effectuer d'application ne peut être formé
qu'en effectuant une application sur la même colonne, une ligne en-dessous.
On peut écrire cela matriciellement :
A=
(T I) T : matrice tridiagonale remplie de 1
(I 0) I : matrice identité 0 : matrice nulle.
Multiplier cette matrice par un vecteur (u v) correspond à effectuer les applications nécessaires sur u, et renvoyez
les deux lignes suivantes, la ligne après u étant v, la ligne après v étant nulle.
On peut ramener le problème initial au même problème avec une image nulle, sauf sur la dernière ligne dont les valeurs
sont les coordonnées du vecteur C et où on n'effectue que des applications sur la première ligne, le reste se
déduisant.
Il faut donc résoudre A\^(m-1)*(x 0)=(C y).
Pour cela, on introduit B tel que B corresponde aux n premières lignes et colonnes de A\^(m-1).
On doit alors résoudre Bx=C. Cela permet de coder un algorithme en O(M(n)*log m) M(n) étant le coût d'une
multiplication de matrice nxn (O(n³) en général).
Soit P(-1)=0,P(0)=1,P(n+1)=P(n)X-P(n-1) (avec P(n) € (Z/2Z)[X])
On peut prouver que A\^m=
(P(m)(T) P(m-1)(T)) T étant la même matrice tridiagonale.
(P(m-1)(T) P(m-2)(T))
Le polynôme caractéristique chi de T est P(n) o (x-1). Par le théorème de Cayley-Hamilton, il est annulateur.
Pour obtenir B, on peut donc calculer (P(m-1) % chi)(T).
On a alors un algorithme en O(n³/log n+nm/log n).
À partir d'ici, on considère que B est inversible et qu'il faut retrouver x (c'est plus simple à expliquer, plus sympa
et il n'est pas très compliqué de répondre au problème d'origine une fois que l'on a compris comme faire ça).
Cherchons un polynôme annulateur non nul q de B.
Pour cela, on calcule vB\^ku pour kSoit ai tel que q=somme pour tout i de ai*x\^i
On a donc somme aiB\^i=0 et somme ai*(vB\^ku)=0, ce qui veut dire que la suite des vB\^ku est une suite récurrente
linéaire à coefficient constant pour tout u et v.
L'algorithme de Berlekamp-Massey (il s'agit essentiellement d'une algorithme d'Euclide étendu sur polynômes) permet,
connaissant les vB\^ku, de connaître le polynôme annulateur de la suite vB\^ku. Avec un peu de chance, celui-ci est
annulateur de B. On prend donc les vecteurs u et v au hasard.
On a alors q(0)=1 et q(B)=0.
On décompose q en 1+x*r.
I+B*r(B)=0 donc I=B*r(B) et B\^-1=r(B).
Enfin x=r(B)C et pour calculer x, il faut calculer des B\^kC.
Il ne reste plus que ce problème : calculer efficacement des vecteurs ayant un rapport avec B\^ku pour 0Soit K un
entier.
On calcule P\^k(m) % chi en effectuant O(K) multiplications (la réduction se fait en effectuant une multiplication par
l'inverse).
On veut calculer maintenant vB\^ku pour tout 0Clairement, vB\^ku=vB\^(k/K)B\^(k%K)u.
On a l'algo suivant qui répond à cette question :
pour chaque k/K
calculer u'=B\^(k/K)u
calculer w, vecteur dont la coordonnée i est vT\^iu', pour 0 pour chaque chaque k%K
vB\^ku=(p\^(k%K) % chi)*w (produit scalaire, le polynôme étant considéré comme un vecteur)
De même, on veut calculer x=somme aiB\^iu.
x=0
pour chaque k/K
calculer u'=B\^(k/K)u
calculer T\^iu' pour 0 v=0
pour chaque chaque k%K
si ai=1
v+=p\^(k%K) % chi (le polynôme étant considéré comme un vecteur)
pour chaque i avec vi=1
x+=T\^iu'
Et voilà :)
Soit MM(n) le temps mis pour effectuer le produit de 2 polynômes de degré n sur Z/2Z[X].
Le temps total est de O(KMM(n)+n/K*n²/log n+nm/log n).
On dérive en fonction de K : MM(n)-n³/log n/K²
On choisit donc K=sqrt(n³/log n/MM(n)).
Une multiplication classique de polynôme se fait O(n²/log n), que l'on peut accélérer facilement en O(n²/log² n).
Cela donne un algorithme en O(n²/log²n*sqrt(n/log n)+nm/log n).
On peut aussi utiliser Karastuba, qui donne des résultats pas mauvais :
(ax\^n+b)(cx\^n+d)=acx\^2n+((a+b)(c+d)-ac-bd)x\^n+bd que l'on applique récursivement.
Enfin, on peut faire une FFT dans l'anneau Z/2Z[X]/(X\^2n+X\^n+1).
Les diviseurs de 0 sont les puissances de x²+x+1
On a :
- pour tout 0 0
- x\^2n=x\^n+1
- pour tout 2n0
- x\^3n=x\^n*x\^2n=x\^n(x\^n+1)=x\^2n+x\^n=1
- x\^(1+2k)=x mod x²+x+1
- x\^(2+2k)=x+1 mod x²+x+1
x est donc une racine primitive d'ordre 3n
Soit f(x)=a0+a1x+a2x\^2...
On cherche F(k)=f(omega\^k)=f(x\^k)=a0+a1x\^k+a2x\^2k...
On définit g(x)=a0+a3x\^k+a6x\^2k..., h(x)=a1+a4x\^k..., i(x)=a2+a5x\^k... et G(k), H(k), I(k) comme F(k).
On a alors
F(k)=a0+a1x\^k+a2\^2k...=(a0+a3x\^3k...)+(a1x\^k+a4x\^4k...)+(a2x\^2k+a5x\^5k...)=(a0+a3x\^3k...)+x\^k(a1+a4x\^3k...)+x\^2k(a2+a5x\^3k...)=G(3k)+x\^kH(3k)+x\^2kI(3k)
Enfin, on a x\^3n=1 donc G(3(k+n))=G(3k+3n)=a0+a3x\^3kx\^3n...=a0+a3x\^3k...=G(3k).
x*x\^(3n-1)=1 donc x\^(3n-1) permet d'obtenir la FFT inverse.
De plus, on a x\^3n=1 donc x\^(3n+k)=x\^k donc on peut calculer l'exposant modulo 3n.
Pour que l'algorithme marche, il y a plusieurs conditions :
- x soit une racine primitive d'ordre 3n pour effectuer une FFT sur 3n coefficients
- les coefficients du produit doivent être de degré - le degré du produit doit lui-même être - 3n soit une puissance de
3
Cela implique que les coefficients des polynômes devant être multipliés soit de degré
L'algorithme pour multiplier 2 polynômes p,q de degré n est donc le suivant :
calculer FFT(p),FFT(q)
calculer le produit de convolution FFT(p)*FFT(q) en utilisant le même algorithme
retourner FFTinv(FFT(p)*FFT(q))
La complexité est C(n)=O(n)+sqrt(n)*C(sqrt(n))=n log log n (O(n) par étage, O(log log n) étages)
K=sqrt(n³/log n/(n log log n))=n/sqrt(log n*log log n)
D'où la complexité finale O(n²sqrt(log log n/log n)+nm/log n).
Sauf que l'on a oublié la lecture des nm entrées, chacune se fait en O(1) et donc on a un temps linéaire.
Note : ça peut se généraliser et résoudre rapidement toute équation du type p(A)x=b où A est creuse.
Exo (pour voir si vous avez bien suivi :) ) : retrouvez l'algorithme de Schönnage-Strassen, l'algorithme de Wiedemann et
l'algorithme pas de bébé-pas de géant dans l'explication.
Code (pour l'algo en O(M(N)+nm/log n)) : http://pastebin.com/mWYURBrh (3 µs/matrice 15x100)
Exemple pour :
0 0 0
0 0 1
1 0 0
n=3, m=3
On appuie sur la dernière ligne, dernière colonne et ça donne comme dernière ligne :
1 1 1
A=
1 1 0 1 0 0
1 1 1 0 1 0
0 1 1 0 0 1
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
A\^3=
0 1 1 1 0 1
1 1 1 0 0 0
1 1 0 1 0 1
1 0 1 1 1 0
0 0 0 1 1 1
1 0 1 0 1 1
Donc B=
0 1 1
1 1 1
1 1 0
On cherche à résoudre
Bx=(1 1 1)
x=(0 1 0) est la seule solution.
Pour aller vite :
P(-1)=0
P(0)=1
P(1)=x
P(2)=x²+1
P(3)=x³
chi=P(3) o (x-1)=(x+1)³=x³+x²+x+1
On peut vérifier que T³+T²+T+T\^0=0
T=
1 1 0
1 1 1
0 1 1
P¹(3)=x³=x²+x+1 : on peut vérifier que B¹=T²+T+1
P²(3)=(x²+x+1)*(x²+x+1)=x⁴+x²+x=x² mod chi : on peut vérifier que B²=T²
P³(3)=x²*(x²+x+1)=x⁴+x³+x²=x mod chi : on peut vérifier que B³=T
P⁴(3)=x*(x²+x+1)=x³+x²+x=1 mod chi : on peut vérifier que B⁴=T⁰=I
P⁵(3)=1*(x²+x+1)=x²+x+1 mod chi : on peut vérifier que B⁵=T²+T+T⁰
P⁶(3)=(x²+x+1)*(x²+x+1)=x² mod chi : on peut vérifier que B⁶=T²
On choisit au hasard u=(0 0 1), v=(1 0 0).
T⁰u=(0 0 1) vT⁰u=0
T¹u=(0 1 1) vT¹u=0
T²u=(1 0 0) vT²u=1
vB⁰u=vT⁰u=0
vB¹u=vT²u+vT¹u+vT⁰u=1
vB²u=vT²u=1
vB³u=vT¹u=0
vB⁴u=vT⁰u=0
vB⁵u=vT²u+vT¹u+vT⁰u=1
vB⁶u=vT²u=1
On a donc la suite récurrente linéaire à coefficients constants 0 1 1 0 0 1 1, on cherche les coefficients.
u(n)=u(n-1)+u(n-2)+u(n-3)+u(n-4)
et u(-3)=1 u(-2)=1 u(-1)=0
Donc on peut espérer que B³+B²+B¹+B⁰=0.
D'où B⁻¹=B²+B¹+B⁰
La dernière ligne était u=(1 1 1)
T⁰u=(1 1 1)
T¹u=(0 1 0)
T²u=(1 1 1)
B⁰=T⁰
B¹=T²+T¹+T⁰
B²=T²
Donc B⁻¹=T¹
On calcule x=T¹u=(0 1 0).
Matrice des applications :
0 1 0 (x)
1 1 1
0 0 1
Résumé :
- on se ramène à la résolution d'un système d'équation du type Ax=B
- A=p(T) où T est une matrice tridiagonale et p un polynôme facilement calculable
- on utilise l'algorithme de Wiedemann pour résoudre p(T)x=B
- on calcule vp\^k(T)u en calculant les p\^k modulo un polynôme annulateur de T en effectuant les multiplications avec
des FFT adaptées à Z/2Z[X]
- utilisation d'un algorithme pas de bébé/pas de géant pour diminuer le nombre de FFT.
Prologin
2024
