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

▍ Введение

Существуют три основных подхода к аппроксимации функций многочленами:

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

2. Метод наименьших квадратов:
Чуть менее известный, более сложный, но и более эффективный. Его идея заключается в минимизации площади между функцией и её аппроксимирующим многочленом. Недостатком является тот факт, что из минимизации площади вовсе не следует минимизация максимальных отклонений — поэтому ошибка аппроксимации хоть и осциллирует, но всё равно имеет тенденцию возрастать на краях диапазона.

3. Многочлены Чебышева:
Ещё более сложный, наименее известный, но часто наиболее оптимальный метод. Оптимальность его проистекает из того, что ошибка распределена по всему диапазону и ограничение максимального отклонения задаётся by design. Сложность вытекает из того, что сначала нужно считать коэффициенты для каждого отдельного многочлена Чебышева через далеко не простые интегралы даже для простых функций (в частности, для синуса получаются функции Бесселя, которые далеко не в каждом калькуляторе присутствуют), а затем их пересчитывать для суммарного итогового многочлена с запасом по точности. К счастью, имея в своём распоряжении современные системы компьютерной алгебры, всё это можно сделать лёгким движением руки, даже не вникая особо глубоко в математические основы происходящего.

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

▍ Подготовка

Используя Wolfram Mathematica, определим функцию, возвращающую аппроксимирующий многочлен заданного порядка для произвольной функции:


Здесь точность полученных коэффициентов равна 20 десятичным знакам — для наших целей этого достаточно. Можно также задать отдельные реализации для чётных/нечётных функций, чтобы ускорить вычисления. Для нечётных функций она будет выглядеть как


а для чётных как:

▍ Начало

Наиболее надёжный способ ускорить вычисление синуса — просто уменьшить степень аппроксимирующего многочлена, пожертвовав точностью. При использовании многочленов Чебышева количество верных знаков ориентировочно соответствует степени многочлена . Для примера возьмём 7-ю степень — её вполне достаточно для задач визуализации:


Получили многочлен:

Реально полученную погрешность можно увидеть, если вывести на график разницу между исходной функцией и этим многочленом:



По графику видно, что на краях диапазона погрешность принимает максимальное значение. Это значит, что sin(90°) с таким многочленом будет равен не единице, а 0.999999414195269. Если важно, чтобы была именно единица — то для этого можно пойти на хитрость и немного расширить диапазон аппроксимации таким образом, чтобы отклонение на краях стало близким к нулю. Значение этого «немного» (как множитель при аргументе синуса) будем находить итеративно через численный поиск корня ошибки (разницы между функцией и аппроксимирующим многочленом), наиболее близким к единице.


Результат:

Вот как теперь выглядит график отклонений:


Погрешность неизбежно выросла, но совершенно незначительно.

Чтобы получить полноценную функцию синуса, необходимы вспомогательные манипуляции для приведения аргумента к диапазону , учёта цикличности и симметрии. Для экспорта полученного многочлена в код нам также пригодятся вспомогательные (уже имеющиеся) функции: HornerForm — для преобразования многочлена к схеме Горнера, и CForm — для преобразования его к виду, традиционно используемого в си-подобных языках программирования.


Сами функции будут выглядеть так:



В версии для C/C++ тип bool заменён на int — как показали измерения, компилятор в этом случае генерирует более эффективный код.

Здесь нужно сделать важное замечание. В целях оптимизации для получения целой части аргумента здесь используется не функция Floor(x), а приведение к целому типу (чтобы избежать вызова внешней библиотечной функции). Это приводит к тому, что аргумент не должен превышать максимальное для int значение, а именно 2147483647. Тип int можно заменить на long (а правильнее — на unsigned long), но в зависимости от платформы и компилятора это также может привести к замедлению кода. На достаточно больших значениях аргумента также возрастает погрешность из-за разности порядков у аргумента и периода . Подобные ограничения являются неизбежной расплатой за более быстрый код.

▍ Продолжение

