«Решить парадокс Грея с дельфинами может любой, а ты попробуй сделать это без дельфинов. »

Вообще то планировал я провести выходные несколько по иному, съездить на Copter Huck (не то, чтобы я был фанатом коптеров, просто посмотреть, что молодежь придумывает, потусоваться типа), но старшая сестра была категорически против. Я, конечно, настаивал (то есть пару раз хмыкнул и сказал" Ну может, все-таки… будет прикольно"), но она была неумолима, а когда супруга приняла ее сторону, шансов на поездку не осталось. Ну и ладно, «не очень то и хотелось», зато немного посидел над забавной задачкой из области программирования, которую сам себе придумал, о чем и докладываю.

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

В одном недавнем посте автор рассматривал задачу ускорения (кроме всего прочего) расчета кривых Безье (КБ) на МК со сравнительно слабыми параметрами. Ну на самом деле эти параметры на уровне среднего мейнфрейма 70х годов, но по нынешним временам считаются явно недостаточными. В результате определенных действий автору удалось несколько ускорить вычисления, на мой взгляд, явно недостаточно, вот и решил написать, как это следует делать в первом приближении. Я прекрасно знаю универсальный рецепт решения проблем с быстродействием — взять МК с частотой повыше или перейти на другое семейство, но я родом из тех времен, когда мы учились обходится тем, что есть, просто потому, что ничего другого не было, от слова совсем. По нынешним временам подход устаревший, но мне показалось, что будет не безынтересен и современным читателям Хабра.

Формулируем задачу — мы хотим как можно быстрее вычислить координаты точек на кривой Безье, задаваемой крайними точками А и Б и мнимым фокусом С. Формула для вычисления точки Р на кривой задается выражением

$$

, где Т меняется от 0 до 1 включительно. (В Вики пишут, что эта формула в свое время была секретной, странно как то, но все возможно). Понятно, что мы не будем считать в комплексном виде, вместо этого будем отдельно искать координаты X и Y. Оценим сложность вычисления по этой формуле, просто подсчитав количество знаков арифметических операций в данном выражении — 7 умножений и 5 сложений (=>7*5+). Не исключено, что хороший компилятор (а сейчас все компиляторы хорошие и отлично оптимизируют, если Вы им явно это не запретите) уменьшит затраты до 7*3+, хотя лучше бы ему помочь, вычислив (1-Т) заранее. Вообще говоря, хороший компилятор может вообще сотворить чудеса, если все значения в формуле представлены константами, но мы предполагаем все величины статически неопределенными.

Часть первая, математическая

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

$$

=>5*5+, что явно лучше исходного значения 7*5+, а вот относительно улучшенного 7*3+ следует еще подумать.

Если мы примем время выполнения операции сложения за единицу, то время выполнения умножения будет точно не меньше единицы, как правило, больше, но вот насколько — зависит от реализации МК (сначала написал — от архитектуры, но это не совсем верно). Когда на кристалле отсутствует аппаратный умножитель, время выполнения умножения будет в десятки (30+) раз превосходить единицу, когда же он присутствует, то в разы (1-6). Поэтому можно с уверенностью полагать, что замена умножения на сложение дает выигрыш (и часто значительный) по времени исполнения практически всегда. Ну и сразу заметим, что переход от чисел с фиксированной запятой к плавающей точке (оставим в стороне доказательство этого факта) приводит к увеличению времени исполнения в 20+ раз для сложения (здесь очень сильно влияет выравнивание), но всего лишь к незначительному увеличению для умножения. Поэтому для чисел с плавающей точкой время сложения и умножения отличаются мало, особенно в относительном выражении (можно ожидать раза в 2 максимум), но все равно отличаются и не в пользу умножения, так что выигрыш и тут присутствует.

Возвращаясь к предыдущему абзацу, получаем, что для ПТ оценка 5*5+ не должна иметь существенного преимущества перед 7*3+, но у нас еще остались резервы. Обратим внимание на то, что мы должны вычислить множество значения точек на кривой Безье, при изменении параметра Т, а все остальные параметры кривой зафиксированы (но не константны, а жаль), тогда остальные члены формулы можно рассчитать заранее и получаем

