Вычисление значения многочлена в точке является одной из простейших классических задач программирования.
При проведении различного рода вычислений часто приходится определять значения многочленов при заданных значениях аргументов. Часто приближенное вычисление функций сводится к вычислению аппроксимирующих многочленов.
Рядового читателя Хабрахабр нельзя назвать неискушенным в применении всяческих извращений. Каждый второй скажет, что многочлен надо вычислять по правилу Горнера. Но всегда есть маленькое «но», всегда ли схема Горнера является самой эффективной?


Я не ставлю цель точно описать алгоритмы для вычисления многочленов, а лишь показать, что в некоторых случаях можно (нужно) применять схемы отличные правила Горнера. Для тех, кого заинтересует материал, в конце статьи приведен список литературы, с которой можно ознакомиться для более детального изучения вопроса.
Кроме того, иногда становиться обидно, что фамилии наших русских математиков остаются малоизвестными. К тому же мне просто приятно рассказать о работах наших математиков.

Схема Горнера


При вычислении значений многочленов очень широкое применение получило правило Горнера. Метод назван в честь британского математика Уильяма Джорджа Горнера.
В соответствии с этим правилом многочлен n-й степени:
{{P}_{n}}(x)={{a}_{0}}{{x}^{n}}+{{a}_{1}}{{x}^{n-1}}+...+{{a}_{n-1}}x+{{a}_{n}}

представляется в виде
{{P}_{n}}(x)=(...(({{a}_{0}}x+{{a}_{1}})x+{{a}_{2}})x+...+{{a}_{n-1}})x+{{a}_{n}}.

Вычисление значения многочлена производится в порядке, определяемом скобками. Что имеем? Чтобы вычислить многочлен {{P}_{n}}(x) по схеме Горнера, надо выполнить n умножений и n-k сложений (здесь k – число коэффициентов многочлена, равных 0). Если {{a}_{0}}=1, то умножений будет n-1.
Можно показать, что для вычисления многочленов, общего вида нельзя построить схему более экономичную по числу операций, чем схема Горнера.
Самая большая привлекательность схемы Горнера состоит в простоте алгоритма для вычисления значения многочлена.

Исключения


При вычислении многочленов специального вида может потребоваться меньшее число операций, чем при применении универсальной схемы Горнера. Например, вычисление степени {{x}^{n}} по схеме Горнера означает последовательное перемножение n множителей и требует n-1 умножение. Однако каждый первый читатель скажет, что для вычисления, например, {{x}^{8}} нужно последовательно вычислить {{x}^{2}}=x\cdot x, {{x}^{4}}={{x}^{2}}\cdot {{x}^{2}}, {{x}^{8}}={{x}^{4}}\cdot {{x}^{4}}, т.е. выполнить всего 3 умножения вместо 7.

А есть что-то еще, ведь схема Горнера самая экономичная?


На самом деле все решают объемы вычислений. Если надо вычислить одно значение многочлена, то лучше схемы Горнера ничего не придумано. Но если значения многочлена вычисляются во многих точках, то появляется возможность сэкономить большое число операций умножения за счет предварительных вычислений, выполняемых ровно один раз. Это может значительно ускорить работу программы.

В некоторых случаях для получения значений полиномов целесообразно использовать двухэтапные схемы. На первом этапе выполняются действия только над коэффициентами многочлена, он преобразуется к специальному виду. На втором же этапе вычисляют уже значение самого многочлена при заданных значениях аргумента. При этом может оказаться, что количество операций, выполняемых на втором этапе будет меньше, чем при вычислениях по схеме Горнера.

Снова замечу, что такие методы вычислений целесообразны при вычислении значений многочлена {{P}_{n}}(x) для большого числа значений x. Выигрыш получается, за счет того, что первый этап для многочлена выполняется лишь один раз. Примером может послужить вычисление элементарных функций, где приближающий многочлен готовиться заранее.

В дальнейших рассуждениях, говоря о количестве операций для вычисления {{P}_{n}}(x), я буду иметь в виду сложность второго этапа вычислений.

Схема Дж.Тодта для многочленов 6 степени


Имеем следующий многочлен:

{{P}_{6}}(x)={{x}^{6}}+A{{x}^{5}}+B{{x}^{4}}+C{{x}^{3}}+D{{x}^{2}}+Ex+F.

Для вычислений используем следующие вспомогательные многочлены:

{{p}_{1}}(x)=x(x+{{b}_{1}}),

{{p}_{2}}(x)=({{p}_{1}}+x+{{b}_{2}})({{p}_{1}}+{{b}_{3}}),

{{p}_{3}}(x)=({{p}_{2}}+{{b}_{4}})({{p}_{1}}+{{b}_{5}})+{{b}_{6}}.

Коэффициенты {{b}_{j}} определяются методом неопределенных коэффициентов исходя из условия {{p}_{3}}(x)\equiv {{P}_{6}}(x). Из последнего условия составляем систему уравнений, приравнивая коэффициенты при равных степенях многочленов.

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

Говорить об универсальности данной схемы не приходится, но зато читатель наглядно может оценить уменьшение числа операций по сравнению со схемой Горнера.

Схема Ю.Л. Кеткова


Наконец-то, добрался и до наших математиков.

Ю.Л. Кетков дал общее представление многочлена n-й степени для n>5, всегда приводящее к действительным выражениям и требующее для вычисления многочлене n-й степени выполнения [(n+1)/2]+[n/4] умножений и n+1 сложений.

Например, при n=2k схема Кеткова сводится к нахождению многочленов:

{{N}_{2}}(x)=x({{b}_{0}}+x),

{{N}_{4}}(x)=({{N}_{2}}+{{b}_{1}}+x)({{N}_{2}}+{{b}_{2}})+{{b}_{3}},

{{N}_{6}}(x)={{N}_{2}}{{N}_{4}}+{{b}_{4}}x+{{b}_{5}},

\cdots\cdots\cdots

{{N}_{2k}}(x)=({{N}_{2}}+{{s}_{k}}{{b}_{2k-2}}){{N}_{2k-2}}+{{r}_{k}}{{b}_{2k-2}}x+{{b}_{2k-1}},

где {{r}_{k}}=0, {{s}_{k}}=1 при k –четном, и {{r}_{k}}=1, {{s}_{k}}=0, если k нечетное (k>2).

Все неизвестные коэффициенты находятся из равенства {{P}_{n}}(x)={{N}_{2k}}. В работах Кеткова для решения получающихся систем дается метод, дающий всегда действительные коэффициенты {{b}_{j}}.

Схемы В.Я. Пана


Э. Белага в своих работах дал строгое доказательство невозможности построения схемы вычисления произвольных многочленов n-й степени, использующей на втором этапе меньше, чем [(n+1)/2]+1 умножений и n сложений.

В.Я. Пан занимался вопросами оптимального вычисления многочленов. В частности, им предложено несколько схем для вычисления действительных многочленов, которые весьма близко подобрались к оценкам Э. Белаги. Приведу некоторые схемы Пана для действительных многочленов.
1. Схема для вычисления многочленов четвертой степени.
Рассматривается многочлен {{P}_{4}}(x)={{a}_{0}}{{x}^{4}}+{{a}_{1}}{{x}^{3}}+{{a}_{2}}{{x}^{2}}+{{a}_{3}}x+{{a}_{4}}.

Представим {{P}_{4}}(x) в виде:

g(x)=x(x+{{\lambda }_{1}}),

{{P}_{4}}(x)\equiv {{a}_{0}}\left( \left( g(x)+{{\lambda }_{2}} \right)\left( g(x)+x+{{\lambda }_{3}} \right)+{{\lambda }_{4}} \right),

где

{{\lambda }_{1}}=\frac{{{a}_{1}}-{{a}_{0}}}{2{{a}_{0}}},{{\lambda }_{2}}=\frac{{{a}_{3}}}{{{a}_{0}}}-{{\lambda }_{1}}\frac{{{a}_{2}}}{{{a}_{0}}}+({{\lambda }_{1}}+1)\lambda _{1}^{2},

{{\lambda }_{3}}=\frac{{{a}_{2}}}{{{a}_{0}}}-{{\lambda }_{1}}({{\lambda }_{1}}+1)-{{\lambda }_{2}},{{\lambda }_{4}}=\frac{{{a}_{4}}}{{{a}_{0}}}-{{\lambda }_{2}}{{\lambda }_{3}}.

2. Схема для вычисления {{P}_{n}}(x), n\ge 5.
Строим вспомогательные многочлены g(x), h(x), {{p}_{s}}(x):

g(x)=x(x+{{\lambda }_{1}}),h(x)=g(x)+x,{{p}_{0}}(x)=x,

