2. Algorithmes et complexité

Cette séance présente quelques techniques ou algorithmes fondamentaux, dont il est bon de connaître l’existence, la complexité ainsi que quelques applications.

Note

Complexité des opérations élémentaires

  • \(A(n)\) : addition dans \(\Z_{\leq 2^n}\) ou \(\K_{\leq n}[x]\)
    • \(O(n)\)
  • \(M(n)\) : multiplication dans \(\Z_{\leq 2^n}\) ou \(\K_{\leq n}[x]\)
    • algorithme naïf \(O(n^2)\)
    • décomposition de Karatsuba \(O(n^{1.58})\)
    • algorithme de Shönhage-Strassen (FFT) \(O(n\log(n)\log\log(n))\)
  • multiplication dans \(M_n(\K)\)
    • algorithme naïf \(O(n^3)\)
    • décomposition de Strassen \(O(n^{2.81})\)

2.1. Exponentiation rapide

Note

Le calcul naïf de \(a^n\) fait intervenir \(n\) multiplications. En écrivant

\[a^n = \bigl( a^{\floor{\frac n2}} \bigr)^2\times a^{n\mod 2}\]

on remarque que l’on peut économiser environ la moitié des multiplications pour le prix d’une mise au carré. On peut appliquer cette idée récursivement.

Pour développer cette idée, décomposons \(n\) en base deux

\[n=e_0+e_1 2+e_2 4+\dots e_k 2^k\]

on peut écrire

\[\begin{split}a^n &= a^{\sum e_i 2^i}\\ &= a^{e_0} \times (a^2)^{e_1} \times \dots (a^{2^k})^{e_k}\end{split}\]

ou bien

\[a^n = a^{e_0} \bigl( a^{e_1} \bigl(\dots a^{e_k}\bigr)^2\dots\bigr)^2\bigr)^2.\]

