Nous tenons à remercier tout particulièrement M. Raymond Seroul, qui a suivit l’élaboration de ce mémoire.

Nous remercions également le personnel de la bibliothèque de l’Institut de Mathématiques de Strasbourg.

Egalement, les salles C12, C15 du même Institut de Mathématiques dont les tableaux nous ont souvent éclairés !

Remerciements *

Sommaire *

Notations, Conventions *

Introduction *

I. Méthode de travail *

II. Conventions *

III. Langages utilisés *

IV. Plan *

Première Partie : Algèbre Linéaire *

I. Bases Mathématiques *

A. Introduction *

B. Représentation matricielle *

II. Algorithmes classiques *

A. Introduction *

B. Pivot de Gauss / Résolution de systèmes linéaires *

C. Déterminant *

D. Noyau et image *

Deuxième Partie : Le vif du sujet *

I. Contexte théorique *

A. Matrices à coefficients entiers *

B. Modules *

1. Z-Modules *

2. Z-Modules libres de rang fini *

C. Algèbre bilinéaire *

1. Forme bilinéaire, forme quadratique *

2. Représentation matricielle *

D. Réseaux *

II. Algorithmes élémentaires *

A. Forme Normale d’Hermite *

B. Procédé d’orthogonalisation de Gram-Schmidt *

III. Algorithme LLL *

A. Réduction de réseaux *

1. Algorithme LLL original *

2. Algorithme LLL avec insertions profondes *

3. Algorithme LLL intégral *

B. Algorithme MLLL pour des vecteurs dépendants *

Troisième Partie : Applications *

I. Dépendance linéaire *

II. Recherche de petits vecteurs d’un réseau *

A. Recherche de petits vecteurs de ⁿ *

B. Algorithme de Cholesky *

C. Recherche de petits vecteurs d’un réseau *

Conclusion *

Bibliographie *

Résumé *

I. Mots Clefs *

II. Condensé *

Le travail que nous proposons est une étude détaillée de l’algorithme de réduction des réseaux proposé par A.K. Lenstra, H.W. Lenstra et L. Lovacz en 1982. Cette étude est essentiellement basée sur l’ouvrage de H. Cohen qui rassemble de nombreux algorithmes mathématiques.

1. Méthode de travail

La méthode de travail utilisée consiste en une description théorique des problèmes et objets, puis une mise en pratique au niveau algorithmique et programmation. La partie théorique tentera de pas être trop exhaustive : il n’est pas ici question d’un cours d’algèbre général. Quant à la partie pratique, elle se veut la plus claire et la plus simple possible.

2. Conventions

Les conventions et notations adoptées dans ce mémoire sont énumérées au fil du besoin, outre les signifiants relevés plus haut dans le lexique.

3. Langages utilisés

Le choix d’un langage de programmation des algorithmes est un problème délicat. Nous avons décidé de retenir le langage Pascal de par son universalité dans le milieu mathématique. La programmation restera simple, voire simpliste dans la plupart des cas. Il ne s’agit pas ici d’un mémoire d’informatique, mais de mathématiques.

Ce choix se justifie également en opposition à celui de H. Cohen : celui-ci implémente en effet ses algorithmes en GP (programme de calcul), ce qui le conduit à une description " logique " (à la Knuth), faite de renvois enchevêtrés et inextricables. La lisibilité informatique de ces algorithmes s’avère très délicate.

Pour des raisons de confort de programmation, les premiers tests sont faits avec Mathematica.

4. Plan

Le travail se décompose en trois grands axes. Dans un premier temps, nous donnerons quelques rappels d’algèbre linéaire : les bases théoriques nécessaires à une bonne compréhension du problème et quelques algorithmes classiques.

Par la suite, nous entrerons dans le vif du sujet : la théorie, essentiellement algébrique, et la mise en pratique de quelques algorithmes qui servirons de base par la suite et une étude plus détaillée des algorithmes LLL.

Enfin, nous envisagerons deux applications parmi les plus importantes de cet algorithme.

1. Bases Mathématiques

1. Introduction

Nous aurons ici besoin des notions de corps, d’espace vectoriels et de base. Nous rappelons a ce propos, que tout espace vectoriel admet une base. Nous travaillerons exclusivement avec des espaces vectoriels de dimension finie.

