Variations sur le schéma de Horner

Partie I - Évaluation d'un polynôme

Question 1 - Évaluation itérative d'un polynôme

L'algorithme de Horner itératif parcourt les coefficients du polynôme selon les puissances décroissantes. À chaque étape, on multiplie le résultat temporaire par t et on y ajoute le coefficient suivant.

Cet algorithme nécessite n multiplications et n additions (où n est le degré du polynôme).

evalh := proc(P::array(numeric), t::numeric)
  local y, k;
  y := 0;
  for k to size(P) do 
    y := y * t + P[k];
  od;
  RETURN(y);
end:

Question 2 - Évaluation récursive d'un polynôme

La formulation récursive exploite le fait que P(t) = t * Q(t) + a0, où Q est le polynôme quotient de P par X. On transmet la longueur du sous-tableau pour traiter les coefficients restants.

evalhr2 := proc(P::array(numeric), n::integer, t::numeric)
  if n <= 1 then 
    RETURN(P[1]);
  else 
    RETURN(evalhr2(P, n - 1, t) * t + P[n]);
  fi;
end:

Question 3 - Évaluation en un point complexe (méthode directe)

En séparant les parties réelle et imaginaire, on évalue P(x + iy) = u + iv. À chaque itération de la boucle de Horner, la multiplication par z = x + iy donne une nouvelle partie réelle (ux - vy) et une nouvelle partie imaginaire (uy + vx).

Cette méthode requiert 3n additions et 4n multiplications.

evalhc := proc(P::array(numeric), x::numeric, y::numeric)
  local u, v, k, t_val; 
  u := 0; 
  v := 0;
  for k to size(P) do
    t_val := u * x - v * y + P[k];
    v := u * y + v * x;
    u := t_val;
  od;
  RETURN(array([u, v]));
end:

Question 4 - Évaluation complexe optimisée (division euclidienne)

Pour réduire le nombre d'opérations, on utilise la division euclidienne par le trinôme réel (X - z)(X - z̅) = X² - 2x X + (x² + y²). Le reste de cette division donne la valeur cherchée avec seulement 2n + 3 additions et 2n + 4 multiplications.

evalhc2 := proc(P::array(numeric), x::numeric, y::numeric)
  local r, m, b, c, k, t_val;
  r := 2 * x; 
  m := x * x + y * y; 
  b := 0; 
  c := 0;
  for k to size(P) do 
    t_val := P[k] + r * b - m * c; 
    c := b; 
    b := t_val; 
  od;
  RETURN(array([b - c * x, c * y]));
end:

Partie II - Division synthétique

Question 1 - Quotient par un monôme

La division de P par (X - t) permet de trouver le quotient Q. L'algorithme est identique à l'évaluation, mais on mémorise les valeurs intermédiaires.

quo1 := proc(P::array(numeric), t::numeric)
  local n, Q, k; 
  n := size(P);
  if n = 1 then RETURN(array([0])) fi;
  Q := array(1..n-1); 
  Q[1] := P[1];
  for k from 2 to n-1 do
    Q[k] := Q[k-1] * t + P[k];
  od; 
  RETURN(eval(Q));
end:

Question 2 - Division avec modification sur place

Pour économiser de la mémoire, on peut inscrire directement les coefficients du quotient (et le reste final) dans le tableau initial P. L'argument n précise la portée de l'opération.

Quo1 := proc(P::array(numeric), n::integer, t::numeric)
  local k;
  for k from 2 to n do
    P[k] := P[k-1] * t + P[k];
  od;
  RETURN();
end:

Question 3 - Quotient par un trinôme

On divise P par (X² + αX + β). Les coefficients du quotient et le reste de degré 1 remplacent progressivement les valeurs du tableau P.

Quo2 := proc(P::array(numeric), a::numeric, b::numeric)
  local m, u, v, k, t_val;
  m := size(P); 
  u := 0; 
  v := 0;
  for k to m do
    t_val := P[k] - a * u - b * v;
    v := u;
    u := P[k];
    P[k] := t_val;
  od;
  P[m] := u + a * v;
  RETURN();
end:

Partie III - Translatés d'un polynôme, dérivées successives

Question 1 - Translation avec algorithme de Horner répété (sur place)