$$

=> 3*2+, где $$ и $$ уже неплохо, а если вспомнить схему Горнера и записать

$$

=>2*2+, то по сравнению с решением «в лоб» мы должны выиграть более, чем в 2 раза, почти в 3, причем эти оптимизации совершенно очевидны.

Проверим теорию практикой (хотя это совершенно излишне, мы уверены в своих оценках, но вдруг я недооценил компилятор), для чего надо замерить реальное время исполнения разных вариантов на настоящем железе. Ну так уж получилось, что у меня дома есть много всяких отладочных плат для МК от различных фирм (в том числе раритеты вроде отладок от Luminary Micro или Intel Edisson, попробуйте-ка нынче купить такое), но нет ни одной платы Ардуино («Ну нет у нас ананасов»). Казалось бы, тупик, но есть варианты — нам на помощь приходит весьма интересный сайт tinkercad.com, на котором Вы можете построить свою схему на макетной плате с использованием Ардуино модуля, написать скетч и тут же его выполнить. При этом Вы можете задавать точки останова, выполнять программу по шагам и даже (невиданная вещь для настоящей Ардуино) просматривать значения переменных в моменты останова.

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

8+8=>8 — 1 такт, 16+16=>16 — 2,
8*8=>16 — 2, 16*16=>16 — 14 (единственное, что оказалось неожиданным, я думал получить 4*2+4*2=16, там интересные оптимизации),
8/8=>8 — 230, 16/16=>16 — 230.

Обратите внимание на последние две цифры, из них понятно, что операция деления относится к запрещенным, если мы действительно хотим считать быстро. Теперь (наконец то) измеряем время выполнения операций над числами ПТ с мантиссой 24 разряда
а+б — 126 (и сильно зависит от операндов), а*б — 140, а/б — 482.
Полученные данные хорошо коррелируют с нашими теоретическими предположениями, видно, что аппаратная реализация на борту данного МК: для умножения есть, для деления нет, для операций ПТ нет.

Теперь начинаем измерять время полного расчета. Задаемся значениями А=140, Б=120, С=70 и строим 170 равномерно распределенных точек на КБ. Почему именно этими значениями — они были даны в указанном посте при оценке производительности. Ниже даны алгоритмы и соответствующее время исполнения теста.

Формула (1) => 20мсек или 1900 тактов на один отсчет
Формула (1) => 18мсек или 1660 тактов на один отсчет (отдельно считаем 1-T)
Формула (2) => 16мсек или 1540 тактов на один отсчет
Формула (3) => 10мсек или 923 такта на один отсчет
Формула (4) => 8мсек или 762 такта на один отсчет

Видно, что полученное сокращение времени исполнения (с 20мсек до 8мсек) хорошо коррелирует с ожидаемым и нам удалось ускорить вычисления более, чем в 2 раза. Обратим внимание, что кроме совершенно очевидных соображений и математики, не выходящей за курс средней школы, нам не потребовалось.

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

Часть вторая, вычислительная

Раз модификации алгоритма исчерпаны, остаются только структуры данных и мы вступаем на почву чисел с фиксированной точкой. Здесь нас ждет множество подводных камней, о которых мы не задумывались для ПТ — диапазон и точность (вообще то и для ПТ над этими вопросами следует задумываться, но здесь все попроще, многое сделано за нас). Надо провести небольшое исследование задачи для определения необходимого представления ФТ (выбраный в упомянутом посте 9.7, судя по результатам, явно недостаточен) но я предлагаю пойти несколько иным путем. Кстати, если принять не 170 шагов на интервале, а 128 (не вижу причины, запрещающей нам этот шаг), то данное представление нас вполне бы устроило. Если мы примем во внимание то факт, что константы, задающие КБ, заданы целыми числами, а единственный параметр Т может быть представлен дробью вида и/И и результат мы будем использовать для отрисовки на экране, то есть переводить в целые координаты, то мы можем просто все делать в целых числах, которые обрабатываются намного быстрее.

Используем только последнюю формулу и перепишем ее в новых обозначениях

$$