Nous nous intéresserons également à la structure d’espace euclidien, c’est-à-dire, d’espace vectoriel muni d’un produit scalaire. La plupart du temps, nous considérerons des espaces réels muni du produit scalaire usuel, bien que cette théorie s’étende à des espaces plus généraux.

2. Représentation matricielle

Un K espace vectoriel V étant un objet abstrait, on le modélise par la donnée d’une base sous forme de matrice à coefficients dans K. On peut alors également représenter un sous-espace (sous-ensemble stable) W de V par une matrice.

Une fonction linéaire f d’un espace vectoriel V dans un espace vectoriel W est représentée par la donnée d’une matrice A dont les colonnes représentent les coordonnées des images des vecteurs de base de V dans W.

Dans la pratique, nous ne manipulerons donc que des matrices, mais il ne faudra pas perdre de vue l’interprétation des calculs effectués.

On notera M(V) la matrice représentant V dans la base canonique, M_n,m(K) l’ensemble des matrices n ´ m à coefficients dans K et GL(n,K) l’ensemble des matrices à coefficients dans K carrées de dimension n inversibles.

2. Algorithmes classiques

1. Introduction

Les quelques algorithmes qui suivent donnent des outils de base sur les objets introduits plus haut.

2. Pivot de Gauss / Résolution de systèmes linéaires

On cherche à résoudre le problème :

trouver x tel que Ax = b où A Î GL(n,K), x et b sont deux vecteurs colonnes.

On recherche dans la première colonne un " pivot " non nul. On échange les lignes pour faire venir le pivot en première position et on normalise la ligne. Par combinaisons linéaires, on annule les coefficients de la ligne autres que le pivot qui vaut alors 1. On effectue les opérations simultanément sur une matrice identité. Et on passe à la ligne suivante. Finalement, on aura par des opérations élémentaires transformé la matrice A en identité et la matrice identité de départ en A^-1. Ainsi, on a x = A^-1.b .

Ce qui donne un code pascal semblable à ce qui suit. Les fonctions asinitrotantes ne figurent pas ici.

 