Développer P(X + t) revient à exprimer P dans la base des (X - t)^k. On utilise des divisions successives par (X - t) via Quo1.

Translat := proc(P::array(numeric), t::numeric)
  local k;
  for k from size(P) to 2 by -1 do
    Quo1(P, k, t);
  od;
  RETURN();
end:

Question 2 - Translation renvoyant un nouveau tableau

Cette version crée une copie pour ne pas altérer le polynôme d'origine et effectue directement les calculs.

translat := proc(P::array(numeric), t::numeric)
  local n, Q, j, k;
  n := size(P); 
  Q := copy(P);
  for k from n to 2 by -1 do
    for j from 2 to k do
      Q[j] := Q[j-1] * t + Q[j];
    od;
  od;
  RETURN(eval(Q));
end:

Question 3 - Translation optimisée avec homothéties

Cette méthode réduit drastiquement les multiplications en passant par Q1(X) = P(tX), Q2(X) = Q1(X+1), puis Q(X) = Q2(X/t).

translat2 := proc(P::array(numeric), t::numeric)
  local n, Q, T, j, k;
  if t = 0 then RETURN(eval(P)) fi;
  n := size(P); 
  Q := copy(P); 
  T := array(1..n); 
  T[n] := 1;
  for k from n-1 to 1 by -1 do
    T[k] := t * T[k+1]; 
    Q[k] := Q[k] * T[k];
  od;
  for k from n to 2 by -1 do
    for j from 2 to k do 
      Q[j] := Q[j-1] + Q[j]; 
    od;
  od;
  for k from n-1 to 1 by -1 do
    Q[k] := Q[k] / T[k];
  od;
  RETURN(eval(Q));
end:

Question 4 - Calcul des dérivées successives

Les coefficients obtenus par translation donnent directement P^(k)(t) / k!. Il suffit de multiplier par les factorielles successives.

derivs := proc(P::array(numeric), t::numeric)
  local n, Q, k, f;
  n := size(P); 
  f := 1; 
  Q := translat(P, t);
  for k from 2 to n-1 do
    f := f * k;
    Q[n-k] := f * Q[n-k];
  od;
  RETURN(eval(Q));
end:

Partie IV - Règle des signes de Descartes

Question 1 - Nombre de changements de signe

On parcourt les coefficients non nuls pour repérer les alternances de signe.

varsign := proc(P::array(numeric))
  local k, s, f;
  f := 0; 
  s := 0;
  for k to size(P) do 
    if f * P[k] < 0 then
      s := s + 1;
    fi;
    if P[k] <> 0 then
      f := P[k];
    fi;
  od; 
  RETURN(s);
end:

Question 2 - Majorant des racines supérieures à a

En étudiant les changements de signe de P(X + a), on majore le nombre de racines réelles strictement supérieures à a.

rootsup := proc(P::array(numeric), a::numeric)
  RETURN(varsign(translat(P, a)));
end:

Question 3 - Majorant des racines inférieures à a

On étudie les racines de P(-X) supérieures à -a.

rootinf := proc(P::array(numeric), a::numeric)
  local Q, k; 
  Q := copy(P);
  for k to size(Q) by 2 do 
    Q[k] := -Q[k];
  od;
  RETURN(rootsup(Q, -a));
end:

Partie V - Méthode de Newton de résolution de P(x) = 0

Question 1 - Convergence classique

On utilise la formule x_{n+1} = x_n - P(x_n) / P'(x_n). On peut calculer P et P' simultanément.

newton := proc(P::array(numeric), x::numeric)
  local n, y, z, k;
  n := size(P); 
  z := 0; 
  y := 0;
  for k to n do
    z := x * z + y; 
    y := x * y + P[k];
  od;
  RETURN(x - y / z);
end:

Question 2 - Convergence quadratique accélérée (racines multiples)

Pour forcer une convergence quadratique même sur les racines multiples, on utilise P(x)P'(x) / (P'(x)² - P(x)P''(x)). Un triple algorithme de Horner calcule P, P' et P''/2.

newton2 := proc(P::array(numeric), x::numeric)
  local n, y, z, t_val, k;
  n := size(P);
  y := 0; z := 0; t_val := 0;
  for k to n do
    t_val := x * t_val + z;
    z := x * z + y;
    y := x * y + P[k];
  od;
  RETURN(evalf(x - y * z / (z^2 - 2 * y * t_val)));
