Algorithmique de base

Entiers

Note

Théorème :

Soit \(b\) un entier supérieur ou égal à deux, pour tout entier \(n\geq0\) il existe une unique suite d’entiers \(a_k,\dots a_0\) tels que \(0\leq a_i<b\), \(a_k\neq 0\) et

\[n = \sum_{i=0}^k a_i b^i\]

De plus, on a \(k = \floor{ \log_b n }\).

Evaluation :

  1. Ecrire une fonction evalue_base(a,b) qui renvoie l’entier \(n\) dont a est la suite de chiffres en base b.

    Par exemple, evalue_base([1,1,0,1],2) retourne 13 (on veut que la liste donnée soit \([a_k,\dots a_0]\)).

    solution xcas ⇲ ⇱ solution xcas
    evalue_base(a,b) := {
      local(j);
      sum(a[j]*b^j,j,0,length(a)-1);
    }
    

Ecriture en base b

  1. Ecrire une fonction ecriture(n,b) qui renvoie la liste \([a_k,\dots a_0 ]\).

    solution xcas ⇲ ⇱ solution xcas

    \(a_0\) est le reste de la division euclidienne de \(n\) par \(b\), et les chiffres \(a_k,\dots a_1\) forment l’écriture en base \(b\) du quotient.

    On a donc sous forme récursive

    ecriture(n,b) := {
      local q,r;
      if(n=0) return [];
      q := iquo(n,b);
      r := irem(n,b);
      return concat(ecriture(q,b),[r]);
    }
    

    ou sous forme itérative, dans laquelle on initialise dès le début une liste de la bonne taille

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    ecriture(n,b) := {
      local a,k,j;
      k := floor(log(n)/log(b));
      a := newList(k+1);
      for(j:=0;j<=k;j++) {
        a[k-j] := irem(n,b);
        n := iquo(n,b);
      }
      return a;
      }
    

Méthode d’évaluation de Hörner

Note

Soit \(P(x)\in\K[x]\) un polynôme de degré \(r\). Le schéma d’évaluation de Hörner consiste à écrire

\[\begin{split}P(x) &= \sum_{i=0}^r a_i x^i\\ &= \bigl(\dots (a_r x + a_{r-1})x+\dots +a_1\bigr)x+a_0\end{split}\]
  • La première expression fait intervenir \(r\) additions, \(r\) multiplications, et \(r\) exponentiations (que l’on peut réduire à \(r\) multiplications si l’on calcule les \(x^i\) de proche en proche)
  • La seconde fait intervenir \(r\) additions et \(r\) multiplications : il s’agit du schéma de Hörner.
  1. Ecrire une fonction evalue_horner(a,b) qui calcule \(n\) en utilisant le schéma de Hörner.

    solution xcas ⇲ ⇱ solution xcas
    evalue_horner(a,b) := {
      local(j);local(s);
      s := 0;
      for(j:=0;j<length(a);j++) {
        s := b*s + a[j];
      }
      return s;
      }
    
  2. Déterminer la plus petite base \(b\) telle que le nombre dont l’écriture en base \(b\) est formée de \(157\) chiffres 1 soit un nombre premier (utiliser la fonction is_pseudoprime).

    solution xcas ⇲ ⇱ solution xcas
    v := seq(1, k, 1, 157);
    for(b:=2;b<100;b++) { if(is_pseudoprime(evalue_horner(v,b))) break; }
    b;
    

    On peut aller plus vite : pour toute base \(b\) ce nombre vaut \(\sum_{j=0}^{156}b^j=\frac{b^{157}-1}{b-1}\).

    b:=2; while(is_pseudoprime((b^157-1)/(b-1))=0) {b++}; b;
    

    Au passage, il y a une réponse pour \(l=157\) qui est un nombre premier, si \(l\) n’est pas premier ce n’est jamais le cas à cause de la factorisation cyclotomique suivante

    \[b^l-1=\prod_{d\mid l}\Phi_d(b)\]

Exponentiation rapide

  1. Ecrire une fonction puissance_naive(a,n) qui effectue le calcul via la formule élémentaire

    \[a^n = a\times a\times\dots \times a\]

    Quel est le nombre de multiplications effectuées ?

    solution xcas ⇲ ⇱ solution xcas
    puissance_naive(a,n) := {
      if (n==0) return 1;
      return a*puissance_naive(a,n-1);
      }
    

    Cette fonction effectue \(n-1\) multiplications (en négligeant la dernière multiplication par 1).

On peut faire beaucoup mieux.

Note

Théorème :
Soit \((G,\times)\) un semi-groupe, et \(a\in G\). Alors pour tout entier \(n\geq 1\), on peut calculer \(a^n\) avec au plus \(2r\) multiplications \(\times\) de \(G\), où \(r=\floor{\log_2(n)+1}\).

Pour obtenir ce résultat, on écrit \(n\) en base \(2\)

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

on peut alors écrire \(a^n\) sous deux formes

\[a^n = a^{e_0} \times (a^2)^{e_1} \times \dots (a^{2^r})^{e_r}\]

ou bien

