Схема Горнера, вещественная версия, последовательный вариант

Основные авторы описания: А.В.Фролов, Вад.В.Воеводин (раздел 2.2)

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

1.1.1 Решаемая задача

Схема Горнера решает задачу деления многочлена [math]P_n(x)[/math] с известными коэффициентами на двучлен [math]x - \alpha[/math]. В качестве результатов выступают коэффициенты многочлена [math]Q_{n - 1}(x)[/math] из соотношения

[math]P_n(x) - P_n(\alpha) = (x - \alpha) Q_{n - 1}(x)[/math],

а также [math]P_n(\alpha)[/math] (значение многочлена [math]P_n(x)[/math] в точке [math]\alpha[/math])

К сожалению, зачастую в учебной литературе схему Горнера сводят к вычислению значения многочлена [math]P_n(x)[/math] в точке [math]\alpha[/math].

1.1.2 Общая схема

Фактически схема Горнера реализует «деление уголком» многочлена на двучлен [math]x - \alpha[/math], с нахождением остатка, который, по теореме Безу, равен значению многочлена [math]P_n(x)[/math] в точке [math]\alpha[/math].

1.2 Математическое описание алгоритма

Исходные данные: одномерный массив [math]n + 1[/math] чисел [math]a_k[/math] и скаляр [math]\alpha[/math].

Вычисляемые данные: одномерный массив [math]n + 1[/math] чисел [math]b_k[/math].

Формулы метода с выводом: из соотношения

[math]P_n(x) - P_n(\alpha) = (x - \alpha) Q_{n - 1}(x)[/math]

если записать оба многочлена в каноническом виде

[math] \begin{align} P_n(x) & = a_0 x^n+ a_1 x^{n - 1} + \dots + a_{n - 1} x + a_n, \\ Q_{n - 1}(x) & = b_0 x^{n - 1} + b_1 x^{n - 2} + \dots + b_{n - 2} x + b_{n - 1} \end{align} [/math]

и приписать «свободному» [math]b_n[/math] значение [math]P_n(\alpha)[/math], то, при подстановке многочленов в каноническом виде в исходное соотношение и приравнивании коэффициентов равных степеней, в качестве решения получается

[math] \begin{align} b_0 & = a_0, \\ b_k & = a_k + \alpha b_{k - 1}, \quad k = 1, \dots, n \end{align} [/math]

что и является формулами схемы Горнера.

1.3 Вычислительное ядро алгоритма

Вычислительное ядро схемы Горнера в последовательном варианте можно представить в качестве последовательного набора [math]n[/math] «двойных» операций (умножение элементов получаемого массива на один и тот же скаляр и добавление результата к следующему элементу входного массива).

1.4 Макроструктура алгоритма

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

1.5 Схема реализации последовательного алгоритма

Формулы метода описаны выше. Последовательность исполнения — по возрастанию индекса k.

1.6 Последовательная сложность алгоритма

Для вычисления схемы Горнера для многочлена степени [math]n[/math], количество операций умножения равно количеству операций сложения и равно [math]n[/math]. Поэтому алгоритм должен быть отнесён к алгоритмам линейной сложности по количеству последовательных операций.

1.7 Информационный граф

Опишем граф алгоритма в виде рисунка. Изображён граф для [math]n = 9[/math] (рис. 1). Как видно, граф чисто последовательный. Ввод и вывод обозначен только для начального коэффициента и значения многочлена. Op - операция "умножить входное данное на скаляр и прибавить к другому данному".

Рис. 1. Схема Горнера: последовательная версия

1.8 Ресурс параллелизма алгоритма

Последовательный вариант вычисления схемы Горнера не имеет ресурсов параллелизма. Его ярусно-параллельная форма (ЯПФ) — единственна и совпадает с последовательным алгоритмом. Таким образом, в получившемся алгоритме высота параллельной формы будет равна [math]n[/math] операций умножения плюс [math]n[/math] операций сложения. В таком виде алгоритм должен быть отнесён к алгоритмам линейной сложности по высоте параллельной формы. Ширина параллельной формы равна [math]1[/math], что даёт нам постоянную сложность по ширине параллельной формы.

1.9 Входные и выходные данные алгоритма

Входные данные: массив [math]a[/math] (элементы [math]a_i[/math] с номерами от [math]0[/math] до [math]n[/math]), скаляр [math]\alpha[/math].

Дополнительные ограничения: отсутствуют.

Объём входных данных: [math]n + 2[/math].

Выходные данные: массив [math]b[/math] (элементы [math]b_k[/math] с номерами от [math]0[/math] до [math]n[/math]).

Объём выходных данных: [math]n + 1[/math].

1.10 Свойства алгоритма

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является константой (1). При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных — всего-навсего 1 (входных и выходных данных даже на 1 больше, чем количество операций). При этом алгоритм полностью детерминирован. Дуги информационного графа локальны. По устойчивости схема Горнера оптимальна для вычисления значения многочлена с известными коэффициентами.

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

На Фортране схему Горнера можно записать так:

	b(0) = a(0)
	DO I = 1, N
		b(I) = a(I)+b(I-1)*alpha
	END DO

2.2 Локальность данных и вычислений

2.2.1 Локальность реализации алгоритма

2.2.1.1 Структура обращений в память и качественная оценка локальности

Рис. 2. Реализация схемы Горнера. Общий профиль обращений в память

На рис. 2 представлен профиль обращений в память последовательной вещественной реализации схемы Горнера. Исходя из рисунка 2, данный профиль устроен очень просто и состоит из двух последовательных переборов элементов, выполняемых параллельно для двух массивов.

Однако даже в таком простом примере для полного понимания структуры обращений в память нужен более детальный анализ на уровне отдельных обращений. На рис. 3 рассмотрим подробнее фрагмент 1, выделенный на рис. 2 зеленым. Можно увидеть, что верхний последовательный перебор несколько отличается от нижнего – в первом случае к каждому элементу выполняется по два обращения подряд. Отметим, однако, что данное уточнение строения профиля слабо сказывается на локальности всего профиля.

В целом, общий профиль обладает высокой пространственной локальностью, поскольку перебор элементов массивов осуществляется последовательно, однако очень низкой временно́й локальностью – в одном из массивов к каждому элементу выполняется обращение только дважды, во втором массиве повторные обращения вообще отсутствуют.

Рис. 3. Фрагмент 1 (первые 100 обращений общего профиля)

2.2.1.2 Количественная оценка локальности

Основной фрагмент реализации, на основе которого были получены количественные оценки, приведен здесь (функция KernelHorner). Условия запуска описаны здесь.

Первая оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.

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

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

for (int i = 1; i < size; i++) {
    c[i] = a[i] + c[i - 1] * x;
}

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

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

Рис. 4. Сравнение значений оценки daps

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

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

Рис. 5. Сравнение значений оценки cvg

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

2.3 Возможные способы и особенности параллельной реализации алгоритма

Описываемый алгоритм не предполагает параллельной реализации.

2.4 Масштабируемость алгоритма и его реализации

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

2.7 Существующие реализации алгоритма

Помимо выписанной выше простейшей реализации, существуют более примитивные коды, реализующие часть функционала схемы Горнера. Это объясняется тем, что для многих нужно только значение многочлена (остаток от деления), но не нужен многочлен-частное. В таком случае вместо всех элементов массива [math]b[/math] используется один и тот же скаляр.