end:

Partie VI - Forme de Newton du polynôme interpolateur

Question 1 - Différences divisées récursives

Cette fonction calcule récursivement d_{i,j}, le coefficient correspondant à l'interpolation.

dij := proc(X, Y, i, j)
  if i = j then 
    RETURN(Y[i]);
  else 
    RETURN((dij(X, Y, i+1, j) - dij(X, Y, i, j-1)) / (X[j] - X[i]));
  fi;
end:

Question 2 - Création de la forme de Newton (avec récursion)

On calcule la première ligne du tableau des différences divisées.

pnewton := proc(X, Y)
  local n, P, k;
  n := size(X); 
  P := array(1..n);
  for k to n do
    P[k] := dij(X, Y, 1, k);
  od;
  RETURN(eval(P));
end:

Question 3 - Calcul itératif de la forme de Newton

Pour éviter la redondance des calculs, on actualise le vecteur Y en place, de droite à gauche.

PNewton := proc(X, Y)
  local n, j, k;
  n := size(X);
  for k to n-1 do
    for j from n by -1 to k+1 do
      Y[j] := (Y[j] - Y[j-1]) / (X[j] - X[j-k]);
    od; 
  od; 
  RETURN();
end:

Question 4 - Ajout d'un point à l'interpolation

La forme de Newton permet d'ajouter un point sans tout recalculer, en prolongeant le tableau des différences divisées.

newpoint := proc(X, P, x_val, y_val)
  local n, Q, k; 
  n := size(P);
  Q := array(1..n+1); 
  Q[1] := y_val;
  for k to n do
    Q[k+1] := (P[k] - Q[k]) / (X[k] - x_val);
  od;
  RETURN(eval(Q));
end:

Question 5 - Évaluation du polynôme interpolateur

L'évaluation s'apparente à l'algorithme de Horner, en remplaçant la multiplication par x par une multiplication par (x - X[k]).

evalnewton := proc(X, P, x_val)
  local y_val, k; 
  y_val := 0;
  for k from size(P) by -1 to 1 do 
    y_val := P[k] + (x_val - X[k]) * y_val; 
  od;
  RETURN(y_val);
end:

Partie VII - Autour du théorème chinois

Question 1 - Coefficients de Bezout

a) Version récursive

bezoutr := proc(m::integer, n::integer)
  local C, q, r;
  if n = 0 then 
    C := array([1, 0, m]);
  else
    q := iquo(m, n); 
    r := irem(m, n);
    C := bezoutr(n, r);
    C := array([C[2], C[1] - q * C[2], C[3]]);
  fi;
  RETURN(eval(C));
end:

b) Version itérative

bezout2 := proc(m::integer, n::integer)
  local a, b, C1, C2, q, r_val, T, k;
  a := m; b := n;
  C1 := array([1, 0, m]);
  C2 := array([0, 1, n]);
  while b <> 0 do
    q := iquo(a, b); 
    r_val := irem(a, b);
    a := b; b := r_val;
    T := copy(C1); 
    C1 := copy(C2);
    for k to 3 do
      C2[k] := T[k] - q * C1[k];
    od;
  od;
  RETURN(eval(C1));
end:

Question 2 - Solution d'un système de congruences

a) Solution directe (Modulos globaux) (D'après le corrigé partiel, on calcule un produit Mj et on somme.)

b) Solution optimisée par schéma de Horner (Version finale réparée)

Le second algorithme évite les entiers géants en gardant les calculs confinés dans chaque module successif. Voici le code complété :

chinese := proc(A::array(integer), N::array(integer))
  local r, U, i, j, B, x;
  r := size(N);
  U := array(1..r, 1..r);
  for i to r do 
    for j to i-1 do
      B := bezout(N[i], N[j]);
      if B[3] <> 1 then error "pgcd<>1"; fi;
      U[i,j] := B[1];
      U[j,i] := B[2];
    od; 
  od;
  B := copy(A);
  for i from 2 to r do
    for j from 1 to i-1 do
      B[i] := ((B[i] - B[j]) * U[j,i]) mod N[i];
    od;
  od;
  x := B[r];
  for i from r-1 by -1 to 1 do
    x := x * N[i] + B[i];
  od;
  RETURN(x);
