Metoda Hornera wyliczania wartości wielomianów

Wielomian

Jak wiadomo, wielomian (stopnia $N$ ), to wyrażenie postaci: $w (x) = \sum_{i = 0}^{N} a_{i} x^{I}$

Wyliczanie wartości wielomianu

W naturalny sposób, kod (C) wyliczający jego wartość wyglądać może tak:

double wielomian (int N, double *a, double x)
{
  double w = 0.;
  int i;
  for (i = 0; i < N; i++)
    w += pow (x, i) * a[i];
  return w;
}

Sprytny student zauważy, że kod ten można „zoptymalizować" w sposób następujący:

double wielomian (int N, double *a, double x)
{
  double w = a[0];
  int i;
  for (i = 1; i < N; i++)
    w += pow (x, i) * a[i];
  return w;
}

oszczędzając jeden obieg pętli (i jedno wywołanie funkcji pow). Czy możliwe są dalsze optymalizacje? Oczywiście, można zrezygnować z funkcji pow:

double wielomian (int N, double *a, double x)
{
  double w = 0.;
  int i;
  double pow = 1;
  for (i = 0; i < N; i++)
    {
      w += pow * a[i];
      pow *= x;
    }
  return w;
}

(Prawdę mówiąc warto rezygnować z funkcji pow, gdyż jest to armata wystawiona na wróbelka: podnosimy $x$ do potęgi całkowitej.)

Czy są możliwe dalsze optymalizacje? Czy może to mieć wpływ na dokładność obliczeń? Czy jakąś własność numeryczną algorytmu?

„Tradycyjna" metoda obliczeń wymaga $1 + 2 + \dots + (N - 1) + N = \frac{N (N + 1)}{2}$ mnożeń i $N$ dodawań.

Metoda Hornera

Zapisując wielomian nieco inaczej: $w (x) = a_{0} + x (a_{1} + x (a_{2} + \dots + x (a_{N - 1} + x a_{N}) \dots))$ redukujemy liczbę operacji do $N$ dodawań i $N$ mnożeń.

Co więcej algorytm Hornera jest numerycznie poprawny.

Algorytm Hornera może być również przydatny podczas dzielenia wielomianu $w (y)$ przez dwumian $y - x$ dając w wyniku współczynniki wyniku (wielomianu będącego ilorazem $\frac{w (y)}{y - x}$ ).