Dans un cas comme dans l’autre, l’expression ne fait plus intervenir que \(k\) multiplications au plus et \(k\) carrés, avec \(k=\ceil{\log_2(n)}\).

  1. Écrire la fonction pow_naif(a,n) qui calcule \(a^n\) de manière récursive.

    solution xcas ⇲ ⇱ solution xcas
    pow_naif(a,n):={
      if(n==0) return(1);
      return( a*pow_naif(n-1) );
      }
    
  2. Écrire une fonction pow_rec(a,n) qui calcule \(a^n\) par exponentiation rapide de manière récursive. La tester en calculant pow_rec(a,19).

    solution xcas ⇲ ⇱ solution xcas

    Il y a deux solutions à ce problème, selon l’endroit où l’on prend le carré:

    pow_rec(a,n):= {
      if(n==0) return 1;
      if(irem(n,2)) return a*pow_rec(a,iquo(n,2))^2;
      return pow_rec(a,n/2)^2;
    }
    

    ou bien:

    pow_rec(a,n):= {
      if(n==0) return 1;
      if(irem(n,2)) return a*pow_rec(a^2,iquo(n,2));
      return pow_rec(a^2,n/2);
    }
    
    solution sage ⇲ ⇱ solution sage

    Il y a deux solutions à ce problème, selon l’endroit où l’on prend le carré:

    def pow_rec(a,n):
      if n==0: return 1
      r = pow_rec(a^2,n//2)
      if n%2:
        return r*a
      else:
        return r
    

    ou bien:

    def pow_rec(a,n):
      if n==0: return 1
      r = pow_rec(a,n//2)^2
      if n%2:
        return r*a
      else:
        return r
    
  3. Écrire la fonction pow_iter(a,n) qui effectue le même calcul de manière itérative.

    solution xcas ⇲ ⇱ solution xcas
    pow_iter(a,n):= {
      local(r:=1);
      while(n>0) {
        if (n%2==1) {r := r*a};
        a := a^2;
        n := iquo(n,2);
      };
      return(r);
    }
    
    solution sage ⇲ ⇱ solution sage

    Si l’on ne veut pas calculer a priori la décomposition binaire de \(n\), la solution la plus simple correspond au premier algorithme récursif.

    def pow_iter(a,n):
      r = 1
      while n>0:
        if n%2:
          r *= a
        a = a^2
        n = n//2
      return r
    
  4. Exemple à retenir : quel est le coût du calcul des trois derniers chiffres de \(3^{123456789}\) ?

    solution xcas ⇲ ⇱ solution xcas

    Pour obtenir les trois derniers chiffres, il faut prendre le reste modulo 1000. Et puisque la projection \(\Z\to \Z/1000\Z\) est un morphisme d’anneaux, il faut surtout effectuer tout le calcul modulo 1000 en écrivant:

    pow_rec(3%1000,123456789)
    

    et non pas prendre le reste à la fin, ce qui demande de calculer un nombre gigantesque (ce que Xcas ne peut pas faire):

    pow_rec(3,123456789) % 1000
    

    En faisant bien les choses, la complexité du calcul est réduite à une exponentiation rapide de

2.1.1. Application : suite de Fibonacci

On considère la suite classique \(f_0=f_1=1\), \(f_{n+1}=f_n+f_{n-1}\).

  1. Programmer le calcul récursif naïf de la suite. Quelle est sa complexité ? Calculer \(f_{34}\).

    solution sage ⇲ ⇱ solution sage
    def fib_naif(n):
      if n<2: return 1
      return fib_naif(n-1)+fib_naif(n-2)
    

    La complexité du calcul obéit à la même récurrence, si bien que le calcul de \(f_n\) est de complexité \(O(f_n)\), donc exponentielle en \(n\).

  2. Programmer le calcul des \(f_n\) à l’aide d’une fenêtre glissante. Quelle est la complexité ? Calculer \(f_{10^5}\).

    solution sage ⇲ ⇱ solution sage

    Il s’agit de conserver à chaque étape les deux dernières valeurs pour ne pas les recalculer.

    def fib_glis(n):
      ff = (1,1)
      for i in xrange(n-1):
        ff = (sum(ff),ff[0])
      return ff[0]
    
    fib_glis(10^5)
    

    À présent la complexité fait intervenir \(n\) additions de taille croissant linéairement vers \(n\), ce qui fait une complexité quadratique \(O(n^2)\).

  3. En écrivant la récurrence linéaire

    \[\begin{split}F_n = A F_{n-1}, \text{ où } A = \left(\begin{array}{cc}1 & 1 \\ 1 & 0\end{array}\right) \text{ et } F_n = \left(\begin{array}{c}f_{n+1}\\f_{n}\end{array}\right)\end{split}\]

    et l’expression de \(F_n=A^nF_0\), calculer à nouveau \(f_{10^5}\). Quelle est la complexité du calcul, pourquoi ?

    solution sage ⇲ ⇱ solution sage
    def fib_mat(n):
      A = matrix(2,2,[1,1,1,0])
      F0 = vector([1,1])
      return (A^n*F0)[0]
    
    fib_mat(10^4)
    
    solution xcas ⇲ ⇱ solution xcas
    fib_mat(n):={
      local(A);
      A := matrix(2,2,[1,1,1,0]);
      return( (A^n)[0,0] );
    };
    

    La complexité est à présent de \(O(\log(n))\) multiplications matricielles \(2\times 2\), d’entiers dont la taille croît en \(O(n)\), d’où une complexité dominée par les derniers produits en \(O(n\log(n)\log\log n)\)

  4. Comparer les temps de calcul de \(f_{2^{25}-1}\) et \(f_{2^{25}}\), commenter.

  5. Rappeler comment on obtient l’expression algébrique exacte

    \[f_n = \frac{\phi^{n+1}-(1-\phi)^{n+1}}{\sqrt{5}}, \phi = \frac{1+\sqrt{5}}2\]

    et utilisez-la pour calculer \(f_n\).

    solution xcas ⇲ ⇱ solution xcas
    phi := (1+sqrt(5))/2;
    fib_alg(n) := simplify( (phi^(n+1)-(1-phi)^(n+1))/sqrt(5) );
    

    Pour le logiciel, cela revient à effectuer les mêmes calculs qu’avec la matrice.

  6. Quel est le nombre de chiffres de \(f_{f_{f_{f_{f_4}}}}\) ?

2.2. Descente binaire équilibrée (binary splitting)

Cette technique reprend aussi le principe diviser pour régner : il s’agit de résoudre un problème en le décomposant en deux sous-problèmes de complexité inférieure et comparable.

2.2.1. Exemple de la fonction factorielle

  1. Écrire une fonction fact(n) qui calcule la factorielle de \(n\). Quelle est sa complexité binaire ?

    solution sage ⇲ ⇱ solution sage
    def fact(n):
      if n<=1: return 1
      return n*fact(n-1)
    
    def fact(n):
      return prod( k for k in xrange(1,n) )
    
    solution xcas ⇲ ⇱ solution xcas

    Avec une solution récursive

    fact(n):= {
      if(n==0) return(1);
      return(n*fact(n-1));
    }
    

    ou bien de manière itérative:

    fact(n):= {
      local(r:=1);
      for(k:=2;k<=n;k++) { r:= r*k; }
      return(r);
      }
    
  2. Implanter la variante de binary splitting, qui consiste à faire uniquement des multiplication entre opérandes de même taille en écrivant

    \[n! = \bigl(1\times 2\times\cdots \floor{\frac{n}2}\bigr) \bigl(\floor{\frac{n}2+1}\times \cdots (n-1) \times n\bigr)\]

    Quelle est la complexité ? Comparer les temps d’exécution pour \(200000!\).

    solution sage ⇲ ⇱ solution sage
    def fact_bin(n2,n1=1):
      """ computes recursively n2!/n1! """
      if n1 == n2: return n1
      mid = (n2+n1)>>1
      if mid == n1: return n1*n2
      return fact_bin(n2,mid)*fact_bin(mid-1,n1)
    
    solution xcas ⇲ ⇱ solution xcas

    On écrit une fonction qui calcule récursivement \(\frac{n_2!}{n_1!}\)

    fact_split(n2,n1):= {
      if(n2==n1) return(1);
      if(n2==n1+1) return(n2);
      local(mid:=iquo(n1+n2,2));
      return( fact_split(n2,mid)*fact_split(mid,n1) );
      }
    

2.3. Multiplication de Karatsuba

  1. Vérifier l’identité de Karatsuba

    \[\begin{split}(aX+b)(cX+d) = (ac)X^2\\ + [(a+b)(c+d)-(ac)-(bd)]X+(bd)\end{split}\]

    et compter le nombre de multiplications et d’additions qu’elle fait intervenir par rapport au développement usuel du membre de gauche.

    solution ⇲ ⇱ solution

    Au lieu de 4 multiplications et 1 addition dans le produit classique, cette identité fait intervenir 3 multiplications et 4 additions.

  2. Soient \(P\) et \(Q\) deux polynômes de \(\K[x]\) de degré \(d=2^n\). La multiplication de Karatsuba consiste à écrire l’identité ci-dessus en posant d’abord \(X=x^{d/2}\), puis à réappliquer récursivement la même méthode pour chacun des produits de taille \(d/2=2^{n-1}\).

    Écrire une fonction karatsuba(n) qui compte le nombre de multiplications de \(\K\) nécessaires pour effectuer le produit de \(P\) et \(Q\).

    solution xcas ⇲ ⇱ solution xcas
    karatsuba(n):= {
      if(n==0) return(3);
      return(3*karatsuba(n-1));
    };
    
  3. Quelle est la complexité de la multiplication de Karatsuba ?

    solution ⇲ ⇱ solution

    Pour des entiers de \(n\) bits, on effectue une descente récursive de longueur \(\log_2{n}\), soit \(3^{\log_2 n}\) multiplications élémentaires (et de l’ordre de \(7^{\log_2n}\) additions que l’on négligera), donc la complexité est de \(O(n^{\log_2(3)})=O(n^{1.59})\).

2.4. Schéma de Hörner

Note

Le schéma d’évaluation de Hörner consiste à écrire

\[\begin{split}P(x) &= \sum_{k=0}^n a_k x^k\\ &= \bigl(\dots (a_n x + a_{n-1})x+\dots +a_1\bigr)x+a_0\end{split}\]

ce qui permet de calculer la valeur \(P(x)\) à l’aide de \(n\) multiplications et additions.

Remarquer la similitude avec l’exponentiation rapide.

  1. Écrire une fonction horner(P,a) qui évalue le polynôme \(P\) en \(a\). La tester.

2.5. Méthode de Newton

  1. Pour calculer numériquement la racine carrée de \(a\), on définit la suite

    \[x_0 = 1, x_{n+1} = \frac{x_n+a/x_n}2\]
  2. Écrire une procédure Newton(a) qui effectue dix itérations de cette suite, et affiche à chaque étape \(n\) les décimales communes à \(x_n\) et \(x_{n-1}\). On fera le calcul en précision 200.

    solution xcas ⇲ ⇱ solution xcas

    Le nombre de décimales communes à deux nombres est le log en base 10 de leur différence relative.

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    Digits:=200;
    Newton(a) := {
      local(b:=1);
      local(c:=0);
      local(k);
      local(n);
      for(k:=0;k<5;k++) {
        c := b;
        b := evalf((b+a/b)/2);
        n := floor(1-logb(abs((c-b)/c),10));
        if(n>0) print( left(string(c),n) );
        };
    };
    
  3. Quelle est la vitesse de convergence de la suite ? Combien d’itérations faut-il pour obtenir un million de chiffres de \(\sqrt{5}\) ?

    solution xcas ⇲ ⇱ solution xcas

    La précision double à chaque itération, la suite converge exponentiellement vite en nombre de chiffres corrects.

    Il faut seulement \(\log_2(10^6)\approx20\) itérations pour obtenir un million de décimales de \(\sqrt{5}\).