{{p}_{s}}(x)={{p}_{s-1}}(x)\left( (g(x)+{{\lambda }_{4s-2}})(h(x)+{{\lambda }_{4s-1}})+{{\lambda }_{4s}} \right)+{{\lambda }_{4s+1}}, s=1,2,…,k.

Для вычисления значения многочлена используем выражения:

{{P}_{n}}(x)\equiv {{a}_{0}}{{p}_{k}}(x), при n=4k+1,

{{P}_{n}}(x)\equiv {{a}_{0}}x{{p}_{k}}(x)+{{\lambda }_{n}}, при n=4k+2,

{{P}_{n}}(x)\equiv {{a}_{0}}\left( {{p}_{k}}(x)(g(x)+{{\lambda }_{n-1}})+{{\lambda }_{n}} \right), при n=4k+3,

{{P}_{n}}(x)\equiv {{a}_{0}}x\left( {{p}_{k}}(x)(g(x)+{{\lambda }_{n-2}})+{{\lambda }_{n-1}} \right)+{{\lambda }_{n}}, при n=4k+4,

Эта схема на втором этапе требует [n/2]+2 умножения и n+1 сложения.

Особенностью данной схемы является то, что коэффициенты {{\lambda }_{j}} всегда существуют при n\ge 5 и действительных коэффициентах исходного многочлена.

У В.Я. Пана существуют и другие схемы для вычисления многочленов, в том числе и для комплексных.

Заключение


Резюмируя сказанное, замечу, что вычисление одного или нескольких значений полинома бесспорно нужно проводить с использованием схемы Горнера.

Однако, если число значений полинома, которые потребуется вычислить велико, а производительность очень важна, то имеет смысл рассмотреть применение специальных методов вычисления многочленов.

Некоторые читатели скажут, что возиться с применением схем, отличных от схемы Горнера, сложно, муторно и не стоит с этим связываться. Однако в реальной жизни встречаются задачи, в которых требуется вычислять просто огромное число значений многочленов с большими степенями (например, на их вычисление могут уходить месяцы), и уменьшение числа умножений в два раза даст существенный выигрыш во времени, даже если вам придется потратить пару дней на реализацию конкретной схемы для вычисления многочленов.

Литература


  1. Кетков Ю.Л. Об одном способе вычисления полиномов на математических машинах. // Известия ВУЗ'ов. Радиофизика, т.1., № 4, 1958
  2. В. Я. Пан, “Вычисление многочленов по схемам с предварительной обработкой коэффициентов и программа автоматического нахождения параметров”, Ж. вычисл. матем. и матем. физ., 2:1 (1962), 133–140
  3. В. Я. Пан, “О способах вычисления значений многочленов”, УМН, 21:1(127) (1966), 103–134
  4. В. Я. Пан, “О вычислении многочленов пятой и седьмой степени с вещественными коэффициентами”, Ж. вычисл. матем. и матем. физ., 5:1 (1965), 116–118
  5. Пан В. Я. Некоторые схемы для вычисления значений полиномов с вещественными коэффициентами. Проблемы кибернетики. Вып. 5. М.: Наука, 1961, 17–29.
  6. Белага Э. Г. О вычислении значений многочлена от одного переменного с предварительной обработкой коэффициентов. Проблемы кибернетики. Вып. 5. М.: Физматгиз, 1961, 7–15.