\[a^n = a^{e_0} \times \bigl( a^{e_1} \times \bigl(\dots a^{e_r}\bigr)^2\dots\bigr)^2\bigr)^2.\]
  1. On a écrit la fonction récursive suivante

    puissance_rec(a,n):= {
      local r;
      if(n==0) return 1;
      r := puissance_rec(a^2,iquo(n,2));
      if(irem(n,2))
        { return a*r; }
      else
        { return r; }
      }
    

    À quelle écriture correspond-elle ?

    solution xcas ⇲ ⇱ solution xcas

    Cette fonction effectue le calcul selon la première écriture.

  2. Modifier la fonction puissance_rec pour qu’elle corresponde à l’autre écriture.

    solution xcas ⇲ ⇱ solution xcas

    Dans l’autre sens on «met le carré en dehors»

    puissance_rec(a,n):= {
      local r;
      if(n==0) return 1;
      r := puissance_rec(a,iquo(n,2))^2;
      if(irem(n,2))
        { return a*r; }
      else
        { return r; }
      }
    
  3. Ecrire une fonction puissance_iter(a,n) qui calcule \(a^n\) de manière itérative. On essaiera de ne pas calculer d’abord la décomposition de \(n\) en base \(2\), mais de le faire au cours du calcul (comme dans les fonctions récursives).

    solution xcas ⇲ ⇱ solution xcas

    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.

    puissance_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);
    }
    
  4. Modifier la fonction puissance_iter(a,n) pour qu’elle accepte un paramètre b et fasse le calcul en décomposant selon une base \(b\) (on pourra cette fois commencer par calculer la décomposition).

    Tester votre fonction sur le calcul de Mod(3, 2^1024+1)^(2^(2^15)-1). Pour quelles bases obtient-on les meilleurs résultats ?

    solution xcas ⇲ ⇱ solution xcas

    Cette fois on utilise le schéma de Hörner en précalculant l’écriture de \(n\) en base \(b\), ainsi que les quelques puissances de \(a^j\) pour \(j<b\).

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    puissance_b(x,n,b) := {
      local a,r,xj,j;
      // decomposition de n en base b
      a := ecriture(n,b);
      r := length(a)-1;
      // calcul des a^i pour i < b
      xj := newList(b);
      xj[0] := 1; for(j:=1;j<b;j++){xj[j]:=xj[j-1]*x;}
      // calcul
      p := 1;
      for(j:=0;j<=r;j++){
        p := p^b;
        if(a[j]>0) p := p * xj[a[j]];
      }
      return p;
    }
    

    On obtient de meilleurs résultats sur les puissances de 2, pour ce calcul 256 est une bonne base.

Application aux suites récurrentes

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

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

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

  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}\). Quel est l’avantage par rapport aux méthodes employées précédemment ?

    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. Quel est le nombre de chiffres de \(F_{F_{F_{F_{F_4}}}}\) ? (utiliser l’expression algébrique exacte de \(F_n\)).

    solution xcas ⇲ ⇱ solution xcas

    Ce nombre est trop gros pour être calculé par xcas. En revanche on peut calculer le terme d’avant, \(n=F_{F_{F_{F_{4}}}}\).

    n := (fib_mat@@4)(4)
    

    Il reste à trouver le nombre de chiffres de \(F_n\), soit la valeur de \(\log_{10}F_n\).

    Soient \(\phi=\frac{1+\sqrt5}2\) et \(1-\phi\) les deux racines de \(x^2=x+1\), on a l’expression

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

    En passant au logarithme, on a \(\log F_n\sim n\log\phi\), et plus précisément

    \[\log_{10}F_n = n\log_{10}\phi - \log_{10}\sqrt{5} + \epsilon\]

    avec \(0<\epsilon=\log_{10}(1+(1-\phi)^n)<(1-\phi)^n\).

    phi := (1+sqrt(5))/2;
    r := n * logb(phi,10) - logb(5,10)/2;
    ceil(r);
    

Nombres premiers

  1. Créer la liste L des nombres entiers inférieurs à 100 en utilisant seq, select, et le test de primalité d’Xcas.

    solution xcas ⇲ ⇱ solution xcas

    il faut faire attention à donner une liste en mettant les crochets

    select(n->is_pseudoprime(n),[seq(n,n=2..100)])
    
  2. En déduire le nombre de nombres premiers inférieurs à 10000.

    solution xcas ⇲ ⇱ solution xcas
    length(select(n->is_pseudoprime(n),[seq(n,n=2..100))])
    

    ou bien:

    sum(is_pseudoprime(n),n,1,100)
    
  3. Faire de même en écrivant votre propre fonction de test de primalité est_premier.

On veut calculer plus rapidement la liste des premiers inférieurs à 10000 (ou leur nombre).

  1. L’algorithme de crible d’Eratosthène consiste à partir de la liste des entiers de \(2\) à \(n\), et en partant du premier \(p=2\) rayer tous ses multiples. Le nombre premier suivant est le premier entier \(p>2\) qui n’a pas été rayé, on raye ses multiples, et ainsi de suite.

    Le programmer, et mesurer le temps mis pour trouver les premiers inférieurs à 10000 (utiliser time).

    solution xcas ⇲ ⇱ solution xcas
     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    eratosthene(n):={
      local L,k,p;
      L := [seq(k,k=0..n)];
      L[1]:=0;
      for(p:=2;p<=sqrt(n);p++) {
        if(L[p]) {
          for(k:=p^2; k<=n; k+=p) L[k]:=0;
        }
      }
      select(k->k,L);
    }
    
    time(eratosthene(10000))
    
  2. Un autre algorithme consiste à tenir à jour une liste des nombres premiers déjà découverts, et à tester la primalité d’un entier en ne divisant que par ces nombres premiers.

    Le programmer, et mesurer le temps mis pour trouver les premiers inférieurs à 10000.

    solution xcas ⇲ ⇱ solution xcas
     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    liste_premiers(n):= {
      local k,P,j,p,np;
      P:=[]; np:= 0;
      for(k:=2;k<=n;k++) {
        p := 2;
        for(j:=0;p^2<=k;p:=P[j++]) {
          if(!irem(k,p)) break;
        }
        if(p^2>k) { P:=append(P,k); }
      }
      P;
    }
    

    c’est plus rapide

    time(liste_premiers(10000))
    

Quelques complexités usuelles

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})\)