end:

Partie VIII - Algorithmes "compte-gouttes"

Question 1 - Numération à base variable adaptée à "e"

a) Preuve et algorithme de base Puisque b_k <= k - 1 et qu'il existe un indice m avec b_m < m - 1, la somme r_n = Σ(k=n+1 à ∞) b_k / k! est strictement majorée par Σ(k=n+1 à ∞) (k-1) / k!. Cette somme vaut Σ (1/(k-1)! - 1/k!), qui se télescope exactement à 1/n!. Donc 0 ≤ r_n < 1/n!. En multipliant par n!, on a n! x = n! x_n + n! r_n. L'entier n! x_n correspond à la partie entière [n! x], d'où b_n = [n! x] mod n.

c) Fonction todec

todec := proc(B, n)
  local x_val, k;
  x_val := B[n];
  for k from n-1 by -1 to 1 do
    x_val := B[k] + x_val / (k + 1);
  od;
  RETURN(x_val);
end:

d) Fonction tovar

tovar := proc(xn, n)
  local B, k, y;
  B := array(1..n);
  y := xn;
  B[1] := floor(y);
  for k from 2 to n do
    y := (y - B[k-1]) * k;
    B[k] := floor(y);
  od;
  RETURN(eval(B));
end:

Question 2 - Le calcul des décimales du nombre e

a et b) Divisions successives avec report

La multiplication par 10^m décale le nombre, et les restes successifs forment les nouveaux chiffres de la base.

goutte := proc(X::array(numeric), m::integer)
  local n, k, c, p, carry;
  n := size(X);
  p := 10^m;
  for k from 1 to n do 
    X[k] := X[k] * p; 
  od;
  carry := 0;
  for k from n by -1 to 2 do
    c := X[k] + carry;
    X[k] := c mod k;
    carry := iquo(c, k);
  od;
  X[1] := X[1] + carry;
  RETURN();
end:

Question 3 - Le calcul des décimales de π

b) Évaluation (todec2) Le rapport u_k / u_{k-1} vaut k / (2k + 1). On adapte le schéma de Horner :

todec2 := proc(X)
  local n, k, x_val;
  n := size(X);
  x_val := X[n];
  for k from n-1 by -1 to 1 do
    x_val := X[k] + x_val * k / (2 * k + 1);
  od;
  RETURN(x_val);
end:

d) Algorithme goutte2 pour π Le poids de base n'est plus k! mais dérive de (2k+1). La retenue à l'étape k doit être distribuée avec un coefficient (k-1).

goutte2 := proc(X::array(numeric), m::integer)
  local n, k, c, q, r_val, p;
  n := size(X);
  p := 10^m;
  for k from 1 to n do X[k] := X[k] * p od;
  for k from n by -1 to 2 do
    c := X[k];
    q := iquo(c, 2 * k - 1);
    r_val := irem(c, 2 * k - 1);
    X[k] := r_val;
    X[k-1] := X[k-1] + q * (k - 1);
  od;
  RETURN();
end:

e) Majoration de l'erreur La série des restes vérifie u_k < 1/2 u_{k-1}. La somme géométrique majorante donne π - π_m < u_m Σ (1/2)^i = 2 u_m. Les chiffres calculés à l'indice m assurent que l'on obtient l'approximation décimale désirée.

Méthode

Face à ce type de sujet mêlant informatique et mathématiques :

1. Tracez toujours l'état de la mémoire sur papier pour les algorithmes itératifs complexes (en particulier ceux avec les "reports" comme les algorithmes compte-gouttes).
2. Validez le concept avec le schéma de Horner : l'évaluation naïve O(n²) doit presque toujours être remplacée par O(n) additions et multiplications successives.
3. Repérez les invariances mathématiques : une translation P(X+t) ou une division P(X)/(X²-zX+m) est souvent un "Horner masqué" qui manipule la structure du polynôme.
4. Pour Maple (ou tout langage), traitez avec soin la manipulation en place (Quo1) contre la création de copies (copy(P)). Le concours pénalise lourdement l'écrasement de données lorsqu'une valeur de retour isolée est attendue.