Function FirstPivot(Var A : Matrix; dA : Dimension; j : Byte) : Byte;
(***********************************************************************)
(* Cherche le 1^er pivot dans la colonne j, renvoie 0 si il n'y en a pas *)
(***********************************************************************)
Var
k,FP : Byte;
Begin
FP:=0;
k:=j;
While (FP=0) And (k<=dA) Do
If A[k,j]<>0
Then FP:=k
Else k:=k+1;
FirstPivot:=FP;
End;
 
Procedure Gauss(A : Matrix; dA : Dimension;
Var InvA : Matrix; Var Inversible : Boolean);
(***********************************************************************)
(* Inverse la matrice A dans InvA. Inversible est True si A inversible *)
(***********************************************************************)
Const
Introuvable = 0;
Var
j,LignePivot : Byte;
 
Procedure WorkLine(i,j : Byte; Var A,InvA : Matrix; dA : Dimension);
(*****************************)
(* Travaille avec la ligne j *)
(*****************************)
Var
Coef : Real;
 
Procedure DoTheWash(c : Byte; Var A, InvA : Matrix; dA : Dimension);
(************************************************)
(* Elimine les coef de la colonne c sauf A[c,c] *)
(************************************************)
Var
k : Byte;
Coef : Real;
Begin
For k:=1 To dA Do
If k<>c
Then
Begin
Coef:=-A[k,c];
LinearCombination(c,k,A,dA,Coef);
LinearCombination(c,k,InvA,dA,Coef);
End;
End;
 
Begin
If i<>j
Then
Begin
SwapLines(i,j,A,dA);
SwapLines(i,j,InvA,dA);
End;
Coef:=1/A[j,j];
ScalarLine(j,A,dA,Coef);
ScalarLine(j,InvA,dA,Coef);
DoTheWash(j,A,InvA,dA);
End;
 
Begin
UnitMatrix(InvA,dA);
j:=1;
Inversible:=True;
While Inversible And (j<=dA) Do
Begin
LignePivot:=FirstPivot(A,dA,j);
If LignePivot=Introuvable
Then Inversible:=False
Else
Begin
WorkLine(LignePivot,j,A,InvA,dA);
j:=j+1;
End;
End;
End;

3. Déterminant

On peut toujours rendre triangulaire supérieure, une matrice A à coefficients dans un corps, par la méthode Gauss. Le déterminant est alors le produit des éléments diagonaux au signe près (on multiplie par –1 à chaque échange de colonnes).

4. Noyau et image

Etant donnée une matrice A, on définit son noyau par la donnée des vecteurs x tels que Ax=0.

L’image est le supplémentaire du noyau dans l’espace (en dimension finie, toujours…).

1. Contexte théorique

1. Matrices à coefficients entiers

On note GL(n, Z) l’ensemble des matrices inversibles à coefficients entiers

Théorème : GL(n, Z) = {A Î M_n,n(Z), Det(A)= ±1}

Démonstration : On a A.A^-1=1

Comme la fonction Det est un morphisme, on a :

Det(A.A^-1)=Det(1)

Û Det(A).Det(A^-1)=1

Û Det(A)=Det(A^-1)^-1

Or, les seuls éléments inversibles de Z sont 1 et –1. þ

2. Modules

1. Z-Modules

On appelle Z-Module de rang fini un groupe abélien isomorphe à un produit de quotients Z/n_iZ et de Zⁿ.

2. Z-Modules libres de rang fini

On appelle Z-Module libre de rang fini n un groupe abélien isomorphe à Zⁿ.

On représente un Z-Module par la matrice de représentation d’une base dans la base canonique. On voudrait alors définir les mêmes notions que pour les corps : inversion, noyau, etc.

4. Algèbre bilinéaire

1. Forme bilinéaire, forme quadratique

Définition : Une forme bilinéaire b est une application de V dans K, où V est un espace vectoriel de dimension 2, linéaire en les deux variables.

Définition : Soit K un corps de caractéristique différente de 2, et soit V en K-espace vectoriel. On dit qu’un application q de V dans K est une forme quadratique si elle satisfait les conditions suivantes :

" l Î K, " x Î V, q(l x)=l ²q(x)

En posant b(x,y)= (q(x+y)-q(x)-q(y)) alors b est une forme bilinéaire (symétrique).

On a alors q(x)=b(x,x).

Remarque : Si K= , on dit que q est définie positive (strictement) si " x Î V, q(x)>0.

2. Représentation matricielle

On représente une forme bilinéaire par la donnée de la matrice b(e_i,e_j) où (e_i)_iÎ   une base de V.

5. Réseaux

On appellera réseau, un -module libre de rang fini muni d’une forme quadratique définie positive.

Un cas particulier est de prendre un produit scalaire pour forme bilinéaire et donc la norme associée pour forme quadratique. C’est dans ce cadre essentiellement que nous poursuivrons.

3. Algorithmes élémentaires

1. Forme Normale d’Hermite

Définition : on dit qu’une matrice M est en Forme Normal d’Hermite (ou HNF pour Hermite Normal Form) si elle est de la forme :



Théorème : Soit A Î M_n,m(). Il existe une unique matrice H = (b_i,j) Î M_n,m() en HNF telle que H = AU où U Î GL(n,).

Remarque : U n’est pas forcément unique.

2. Procédé d’orthogonalisation de Gram-Schmidt

Théorème (Gram-Schmidt) : Un espace vectoriel euclidien admet une base orthogonale.

Démonstration : On définit par récurrence la nouvelle base où , avec les notations usuelles.

La nouvelle base ainsi définie est trivialement orthogonale. þ

4. Algorithme LLL

1. Réduction de réseaux

Il est évident que l’algorithme de Gram-Schmidt ne peut se généraliser aux -modules à cause de la division : on ne retombe généralement pas sur un point du réseau… Cependant, on peut espérer approcher une base qui soit le plus orthogonale possible, en un sens qu’il reste définir. C’est cette définition et un algorithme pour trouver cette base qui font l’objet central de ce mémoire.

Définition : En gardant les notation plus haut, une base b₁, b₂, …, b_n d’un réseau est dite LLL-réduite si :

½ m _{i,j ½}   " 1 j < i n

et b  + m _i,i-1b ²   b  ² " 1 < i n

Remarque : la deuxième inégalité est équivalente à b ² ( -m  )b  ²

Cette inégalité est appelée condition de Lovasz.

1. Algorithme LLL original

Une version Pascal de l’algorithme mettant en œuvre cette réduction est proposée ici.

On suppose que (b_k) est une base. On procède alors par récurrence : on suppose que les k-1 premiers vecteurs sont déjà LLL-réduits. On peut initialiser cette récurrence pour k=2. On cherche maintenant à réduire le k-ième vecteur, de sorte que ½ m _{i,j ½}   "  1 j < i n. Ceci peut être fait en remplaçants b_k par , a_j Î  : on suppose ½ m _{i,j ½}   "  l j <k (vrai pour l=k). En posant alors q comme l’arrondi de m _k,l, b_k=b_k-q.b_l, m _k,j n’est pas modifié pour j>l (car b  est orthogonal à b_l pour l<j). On remplace alors m _k,l par _{m k,l}-q et ½ m _k,l-q½  . On a alors ½ m _k,j½   pour l-1<j<k.

Maintenant que la réduction a été effectuée, il faut aussi que b_k vérifie la condition de Lovasz. Si (par miracle !) b_k satisfaisait cette condition, on passe au vecteur suivant, s’il existe. Sinon, on échange les vecteurs b_k et b_k-1 et on décrémente k, considérant que seuls les k-2 premiers vecteurs sont LLL-réduits.

Pour le calcul des coefficients de Gram (m _i,j), la solution retenue est un calcul au fur et à mesure que les coefficients sont nécessaires.

Remarques :

ü l’échange des vecteurs b_k et b_k-1 induit des changements sur les coefficients de Gram.

ü On calcule tout au long de l’algorithme la matrice de passage H.

Ceci donne l’algorithme qui suit

Procedure BasisLLL(b_in : Vector; Var b_out : Vector; Var H : Matrix);
k, kmax : Integer;
B : Vector;
Lovasz : Boolean;
Begin
(* Begin of Reduction *)
(* Variables Initialization *)
k:=2;
kmax:=1;
SetVectors(b_out[1], b_int[1]);
B[1]:=Dot(b_in[1], b_in[1]);
IdentityMatrix(H, n);
(* Main Loop *)
Repeat
If k>kmax Then Gram(k, kmax);
(* Lovasz's Loop *)
Repeat
Red(k, k-1);
Lovasz:=B[k]>=(0.75-Mu[k,k-1]^2)*B[k-1];
If Not Lovasz
Then Begin
Swap(k);
k:=max(2, k-1);
End;
Until Lovasz;
(* End Lovasz's Loop *)
For l:=k-2 DownTo 1 Do Red (k,l)
k:=k+1;
Until k>n;
(* End Main Loop *)
(* End of Reduction *)
End;

Procedure Red(k, l : Integer);
Var i, q : Integer;
Begin
If Mu[k,l]>0.5 Then
Begin
q:=Round(Mu[k,l]);
"b_in[k]:=b_in[k]-q*b_in[l]";
"Hk:=Hk-q*Hl"; (* Hk est la K-ieme colonne de H *)
Mu[k,l]:=Mu[k,l]-q;
For i:=1 To n-1 Do
Mu[k,i]:=Mu[k,i]-q*Mu[l,i];
End;
End;

Procedure Swap(k : Integer);
Var i,j : Integer;
BT, Temp, Mu : Real;
Begin
ExchangeVectors(b_in[k],b_in[k-1]);
ExchangeMatrixCol(H, k, k-1);
If k>2 Then
For j:=1 To k-2 Do
Begin
Temp:=Mu[k,j];
Mu[k,j]:=Mu[k-1,j];
Mu[k-1,j]:=Temp;
End;
Mu:=Mu[k,k-1];
BT:=B[k]+Mu^2*B[k-1];
Mu[k,k-1]:=Mu*B[k-1]/BT;
B[k]:=B[k-1]*B[k]/BT;
B[k-1]:=BT;
For i:=k+1 To kmax Do
Begin
Temp:=Mu[i,k];
Mu[i,k]:=Mu[i,k-1]-Mu*Temp;
Mu[i,k-1]:=Temp+Mu[k,k-1]*Mu[i,k];
End;
End;


Détails :

La procédure Red : Si m _k,l 0,5 , b_k est déjà de norme minimale. Sinon, on va soustraire le multiple maximum de b_l pour obtenir un m _k,l 0,5. On fait les mêmes opérations sur H et il faut rectifier les m _k,i pour i<l.

La procédure Swap : En changeant les vecteurs b_k et b_k-1 (les colonnes k et k-1 de H), il faut modifier les relations de Gram-Schmidt. Pour j<k-2, il suffit d’échanger m _k,j avec m _k-1,j. Pour m _k,k-1, …m _k,kmax, il faut inverser le système :

pour k-1 l kmax. (ce qui induit des changements sur B_k).

C’est-à-dire, en échangeant b_k et b_k-1 dans ce système.

Ce qui donne le code optimisé plus haut.

Remarques :

ü On peut modifier l’algorithme pour récupérer uniquement la base ou la matrice H.

ü Cet algorithme faisant intervenir des quantités réelles, dans la pratique, il se peut que le résultat ne soit pas LLL-réduit, mais en général, ce sera le cas. On pourra se référer à l’algorithme LLL intégral plus loin pour contourner cette difficulté.

Preuve de l’algorithme : On suppose pour la preuve les b  initialisés aux b_j. Il est alors clair que (b ) est une base du réseau, tout au long de l’algorithme.

Soit et . Il est évident que D est à valeurs entières.

L’unique partie de l’algorithme ou les B_j sont modifiés (et donc D) est le sous-algorithme SWAP dans lequel on modifie un seul B_k que l’on diminue au moins d’un facteur   par la condition de Lovasz. Donc D est diminuée d’un facteur   également et reste bien entendu à valeurs entières (les B_j sont toujours des normes au carré de vecteurs à coordonnées entières). Ceci montre que la boucle centrale mettant en jeu l’algorithme SWAP est finie. Comme c’est le seul endroit où k est decrémenté, on a montré que l’algorithme se termine. Il est évident que la base obtenue est LLL-réduite.

þ

2. Algorithme LLL avec insertions profondes

A la place d’échanger simplement b_k et b_k-1, on peut directement échanger b_k avec b_i, i<k tel que : . Ceci donne alors l’algorithme schématique qui suit, dont nous ne donnerons pas d’équivalent Pascal (il ne nous sera pas utile par la suite).



3. Algorithme LLL intégral

Le but de cet algorithme est de proposer une version de l’algotihme LLL original ne faisant intervenir que des calculs sur entiers, ce qui garantit la stabilité de l’algorithme.

Théorème et Corollaire :

On suppose que la matrice de Gram (<b_i,b_j>) est à coefficients entiers (c’est le cas ici).

Soit d_i=det(<b_r, b_s> _{1 r, s i}) = B₁…B_i

Alors pour tout i fixé et pour tout j < i, on a :

d_i-1B_i Î ,

d_{im i,j} Î ,

pour tout m tel que j < m i, Î .

Soit l _i,j = d_{j m i,j} pour j < i et l _i,i=d_i. Alors pour j i fixé, on a :



Alors u_k Î et u_j-1=l _i,j.

Preuve : cf. Ouvrage H. Cohen p. 91-92.

Avec ces notations, on aboutit à l’algorithme suivant :

Procedure IntegralLLL(b_in : Vector; Var b_out : Vector; Var H : Matrix);
Var k, kmax, l : Integer;
Lovasz : Boolean;

Begin
(* Begin Of Reduction *)
(* Gram-Schmidt(b_in, b_out, H); *)
(* ou H=IdentityMatrix(n); b_our:=b_in; *)
(* Variables Initialization *)
k:=2; kmax:=1;
B[0]:=1;
d[1]:=Dot(b_out[1],b_out[1]);
(* Main Loop *)
Repeat
If k>kmax Then
IntegralGram(k, kmax);
(* Lovasz's Loop *)
Repeat
RedI(k, k-1);
Lovasz:=d[k]*d[k-2]>=(0.75*d[k]^2-Lambda[k,k-1]^2);
If Not Lovasz Then
Begin
SwapI(k);
k:=Max(2,k-1);
End;
Until Lovasz;
(* End Of Lovasz's Loop *)
For l:=k-2 DownTo> 1 Do
RedI(k,l);
k:=k+1;
Until k>n;
End;

Procedure RedI(k, l : Integer);
Var i,q : Integer;
Begin
If Abs(2*Lambda[k,l]>d[l]) Then
Begin
q:=Round(Lambda[k,l]/d[l]);
(* Hk:=Hk-q*Hl; *)
Lambda[k,l]:=Lambda[k,l]-q*d[l];
For i:=1 To l-1 Do
Lambda[k,i]:=Lambda[k,i]-q*Lambda[l,i];
End;
End;

Procedure SwapI(k : Integer);
Var i, j, Temp, LambdaT, BT : Integer;
Begin
ExchangeMatrixCol(H, k, k-1);
If k>2 Then
For j:=1 To k-2 Do
Begin
Temp:=Lambda[k,j];
Lambda[k,j]:=Lambda[k-1,j];
Lambda[k-1,j]:=Temp;
End;
LambdaT:=Lambda[k,k-1];
BT:=(d[k-2]*d[k]+LambdaT^2)/d[k-1];
For i:=k+1 To kmax Then
Begin
Temp:=Lambda[i,k];
Lambda[i,k]:=(d[k]*Lambda[i,k-1]-LambdaT*Temp)/d[k-1];
Lambda[i,k-1]:=Temp+LamdaT*Lambda[i,k]/BT;
End;
d[k-1]:=BT;
End;

Procedure IntegralGram(k : Integer; Var kmax : Integer);
Var i,j : Integer;
u : Integer;
Begin
kmax:=k;
For j:=1 To k Do
Begin
u:=Dot(b_in[k],b_in[j]);
For i:=1 To j-1 Do
u:=(d[i]*u-Lambda[k,i]*Lambda[j,i])/d[i-1];
If j<k Then Lambda[k,j]:=u;
If j=k Then d[k]:=u;
End;
If d[k]=0 Then
Begin
Writeln('Les b_in ne forment pas une base.');
Halt(0);
End;
End;

Cet algorithme n’appelle pas de commentaire particulier. Sa preuve est identique à celle du précédent.

2. Algorithme MLLL pour des vecteurs dépendants

1. Dépendance linéaire

On cherche à découvrir une dépendance -linéaire entre les n nombres complexes (z₁,… , z_n), si il en existe une.

Supposons tout d’abord que ces nombres soient réels. Pour un entier N assez grand, on considère la forme quadratique en les a_i suivante :



Cette forme est définie positive car elle est toujours positive (somme de carrés) et est nulles si et seulement si tous les carrés sont nuls, i.e. et donc a₁ = 0 et a est nul. Q définit donc un produit scalaire euclidien sur  ⁿ si z₁ ¹ 0.

Soit a l’ordre des a_i. On peut choisir a ⁿ < N < a ²ⁿ, par exemple : N = [a ^1,5n].

En petit vecteur pour Q vérifie alors : est petit ainsi que les a_i pour i > 1. Si les z_i sont effectivement -linéaire dépendants, alors on aura trouvé une relation si est nul.

On va en fait appliquer l’algorithme LLL a la base canonique (e_i) de ⁿ pour la forme quadratique Q. On a alors le lemme suivant :

Lemme : Avec les notations introduites plus haut, si on exécute l’algorithme d’orthogonalisation de Gram-Schmidt (pour Q) sur la base (e_i), on a :

m _i,1 = z_i/z₁, pour 2 i n

m _i,j = 0, pour 2 j < i n

e  = e_i - (z_i/z₁) e₁

B₁ = Nz  et B_k = 1 pour 2 k n.

Démonstration : soit b la forme bilinéaire définie par Q. On rappelle que

b(x,y) =  (Q(x+y)-Q(x)-Q(y)) pour deux vecteurs x et y.

soit k 2, on a :



Or : Q(e₁) = N z , Q(e_k) = 1 + N z , Q(e₁+e_k) = 1 + N(z₁+z_k)²

D’où : m _k,1 = z_k/z₁.

Soit 2 j < i n :



Q(e_i) = 1 + N z , Q(e ) = 1, Q(e_i+e ) = 2 + N z 

D’où m _i,j = 0.

B₁ = Q(e₁) = N z  et B_j = Q(e ) = 1. þ

 

On peut facilement modifier ce procédé pour travailler avec des nombres complexes. Dans ce cas on prend la forme quadratique:



Q est alors définie positive et définit un produit scalaire euclidien sur R²ⁿ ( ⁿ) si et seulement si z₁/z₂ Ï  . Si jamais z₁/z₂ est réel, on cherche un couple z_i/z_j qui n’est pas réel et on réordonne. Si il n’en existe pas, on utilise le processus pour les réels 1, z₂/z₁, …, z_n/z₁.

On utilise l’algorithme LLL pour trouver la relation de dépendance entre les z_i la plus " petite " possible, c’est-à-dire les coefficients a_i les plus petits possibles.

Ceci donne l’algorithme suivant :

Procedure LinDep(z : TabComplex; BigN : Integer);
Var b_in, b_out, B : Vector;
k, kmax, i, j : Integer;
H, Mu : Matrix;
Begin
IdentityMatrix(b_in);
For i:=3 To n Do
For j:=3 To i-1 Do
Mu[i,j]:=0;
B[1]:=Re(z[1])*Re(z[1])+Im(z[1])*Im(z[1]);
B[2]:=Im(cprod(z[1],conj(z[2])));
For i:=3 To n Do
B[i]:=1;
For i:=2 To n Do
Mu[i,1]:=Re(cprod(z[1],conj(z[i])))/B[1];
If B[2]<>0 Then
Begin
For i:=3 To n Do
Mu[i,2]:=Im(cprod(z[1],conj(z[i])))/B[2];
B[2]:=BigN*B[2]*B[2]/B[1];
End
Else Begin
For i:=3 To n Do
Mu[i,2]:=0;
B[2]:=1;
End;
(* LLL *)
SetVectors(b_out[1], b_int[1]);
B[1]:=BigN*B[1];
k:=2;
kmax:=n;
IdentityMatrix(H, n);
Repeat
Repeat
Red(k, k-1);
Lovasz:=B[k]>=(0.75-Mu[k,k-1]^2)*B[k-1];
If Not Lovasz
Then Begin
Swap(k);
k:=max(2, k-1);
End;
Until Lovasz;
For l:=k-2 DownTo 1 Do Red (k,l)
k:=k+1;
Until k>n;
(* End of LLL *)
Output(b_out, z);
End;

Remarque :

ü On n’a pas besoin ici de Gram, puisqu’on l’effectue d’un seul tenant au début de la procédure.

2. Recherche de petits vecteurs d’un réseau

On veut trouver les petits vecteurs d’un réseau, c’est-à-dire les vecteurs x tels que Q(x) C, où C est une constante positive entière et Q est la forme quadratique du réseau. On va résoudre ce problème en trois temps : grâce à l’algorithme LLL, on va se placer dans une base plus " sympathique ", chercher dans celle-ci les petits vecteurs et redonner les solutions dans la base originale. On a donc besoin de deux outils : l’un qui sache chercher les petits vecteurs de ⁿ et l’autre qui sache transformer la matrice associée à la forme quadratique sous une forme adéquate : c’est l’algorithme classique de Cholesky.

1. Recherche de petits vecteurs de ⁿ

Le but est de trouver tous les vecteurs x non nuls de ⁿ vérifiant Q(x) C, où C est une constante positive entière et Q est une forme quadratique du type :



On va choisir (arbitrairement) une coordonnée du vecteur et lui affecter la plus grande valeur possible. On va ensuite la décrémenter et incrémenter successivement les coordonnées suivantes. C’est cette idée intuitive que l’on va formaliser…

On va choisir et ensuite , et ainsi de suite jusqu’à obtenir le vecteur. On décrémente ensuite x_n et on recommence. Jusqu’à avoir le vecteur nul…

Ceci donne l’algorithme pascal suivant (complet) :

Program ShortVectors;
Uses WinCrt;
Const n=5;
Type Vecteur = Array[1..n] of Integer;
Var
C,i,j : Integer;
q : Array[1..n,1..n] Of Integer;
Function IsNul(x : Vecteur) : boolean;
Var i : Integer;
IsN : Boolean;
Begin
IsN:=True;
For i:=1 To n Do
If x[i]<>0 Then IsN:=False;
IsNul:=IsN;
End;

Procedure SV(C : Integer);
Var i,j : Integer;
count : Integer;
Z : Real;
T,L,x,U : Vecteur;
NewBounds : Boolean;
Begin
Count:=0;
i:=n;
T[i]:=C;
U[i]:=0;
NewBounds:=True;
Repeat
If NewBounds Then
Begin
Z:=Sqrt(T[i]/q[i,i]);
L[i]:=Trunc(Z-U[i]);
x[i]:=Trunc(-Z-U[i])-1;
End;
Repeat x[i]:=x[i]+1;
i:=i+1;
Until x[i-1]<=L[i-1];
i:=i-1;
If i>1 Then
Begin
T[i-1]:=T[i]-q[i,i]*(x[i]+U[i])*(x[i]+U[i]);
i:=i-1;
U[i]:=0;
For j:=i+1 To n Do U[i]:=U[i]+q[i,j]*x[j];
NewBounds:=True;
End
Else
If Not IsNul(x) Then
Begin
Count:=Count+2;
Write('x = (');
For j:=1 To n-1 Do Write(x[j],', ');
Writeln(x[n],')');
Write('-x = (');
For j:=1 To n-1 Do Write(-x[j],', ');
Writeln(-x[n],')');
Writeln('Q(x) = ',C-T[1]+q[1,1]*(x[1]+U[1])*(x[1]+U[1]));
Writeln;
NewBounds:=False;
End;
Until IsNul(x);
Writeln('Nb de vecteurs trouves : ',Count);
End;


Begin
For i:=1 To n Do
For j:=1 To n Do
If i=j Then q[i,j]:=1 Else q[i,j]:=0;
C:=7;
InitWinCrt;
SV(C);
Writeln('Appuyez sur Return...');
Readln;
DoneWinCrt;
End.

 

La preuve de cet algorithme est triviale : si celui-ci ce termine, par les calculs effectués, on est assuré que les vecteurs x sont solutions du problème. De plus, l’algorithme se termine car la norme (forme quadratique) des vecteurs x diminue dans la boucle principale en restant entière.

2. Algorithme de Cholesky

Le problème est le suivant : soit A une matrice réelle, symétrique définissant la forme quadratique Q. On va chercher une matrice R triangulaire supérieure telle que A= ^tR.R. En fait, on va mettre Q sous la forme



L’algorithme est classique pour trouver R, on pourra se référer pour de plus amples détails (par exemple) au cours d’Analyse Numérique de Licence de M. Komornik (Strasbourg) p. 55 ou de DEUG de M. Akesby (Mulhouse) p. 34. La version qui calcule les coefficients de Q est moins usuelle : on pourra se référer à l’article de U. Fincke et M. Pohst paru dans Mathematics of Computation.

Procedure Cholesky(A : Matrix; Var Q, R : Matrix);
Var i,j,k,l : Integer;
Begin
For j:=1 To n Do
For i:=1 To j Do
Q[i,j]:=A[i,j];
i:=0;
Repeat
i:=i+1;
If i<n Then
Begin
For j:=i+1 To n Do
Begin
Q[j,i]:=Q[i,j];
Q[i,j]:=Q[i,j]/Q[i,i];
End;
For l:=i+1 To n Do
For k:=i+1 To l Do
Q[k,l]:=Q[k,l]-Q[k,i]*Q[i,l];
End;
Until i=n;
For i:=1 To n Do
Begin
R[i,i]:=Sqrt(Q[i,i]);
For j:=1 To i-1 Do R[i,j]:=0;
For j:=i+1 To n Do R[i,j]:=R[i,i]*Q[i,j];
End;
End;


 

3. Recherche de petits vecteurs d’un réseau