Комментарии (10)


  1. Uranix
    12.03.2016 13:45
    +2

    А как обстоит дело с векторизованными операциями, когда одновременно можно проводить несколько умножений и сложений?


  1. Moric
    12.03.2016 16:00
    +6

    На практике ещё нужно учитывать устойчивость схемы вычисления к ошибкам округления. Иначе можно получить быстрый, но не точный алгоритм вычисления полинома.


  1. Mercury13
    12.03.2016 20:22
    +2

    Можно вопрос по поводу спирта? Мне удалось записать вот такую шутку с нашего факультета.

    В 1999-2000 учебном году ему [преподавателю] впервые поручили вести лекции. Он вёл их препаршиво. На пресс-конференции, связанной с днём факультета:
    Студент: Вы видели когда-нибудь компьютер?
    И.И.: Конечно. Вот я помню машину XX-XX, 2-го поколения. Я бы её не променял даже на вагон персоналок.
    Вмешивается декан: Ведь эта машина работала на спирту.

    Что это за машина и почему «работала на спирту»?


    1. Mrrl
      12.03.2016 22:50
      +1

      Судя по сигнатуре, EC-40. Хотя могла быть и СМ-1420. Спирт использовался для протирки поверхностей накопителя на магнитных дисках и головок самого накопителя.


      1. 32bit_me
        13.03.2016 10:39

        У нас на факультете тоже была СМ-ка. Спирт на неё выдавался (начало 1990-х) и использовался для очевидных целей. Поэтому её и не списывали долго, наверное.


    1. aGRa
      13.03.2016 06:23
      +2

      Я припоминаю компьютер Искра 1030М Курского завода Счётмаш, на котором в детстве писал свои первые программы на Бейсике и Паскале. К нему прилагалось руководство, в котором, в лучших советских традициях, указывались и нормы расходных материалов, в частности, спирта (на протирку контактов, видимо). Так вот, количество потребного спирта составляло 2,5 литра в месяц. Это на персональный компьютер размером чуть больше современного системника. Естественно, в руководстве на более новые модели (выпущенные после того, как перепились все ответственные за обслуживание 1030М — например, 1031) норма расхода спирта была снижена более чем на порядок. Так что на спирту работали скорее не машины, а их операторы.


  1. MichaelBorisov
    12.03.2016 23:26
    +2

    Спасибо, автор, за интересный обзор.

    А вот я еще недавно наткнулся на следующую схему вычисления полиномов, не требующую умножений, если полином требуется вычислить в равноотстоящих точках.

    Идея — в составлении системы разностных уравнений и их решении в прямом времени. Для этого требуется только сложение. Так, для полинома первой степени y = a1x+a0
    y(x+dx)-y(x) = a1
    dx, где dx — интервал между точками, в которых требуется вычислить полином. Так как a1 и dx — константы, то умножение требуется выполнять только один раз.

    Для многочлена второй степени y0 = a2x^2 + a1x + a0 имеем:
    y0(x+dx)-y0(x) = a2(x+dx)^2 + a1(x+dx) — a2x^2 — a1x = 2a2xdx + a2dx^2 + a1dx
    Здесь в правой части имеется константа и функция 2
    a2xdx. Но это многочлен первой степени, который мы уже умеем вычислять в равноотстоящих точках. Назовем его y1(x) и составим для него разностное уравнение:
    y1(x+dx)-y1(x) = 2a2dx.
    В конечном итоге получаем систему из двух разностных уравнений:
    y0(x+dx)-y0(x) = y1(x) + a2dx^2 + a1dx
    y1(x+dx)-y1(x) = 2a2dx
    Эту систему можно легко решать в прямом времени. В правой части каждого уравнения у нас имеются только константы, которые надо прибавлять, и значения других переменных, которые мы получаем в ходе решения системы.

    Аналогичным образом для многочлена третьей степени можно составить систему из трех линейных разностных уравнений, для четвертой степени — систему из четырех уравнений и т.д. Один шаг решения каждого из уравнений — добавление константы, поэтому решение всей системы (вычисление всего полинома в следующей точке) нам обойдется во столько сложений, какова степень этого полинома.

    Если все константы и начальные условия выражаются целыми числами — то решение будет точным.

    Именно таким образом был получен лучший алгоритм для вычисления таблицы синусов на ZX Spectrum на форуме zx-pk.ru. Синус был приближен полиномом 3-й степени, который вычислялся по приведенному выше алгоритму. Размер программы и точность оказались рекордными среди предложенных вариантов.


    1. Uranix
      13.03.2016 02:32
      +1

      Этот принцип был заложен еще в разностной машине Бэббиджа


      1. MichaelBorisov
        13.03.2016 14:08

        Красота. Спасибо за наводку!


  1. avsmal
    16.03.2016 16:11

    Недавно разбирал статьи про вычисление полиномов в нескольких точках, правда меня интересовала асимптотическая сложность, а не точно количество умножений.

    В статье Charles M. Fiduccia. 1972. Polynomial evaluation via the division algorithm the fast Fourier transform revisited. In Proceedings of the fourth annual ACM symposium on Theory of computing (STOC '72). ACM, New York, NY, USA, 88-93 рассказывается как вычислить полином степени не более N в N точках за время O(N log^3 N).

    В статье самой разобраться сложно, советую смотреть Ellis Horowitz. 1974. A unified view of the complexity of evaluation and interpolation. Acta Inf. 3, 2 (June 1974), 123-133 или более поздние источники.