n-мерный полином, способ задания и его производная

В программировании довольно часто приходится работать с полиномами. Будь то простое квадратное уравнение или многочлен n-й степени его все-равно нужно как-то програмно задать. Рассмотрим один из способов задания полиномов:

Общий вид: P_n(x)=a[n]*xⁿ+a[n-1]*x^n-1+…+a[1]*x+a[0]

Если вы захотите его использовать в этом виде то ничего хорошего у вас не выйдет т.к. прийдется использовать функцию возведения в степень, что в свою очередь потянет за собой огромное кол-во умножений (2n-1), а операция умножения, как известно, далеко не самая быстрая.
Но это еще не все… Как говорит теория ошибок, погрешность результата прямопропорциональна кол-ву сомножителей т.е. чем больше вы делаете умножений, тем выше погрешность результата.
Что же делать? Как уменьшить число умножений? На самом деле все просто. Вспомним школьную математику и вынесем x за скобки:

Pn(x)=(…(a[n]*x+a[n-1])*x+…+a[1])*x+a[0]

Просмотрев это выражение видно, что мы избавились от операции возведения в степень и сократили число умножений до n раз, а следовательно, повыслили точность и скорость вычислений.
Как же его посчитать? Тут тоже все просто: нужно представить коэффициенты полинома в виде массива и организовать примитивный цикл:

     p:=0;
     p:=p*x+a[n];
     p:=p*x+a[n-1];
     . . .
     p:=p*x+a[1];
     p:=p*x+a[0];

Посмотрим как это можно реализовать на Паскале:

const
  anMax = 3;
  An: array[0..anMax] of double = (-5.372,1.2493,0.559,-0.13);
{ Тут описан следующий полином: Pn(x)=-5.372 + 1.2493*x + 0.559*x^2 - 0.13*x^3 }
function Pn(x: double): double;
var
  i: byte;
  Res: double;
begin
  Res:=0;
  for i:=anMax downto 0 do
    Res:=Res*x+An[i];
  Pn:=Res;
end;

Это конечно все хорошо, скажете вы, но что-то все равно это не «греет».
Посмотрим теперь пример, как можно легко и просто посчитать n-ю производную заданного полинома.

Вспомним правила вычисления производных:

1-я (xⁿ)' = n*x^(n-1) 2-я (xⁿ) = n*(n-1)*x(n-2) 3-я (xn)' = n*(n-1)*(n-2)*x^(n-3) . . . i-я (xⁿ)⁽ⁱ⁾ = n*(n-1)*(n-2)*…*(n-i+1)*x^(n-i)

Или на конкретном примере:

Pn(x) = a[0] + a[1]*x + a[2]*x² + a[3]*x³ Pn'(x) = a[1] + a[2]*2*x + a[3]*3*x² Pn(x) = a[2]*2+ a[3]*3*2*x Pn'(x) = a[3]*3*2 Pn(x) = 0

Как видно из примера, чем выше степень производной, тем меньше остается сомножителей и тем меньше мы используем коффициентов из массива, сдвигаясь каждый раз на позицию вправо. Кроме того, не сложно заметить, что ростом степени производной падает степень и увеличивается коэффициен перед x.
Если со степенью все просто и понятно с первого взгляда: из начальной степени вычитаем степень производной и если результат меньше нуля, то выбрасываем этот сомножитель вместе с его коэффициентом.
А как же расчитать коэффициент? Если немного присмотреться, то можно заметить что коэффициет равен частному двух факториалов:

k = a! / (a-i)! , где:

k - коэффициент a - начальная степень i - степень производной

Таким образом, получаем формулу:

i-я (xⁿ)⁽ⁱ⁾ = n!/(n-i)! * x^(n-i)

const
  anMax = 3;
  An: array[0..anMax] of double = (-5.372,1.2493,0.559,-0.13);
 
function Func(x: double; d: byte):double; far;
{ где d - степень производной }
var
  Res: double;
  i: byte;
function frac(a,b: byte): longint;
  var
    res: longint;
  begin
    res:= 1;
    while a > b do
      begin
        res:= res*a;
        dec(a);
      end;
    frac:= res;
  end;
begin
  Res:=0;
  for i:= anMax-d downto 0 do
    Res:=res*x+An[i+d]*frac(i+d,i);
  Func:= Res;
end;