(=>2*2+2/), где А1 и Б1 вычисляем так же, как и для ПТ. Очевидно, что все числа целые и соответствующие операции должны выполнятся намного быстрее. Для того, чтобы не потерять точность при операции деления целых чисел (2/3 = 1 != 1.5) и проводить деление в самый последний момент, немного преобразуем формулу к виду

$$

(=>4*2+2/). Все числа ФТ, так что реализуем данный алгоритм и получаем… вот тебе, бабушка, и Юрьев день… 1869 тактов, да ведь это намного хуже, чем для ПТ, мы с этого начинали, фигня какая то, ведь целые числа намного быстрее.

Начинаем разбор полетов и оказывается, что недостаточно просто поменять тип переменных. Во-первых, мы должны использовать числа не 8 и даже не 16, а 32х разрядные, иначе наступит переполнение, а длинные числа, хоть и быстрее ПТ, но не настолько, чтобы компенсировать огрехи алгоритма.Во-вторых, вот эти огрехи — у нас опять появились константы, вычисляемые на каждом такте — убираем их предварительным вычислением Б2=Б1*И, А2=А*И*И. Тогда получим

$$

(=>2*2+2/) с результатом 1684 — лучше, чем предыдущее, но все равно не за этим мы уходили от ПТ.

Исключаем вычисление еще одной константы И2=И*И и получаем

$$

(=>2*2+1/), со временем исполнения 956 тактов — а вот это интересно, исключение одной операции привело к значительному росту производительности.

Вот что нас тормозит — деление, поскольку это очень затратная по времени операция, но для борьбы с ним у нас есть интересный трюк. Для вычисления выражения 1/И мы можем провести элементарные преобразования 1/И=1/И*(Н/Н)=1*(Н/И)/Н. Если мы в качестве Н выберем степень двойки, то деление на Н можно заменить сдвигами, а если показатель степени будет кратен 8, то даже сдвиги не понадобятся. А значение Н/И придется вычислять по честному, но только один раз, после чего в цикле вычисления остается только умножение.

Обратим внимание на то, что мы сделали не вполне корректное преобразование и заменили Н/И его округленным значением К, чтобы перейти к операциям исключительно с целыми числами. Некорректность заключается в потере точности и следует провести дополнительное исследование, чтобы доказать применимость данного подхода к нашему случаю. Запишем Н/И в виде (К*И+д)/И=К+(д/И), где д меньше И. Тогда абсолютная погрешность при переходе к Н/И к К составит д/И, а относительная — д/И/(К+д/И) >=д/И/(К+1) ~ д/И/К, при условии К>>1 (это не сдвиг). Отсюда следует, что значение Н следует выбирать как можно большим, поскольку абсолютная погрешность вычисления равна А*д/И/К >= А*1/Н/И. Если мы хотим, чтобы погрешность составила не более единицы, мы должны выдержать условие А/К<=1, тогда К>=А, преобразуем К*И>=А*И, что означает Н>=А*И, тогда мы не теряем в точности. Для нашего случая А<=256 и И<=256 получаем Н>=2**16, что вполне приемлемо. Очевидно, что в вышеприведенных формулах следует использовать модули исходных чисел.

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

В любом случае, требуемую точность мы сможем обеспечить и получаем следующий алгоритм: Н=2**16; К=[Н/И] (И<256); 0 <= и<=И;

$$

(=>4*2+2>>16) где все операции проводятся над длинными целыми числами. Реализуем данный алгоритм и получаем 583 такта… а вот это уже близко к идеалу, но пока еще не идеал.

Дальше идут мелкие настройки на конкретный МК — работа с глобальными переменными осуществляется быстрее. чем с локальными, но еще быстрее с регистровыми локальными, что приводит к снижению времени до 506 тактов.

Далее заметим, что последнее умножение перед сдвигом можно провести с 16 разрядными числами, что даст 504 — мелочь, а приятно.

Итого мы ускорили расчеты по сравнению с реализацией «в лоб» в 1900/504 — более, чем в 3 раза, причем не потеряли в точности от слова совсем. Вот такой результат я называю оптимизацией по времени, а не полученные в исходном посте 20%.

Можно ли добиться еще лучших показателей — можно, но это тема следующего поста.