При увеличении степени многочлена увеличивается количество вычислений и как следствие, накапливается дополнительная погрешность. Поэтому для близкой к double точности разумно разбить диапазон на два и для каждого рассчитывать отдельный многочлен. Для левой части будем считать многочлен для , а для правой для (косинус «переворачивается» для того, чтобы выкинуть постоянную константу, смещающую ошибку к центру, обеспечив таким образом равенство ). После получения таким образом многочленов нужно вернуть аргумент к исходному диапазону, пересчитав коэффициенты заменой . Опытным путём было установлено, что для достижения необходимой точности достаточно 13-ой и 12-ой степени (для синуса и косинуса соответственно):


Графики ошибок:


Здесь ошибка для косинуса получилась чуть больше, чем у синуса — но это не имеет значения, потому что они всё равно замаскируются накопленными погрешностями при вычислении многочлена.

▍ Табличный метод

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

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

Оптимальный размер таблицы n выбирается экспериментально, чтобы разница между функциями , и их аппроксимирующими многочленами выбранной степени не превышала заданной погрешности. Для 5-ой степени достаточно таблицы из 128-и значений, для 3-й потребуется уже 1024.


Графики ошибок для 5-ой (синус) и 6-ой (косинус) степеней:


Код для получения таблицы тоже предельно прост:

With[{n = 128}, Table[N[Sin[x Pi/2], 18], {x, 0, 1, 1/n}]]


Значения синуса и косинуса из неё получаются в зависимости от того, индексировать её с начала в конец или наоборот. Стоит уточнить, что здесь по факту размер получается на единицу больше, то есть 129 (из-за замыкающей единицы для ). Сам же размер как степень двойки выбран для того, чтобы избежать потери точности при умножении на него аргумента и обеспечить табличную точность в значениях (что может быть полезным в алгоритмах быстрого преобразования Фурье). В случае если аргумент синуса удобнее задавать в других единицах — скажем, в градусах — то размер таблицы разумно делать кратным 90, чтобы получать максимально точные значения при целых значениях аргумента.



Табличный метод можно использовать и по-другому — только для синуса, если разбить интервал на несколько и хранить в ней коэффициенты аппроксимирующих многочленов небольшой степени для каждого из них по отдельности. Но таким образом не получится сэкономить за счёт чётности/нечётности функций — и для многочлена скажем 5-ой степени придётся хранить и вычислять все 6 коэффициентов — что не даст прироста в производительности, но зато может дать несколько дополнительных бит точности.

▍ Как тестировать

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

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

Для косинуса получается похожий результат:

Также нужно учитывать, что некоторые современные оптимизирующие компиляторы (в частности C++ от Intel и Microsoft) могут распознать сумму и вычислять её в 4 параллельных потоках на регистрах AVX, что при замерах производительности может показать значительное превосходство библиотечных реализаций (да, я столкнулся с этим). Чтобы этого избежать, можно на каждом этапе менять знак (то есть считать функцию ), аналитическим результатом которой будет ноль (для ). Но такой подход подойдёт только для замера производительности и оценить погрешность по отклонению результата от нуля уже не получится — поскольку будет происходить взаимная компенсация погрешностей.

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



В варианте для C/C++ добавлено немного бесполезных вычислений во избежание автоматической векторизации. Также по-хорошему его можно реализовать на шаблонах, чтобы не дублировать одинаковые куски кода для разных функций.

В обоих случаях время измеряется в секундах с точностью до миллисекунд. Сделано это для того, чтобы привести результаты к общему знаменателю и не зависеть от платформо-зависимых средств и отсутствующем в си std::chrono. Ну и заодно стимулировать выполнять код достаточно долго, чтобы накопилось хотя бы несколько секунд.

▍ Заключение и результаты

При замерах производительности все рассмотренные реализации оказались быстрее библиотечных и на C#, и на C++ с компилятором от Microsoft. На C++ была дополнительно протестирована функция fsincos из математического сопроцессора FPU процессоров x86/x64, которая тоже умеет считать и синус, и косинус одновременно.

Сводная таблица по результатам замеров на Core i7:
Последний столбец показывает отношение между одними и теми же функциями — как видно, разница не такая уж и значительная (не считая библиотечных реализаций). Также можно заметить, что С++ более эффективно работает с таблицами (что тоже вполне ожидаемо). А вот ассемблерная fsincos показала наихудший результат — даже несмотря на её аппаратную реализацию.

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

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

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