К вопросу о кривых Безье и быстродействии Ардуино, часть вторая

Мы пойдем мимо — и дальше

В своем предыдущем посте я показал, как можно улучшить быстродействие расчета точек на кривой Безье (КБ) путем:

1. Преобразования расчетных формул — ускорение в ~3 раза.
2. Перехода от чисел ПТ к ФТ — ускорения почти нет, но позволяет провести 3.
3. Заменой операции деления умножением и сдвигом — ускорение еще на 40%.

В конце упомянутого поста я заявил, что можно вычисления сделать еще быстрее, ну что же «пацан сказал — пацан сделал».

Напомню еще раз полученную формулу для вычисления точки на КБ

.Очередное увеличение скорости связано с особенностью задачи — нам нужно не И раз рассчитать значение КБ при различном значении параметра «и», а найти ряд значений при изменении (в данном случае увеличении) данного параметра на известную, более того, фиксированную величину (в нашем случае единицу), что позволяет воспользоваться описанным далее приемом. В мое время это называлось разностным методом вычисления (если память мне не изменяет, по крайней мере, так его называли на лекциях), в Вашем распоряжении весь Инет, может быть (даже наверняка), есть и другое название.

Рассмотрим случай Р=А*и (=>1*), и предположим, что мы знаем значение Р0 при некотором аргументе и0. Тогда значение при следующем аргументе и0+1 можно рассчитать, как Р=А*(и0+1)=А*и0+А=Р0+А (=>1+). Нисколько не потеряв в точности, мы смогли заменить операцию умножения на операцию сложения, которая намного быстрее.

Теперь более сложный пример Р=А*и*и (=>2*), рассматриваем по аналогии Р=А*(и+1)*(и+1)=А*(и*и+2*и+1)=А*и*и+2*А*и+А=Р0+2*А*и+А(=>2*2+). На первый взгляд, мы ничего не выиграли, но, если заранее вычислить А'=2*А, то получим (=>1*2+), вполне возможен и выигрыш. Но мы не остановимся на достигнутом и к полученному выражению А'*и применим уже известную нам методику, тогда получим две операции над двумя переменными: Р=Р+А'+А; А'=А'+А (=>3+). Если мы еще примем во внимание, что начальное значение А' можно сделать больше на А, то вообще остается только два сложения вместо двух умножений. Да, нам пришлось завести две дополнительные переменные, но это класический размен — расплачиваемся памятью за время.

Отается только вычислить правильные начальные значения, но это делается только один раз в начале итераций, а если начальное значение параметра и=0, то и вообще тривиально Р=0, А'=А. Проверим наш метод на нескольких итерация (это совершенно излишне, правильно примененная математика не может дать неправильного ответа), но позволит лучше понять происходящее. Итак
итерация 0: и=0; Р=0 (верно); А'=А; А''=2*А;
итерация 1: и=1; Р=0+А'=0+А=А (верно); А'=А'+А''=А+2*А=3*А;
итерация 2: и=2; Р=А+3*А=4*А (верно); А'=3*А+2*А=5*А;
итерация 3: и=3; Р=9*А (верно); А'=7*А и так далее.

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

Перепишем наше исходное выражение для КБ в развернутой форме

, тогда для вычисления с применением данного метода нам потребуется 2+ для первого слагаемого (и две переменные), 1+ для второго (и одна переменная) и 2+, чтобы сложить все вместе (=>5+). Можно ожидать, что даже такой (неправильный) подход даст выигрыш по сравнению с 2*2+, но все намного лучше. Очевидно, что операция сложения является аддитивной (спасибо, капитан), поэтому можно сгруппировать слагаемые и заменить константные члены на новые выражения, что дает следующий алгоритм:

1. начальные значения Р=С; А'=А1+Б1; А''=2*А1;
2. на очередном шаге Р=Р+А'; А'=А'+А'' (=>2+), что несомненно быстрее схемы Горнера.

Реализуем наш алгоритм в виде программы (это не столь тривиально, как может показаться, поскольку я для упрощения забыл о необходимость сдвигов, но и не слишком сложно… минут на 20), получаем вычислительную сложность (=>2+1>>) и замеряем быстродействие — получилось 140 (увеличение быстродействия еще в 3 раза) тактов на цикл, а вот это уже почти идеальный результат. Что нам осталось сделать, чтобы получить идеальный (в некотором смысле) вариант — обратить внимание на размерность операндов в формулах. Сейчас мы проводим все операции над длинными (32х разрядными) целыми, а это в некоторых местах может быть излишне. Если сократить разрядность операндов до минимально необходимой, то можем получить выигрыш процентов на 20-25, но это потребует от нас перехода к ассемблеру (язык С не располагает подходящими для таких операций средствами) и внимательного анализа данных исходной задачи. Делать ли это — решать читателю, мы и так уже ускорили вычисления более, чем в 1900/140=13 раз по сравнению с «наивным» подходом.

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

Однако совершенно неожиданно возникла одна проблема — мы ясно видим две операции сложения над 32 битовыми числами, которые должны занять 4+4=8 тактов, положим еще 8*3*2=48 тактов на пересылки операндов, 4 такта на получение результата сдвига, 4 такта на вызов процедуры вычисления и возврат и 4 такта на организацию цикла — откуда еще 60 тактов, непонятно. Раньше мы этого просто не замечали на фоне сотен тактов вычисления, но теперь то можно и обратить внимание. Избыточные такты легко нашлись — в цикле оставались ненужные отладочные операции, если все аккуратно подчистить, то остается только 120 тактов и назначение каждого вполне понятно (ну не так, что бы уж совсем понятно, но все таки). Далее выясним время выполнения пустого цикла — 22 такта, интересно, чем они там занимаются все это время, ну да ладно, и очистим собственно время вычисления, которое составит 98 тактов. Обратим внимание, что изменяется оценка полученного ускорения работы: вместо 1900/140=13 получаем (1900-50)/(140-40)=19 раз, что не имеет никакого практического смысла, поскольку цикл все равно необходим, но позволяет еще больше поднять самооценку.

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

Ну а в заключение о примечании «в некотором смысле» — если речь идет о последовательном вычислении координат очередной точки на КБ при изменении параметра (олицетворяющего собой время), то предложенный алгоритм (после всех оптимизаций) улучшить уже не удается. Но если задачу переформулировать и, к примеру, задаться целью просто построить КБ на экране (без привязки ко времени), то тут есть весьма перспективный метод, ключевое слово «Брезенхайм», но «это уже совсем другая история», которой я посвящу отдельный пост (может быть, когда-нибудь, если старшая сестра не будет против).