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+\cdots+(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 +\cdots+x(a_{N-1}+ x a_N)\ldots))$$ 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}$).