Algorithmique de base, suite

Les 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.

    solution xcas ⇲ ⇱ solution xcas
    evalue_base(a,b) := {
      local(j);
      sum(v[j]*b^(j),j,0,length(v)-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 lsite 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):= {
      if(n==0) return 1;
      if(irem(n,2)) {
        return a*puissance_rec(a^2,iquo(n,2));
        }
      else {
        return puissance_rec(a^2,n/2);
        }
      }
    

    À 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_rec2(a,n):= {
      if(n==0) return 1;
      if(irem(n,2)) {
        return a*puissance_rec(a,iquo(n,2))^2;
        }
      else {
        return puissance_rec(a^2,n/2)^2;
        }
      }
    
  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

  1. Ecrire une fonction fib_rapide(n) qui calcule le \(n\)-ième terme de la suite de Fibonacci par puissance rapide de la matrice A:=[[0,1],[1,1]].

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

  2. 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);
    

Exponentiation modulaire

  1. Exemple à retenir : calculer les 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
    
  2. Déterminer les fonctions xcas qui effectuent les exponentiations modulaires \(a^k\bmod n\), ou \(P^n\bmod Q\) pour des polynômes \(P\) et \(Q\).

RSA

  1. Soit \(n=pq\)\(p\) et \(q\) sont des nombres premiers distincts. A quelle condition sur \(e\) la fonction \(f_e:a\to a^e\) est-elle une bijection de \(\Z/n\Z\) ?

    solution xcas ⇲ ⇱ solution xcas

    Il faut et il suffit que \(e\) soit premier à \(\phi(n)=(p-1)(q-1)\). Dans ce cas, l’inverse est \(f_d\)\(d\) est l’inverse de \(e\) modulo \(\phi(n)\).

    remarque: on vérifie que \(f_d\) et \(f_e\) sont inverses l’une de l’autre, cela découle du théorème de Lagrange, de sorte que la condition sur \(e\) est suffisante. Elle est aussi nécessaire si l’on considère un élément d’ordre un éventuel diviseur premier commun à \(e\) et \(\phi(n)\), lequel élément existe d’après le lemme de Cauchy (ou les théorèmes de Sylow): cet élément serait dans le noyau de \(f_e\).

  2. On pose \(e=65537\). Programmer efficacement l’application \(f_e\). Programmer également l’application inverse de manière efficace.

    solution xcas ⇲ ⇱ solution xcas

    \(e\) vaut \(2^{16}+1\), si bien que son écriture binaire est très simple. Il suffit d’élever 16 fois au carré, puis de multiplier par \(x\).

    fe(x,n) := {
      local j, r;
      r := x;
      for(j:=0;j<16;j++) r := irem(r^2,n);
      r := irem(r*x,n);
    }
    

    Pour \(f_d\), on se contente de la puissance rapide générique (voir question suivante).

  3. On pose p:=nextprime(2^255) et q:=nextprime(2^256). En utilisant les fonctions asc et char, ainsi que les fonctions de représentation des nombres en base \(b=256\), codez un message sous forme d’un entier inférieur à \(n=pq\), puis envoyez le chiffré par RSA à votre voisin et déchiffrez le sien

    solution xcas ⇲ ⇱ solution xcas

    On met en place le système RSA

    p := nextprime(2^255);
    q := nextprime(2^256);
    n := p*q;
    phin := (p-1)*(q-1);
    e := 2^16+1;
    d := ( (e mod phi)^-1 ) % 0
    fe(m) := powermod(m, e, n);
    fd(m) := powermod(m, d, n);
    

    Cornélius peut chiffrer le message suivant

    M := asc("Zephir va trahir Babar");
    m := evalue_base(M,256);
    c := fe(m)
    

    et Céleste le déchiffrer

    m := fd(c);
    M := ecriture(m, 256)
    char(M);
    
  4. Pour quelles raisons est-il d’usage de prendre cette valeur de \(e\) ?

    solution xcas ⇲ ⇱ solution xcas

    \(e\) est un nombre de Fermat (d’où une exponentiation rapide très simple) et le plus grand (connu) qui soit premier (ce qui diminue les chances qu’il soit non premier à \(\varphi(n)\)).

Symbole quadratique

  1. A l’aide de la commande factor, déterminer un polynôme \(P\in\F_{101}[X]\) irréductible de degré 4.

    solution xcas ⇲ ⇱ solution xcas

    On cherche un polynôme au hasard, et on essaie de le factoriser.

    p := 101;
    P := ( x^4+rand(p)*x+rand(p) ) % p;
    factor(P);
    
  2. On note \(K\) le corps fini \(\F_{101}[X]/(P)\). Quel est le cardinal de \(K^\ast\) ? Quel est le cardinal de \(\{u^2, u\in K^\ast\}\) ?

    solution xcas ⇲ ⇱ solution xcas

    On pose \(p=101\) : le cardinal de \(K\) vaut \(p^4\), donc \(\card{ K^\ast}=p^4-1\). Les carrés sont l’image de l’application carré, qui a pour noyau \(\pm1\), on en a donc \(\frac{p^4-1}2\).

  3. En déduire un test, à l’aide d’un calcul de puissance rapide, déterminant si la classe de \(x\) est un carré dans \(K\).

    solution xcas ⇲ ⇱ solution xcas

    Les carrés forment un sous-groupe, donc d’après le théorème de Lagrange tout carré \(c\) vérifie \(c^{\frac{p^4-1}2}=1\). Réciproquement, cette équation est polynomiale et n’admet pas plus de solutions que son degré, donc elle caractérise exactement les carrés. Pour les \(x\) non-carrés, on a \(r=x^{\frac{p^4-1}2}=-1\) (car \(r^2=1\) et \(r\neq 1\)).

    carremodP(x) := {
      return powermod(x,(p^4-1)/2,P)==1;
      }
    

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