Как посчитать синус быстро

и точно. Точнее, с заданной точностью, простите за каламбур.

Под катом я расскажу, как сделать это с использованием школьного курса алгебры и целочисленной арифметики, при чём здесь полиномы Чебышёва I-го рода, и дам ссылки на примеры реализаций для ПК и Cortex-M3.

Если вам лень вникать в теорию, и нужны только примеры реализации - сразу переходите в конец статьи.

Поехали.
Если нужно найти синус, или другую тригонометрическую функцию на ПК, это делается просто - сейчас в процессоры семейства x86/x64 встроен FPU, а он знает инструкцию fsin, при помощи которой вычисляется искомое. Вычисляется, если говорить с натяжкой, быстро.

Если это надо сделать это на МК без плавающей точки - то возникают проблемы. Можно использовать функцию из поставляемой вместе с компилятором библиотеки, будет точно, но очень медленно. Если надо быстрее - то первое, что приходит в голову, заранее посчитать таблицу со значениями, но точность при этом сильно упадёт, и будет зависеть от шага аргумента между смежными значениями. Следующий интуитивно понятный шаг - использовать интерполяцию по двум соседним точкам таблицы. Это поможет поднять точность, но несильно. Иногда для достижения нужной точности размер таблицы всё равно превосходит разумные пределы.

Что же тогда делать? Увеличивать степень аппроксимации. Это позволит увеличить точность вычислений и (или) уменьшить размер таблиц. И сделать это совсем несложно.

Аппроксимация полиномами

Кусочно-линейная аппроксимация - это, по сути, аппроксимация полиномом первой степени. Выходное значение вычисляется по формуле y = A1·x + A0, где x - это входной аргумент, A0, A1 - некие коэффициенты.

Для примера разобьём весь период от 0 до 2π на 8 равных интервалов, и для каждого подберём такие коэффициенты A0 и A1, что бы максимум ошибки при любом аргументе был минимален (как это сделать, я расскажу ниже). Если мы построим график этого "синуса", он будет выглядеть примерно так (сплошная линия):

Штриховой линией нарисовано более точное значение синуса. А график ошибки (разницы между "настоящим" и аппроксимированным значениями будет таким:

На нём, если приглядеться, не выполняется условие "максимум ошибки при любом аргументе был минимален", но сейчас это неважно. Максимальная ошибка составляет 0.03684497.

Небольшое отступление. Числа вроде 0.03684497 или 2.448728e-09 слишком абстрактны для понимания, их трудно сравнивать. Мне, как человеку, долгое время работающему в двоичной системе, гораздо ближе биты. Поэтому к ним я и буду приводить величину ошибки, беря от неё минус логарифм по основанию 2, тогда она превращаются соответственно в 4.76 или 28.6 бит. По мне, такие логарифмические "попугаи" проще для восприятия и нагляднее. По итогам предыдущего абзаца можно сказать, что "значение синуса вычислено с точностью 4.76239 бита".

Теперь давайте разобьём один период (от 0 до 2π) не на 8, а на 64 интервала, по 16 на каждый квадрант. И рассмотрим один интервал в районе, например, 45°. График ошибки (сплошная чёрная линия) выглядит так:

Он очень похож на график параболы y = x², который нарисован там же красным штрих-пунктиром. Так же, как парабола, выглядит график полинома Чебышёва I-рода 2-й степени (более того, он ей и является). На КДПВ, он изображён голубой линией.

А что если мы учтём это, и попробуем на этом же отрезке вычислять синус по формуле y = A2·x² + A1·x + A0?

Считаем тройку коэффициентов А0, А1, А2 (как - об этом ниже), по ним считаем целевой синус, и находим ошибку (сплошная чёрная линия):

Красным штрих-пунктиром нарисован полином Чебышёва3-й степени, и снова оба графика очень похожи.

Здесь я открою маленький секрет. График синуса, найденный аппроксимацией функцией 2-степени по 8 интервалам, с точностью до долей пикселя совпадает с "настоящим", и, будучи нарисованными вместе, они полностью перекрывали бы друг друга. Тремя картинками выше, где нарисованы графики, сиреневым пунктиром изображён не "настоящий" синус, а полученный функцией 2-й степени.

Как сильно увеличилась точность при переходе от 1-й ко 2-й степени? Для 64 интервалов на период она возросла с 10.7 до 17.63 бит. Откуда взялись эти числа, я расскажу ниже.

Если мы дальше продолжим увеличивать степень полинома, то ошибки будут ещё меньше, и будут приобретать формы полиномов Чебышёва разных степеней:

Точность аппроксимации на 64 отрезках при этом возрастёт до 24.980 бит (3-я степень) и 32.651 бит (4-я степень).
Если дальше повышать степень - тенденции сохранятся. Думаю, смысл применения полиномов для аппроксимации синуса понятен.

Как считать коэффициенты

В википедии, в статье про полиномы Чебышёва есть фраза: "Многочлены Чебышёва играют важную роль в теории приближений, поскольку корни многочленов Чебышёва первого рода используются в качестве узлов в интерполяции алгебраическими многочленами". Это как раз наш случай.

Посмотрим ещё раз на их графики:

Голубая линия (полином 2-й степени) пересекает ось абсцисс (y=0) в двух точках x = ±0.7071. Оранжевая (3-й степени) - в трёх точках: x = 0 и x = ±0.866. Это и есть корни полиномов, они будут использованы "в качестве узлов в интерполяции".

Ещё небольшое отступление, как представлять аргументы. На 32-битных процессорах полный период (круг) удобно представить как 2³² (0x100000000), а угол, соответственно, в диапазоне от 0 до 0xffffffff. Если количество интервалов аппроксимации, на которые разбит круг, равно 2ⁿ, то угол можно интерпретировать как двоичное число с фиксированной точкой. Например, при разбиении на 64 интервала (2⁶), угол будет интерпретироваться как число размерности 6.26. Здесь старшие 6 бит - это номер интервала (считая от нуля), а младшие 26 бит - смещение внутри него.

Возьмём, для примера, угол 15°. Если вычислить 15°/360° · 2³², то получим 0000_1010_1010_1010_1010_1010_1010_1011 в двоичном виде. В размерности 6.26 будет выглядеть как 000010.10101010101010101010101011. Здесь первые 6 бит 000010 - это десятичное 2 (2-й интервал), а оставшиеся 26 10101010101010101010101011, будучи поделены на 2²⁶, дадут 0.66666667 - это смещение внутри интервала. То же самое получим: 2.66666667 = 64 * 15/360.

Смещение внутри интервала лежит в диапазоне 0.0 ≤ x < 1.0. К нему же надо привести корни многочленов Чебышёва, которые находятся между -1 и +1. Корни полинома 2-й степени превратятся 0.14645 и 0.85355, 3-й степени - 0.066987, 0.5 и 0.933013.

Теперь давайте найдём коэффициенты A0, A1 и т.д. При аппроксимации полиномом 1-й степени y = A1·x + A0 нам нужно, что бы на заданном интервале с номером N в точке со смещением 0.14645 целевое значение было равно y=sin((N + 0.14645)/64 * 2π), а в точке 0.85355 - y=sin((N + 0.85355)/64 * 2π). Это делается при помощи системы из двух линейных уравнений с 2 неизвестными A1 и A0:

A1·0.14645 + A0 = sin((N + 0.14645)/64 * 2π)
A1·0.85355 + A0 = sin((N + 0.85355)/64 * 2π)

Для интервала N=2 (где находится 15°), получаем: A1 = 0.09521 и A0 = 0.19523.

Пройдясь по всем N, от 0 до 63, получим таблицу с наборами коэффициентов A1 и A0. С ней уже можно считать синус с точностью 10.7 бит. Как это делать, расскажу ниже (если кто до сих пор не понял сам).

Перейдём ко 2-й степени y = A2·x² + A1·x + A0. В качестве аргумента x подставим, соответственно, корни полинома 3-й степени 0.066987, 0.5 и 0.933013. Напишем систему из 3 уравнений с 3 неизвестными A2, A1 и A0:

A2·0.066987² + A1·0.066987 + A0 = sin((N + 0.066987)/64 * 2π)
A2·0.5² + A1·0.5 + A0 = sin((N + 0.5)/64 * 2π)
A2·0.933013² + A1·0.933013 + A0 = sin((N + 0.933013)/64 * 2π)

Решения для интервала N=15 будут следующие:

A2 = -0.004812613
A1 = +0.009628370
A0 = +0.995184425

Обратите внимание на A2 и A1. Если их умножить 2⁷ и 2⁶ соответственно, то их значения всё равно будет лежать в пределах от -1 до +1. Интервал №15 я выбрал не случайно - на нём значение A0 максимально и близко к 1.

Вообще, в большинстве случаев коэффициенты при больших степенях можно увеличить на некий коэффициент. При целочисленных операциях это уменьшит погрешность вычислений, а на Cortex-M3 к тому же сократит их время - об этом я расскажу ниже.

Для вычисления таблиц других размеров в предыдущую систему уравнений вместо 64 нужно подставить нужный размер. Для аппроксимации полиномом степени P нужно найти корни полинома Чебышёва степени P+1, и записать систему из P+1 уравнений с P+1 неизвестными, не забывая возводить корень многочлена в нужную степень 'n' при каждом An. (Если предыдущее предложение непонятно, то ничего страшного. Ближе к концу статьи будет ссылка на готовый генератор таблиц и краткая инструкция к нему.)

Точность вычислений

Разные комбинации размера таблицы и степени полинома дают разную погрешность. Не имея представления, как определить её аналитически, я решил считать "в лоб", перебирая все возможные комбинации, и прогоняя через каждую все углы от 0 до 0xffffffff. Для этого был создан проект accuracy_test. Суть в следующем:

• в проект включены таблицы, степенью полинома от 1 до 6 и размером от 4 до 65536, всего 90 штук. Какие таблицы проверяем, указывается в аргументах при запуске. Можно указать набор или диапазон проверяемых таблиц;

• с помощью заданной таблицы производится аппроксимация синуса по всем возможным 4294967296 аргументам. Результат аппроксимации сравнивается со значением, полученным стандартной функцией sinl(), получается разница-ошибка;

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

Количество аргументов от 0 до 2³²-1 велико, но не бесконечно, и их перебор занимает какое-то конечное время. Ноутбук с Intel Core-I7 перебирает их примерно за 6, а одноплатный компьютер на Allwinner A20 за 25 минут. Причём на ПК, как показал один эксперимент, львиную часть времени занимает вычисление "эталонного" значения функцией sinl.

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

На ней видно, что с увеличением размера таблицы точность растёт, и тем быстрее, чем выше степень полинома, при этом не поднимаясь выше 54 бит. Это самая большая точность, которую мне удалось достичь, при этом вычисления производились на ПК, с помощью 80-битных чисел с плавающей точкой (long double). Почему я заострил внимание "на ПК", будет чуть дальше.

При использовании 32-битных целочисленных вычислений результаты скромнее:

Здесь точность не превышает 30.36 бит. Но с учётом того, что старший бит отведён под знак числа, ошибка получается менее самого младшего бита.
На этой диаграмме, как видно, присутствуютне все комбинации размера таблицы и степени полинома. Это связано с тем, что некоторые коэффициенты при малом количестве интервалов превышают 1, и при умножении на 0x7f000000 выходят за пределы знакового 32-битного числа. Компилятор предупреждает об этом, а попытки определить точность дают результат в районе 0 бит.
Максимальная точность на 32-битных целых числах, которую я получил, составляет 30.37 бит при множителе 0x7fffff00.

Несколько слов о вычислениях на ARM со встроенным АРМовским FPU (armhf)

Результаты неприятно удивили.

Во-первых, выяснилось, что он может работать только с 64-битными числами с плавющей точкой (в отличии от 80-битных long double на x86/x64), и нет разницы, переменные какого типа пытаться использовать - double или long double, компилятор всё сводит к double.

Во-вторых, вычисления на ARM по точности немного отличаются от таких на x86 в худшую сторону - при долгой работе с числами набегают ошибки округления. Это можно увидеть, если посмотреть на соответствующие диаграммы для 64-битных float (файл формата ods, xls). ARM хуже вычисляет как и сами таблицы, так и аппроксимацию синуса на их основе. Причём худший вариант - когда он вычисляет на таблицах, созданных им же. Использование таблиц от x86 дают немного лучший результат, значит таблицы лучше не создавать на ARMах.

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

- для чисел с плавающей точкой - ods, xls;
- для целых 32-битных чисел - ods, xls;
- для чисел с плавающей точкой (ods, xls), когда в качестве эталона использована не sinl, а аппроксимация теми же таблицами.

О последнем надо сказать особо.

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

- float 32-бит, мантисса 23+1 бит, точность аппроксимации - 24 бита.
- float 64-бит, мантисса 52+1 бит, точность аппроксимации - 53 бита.
- float 80-бит, мантисса 63+1 бит, точность аппроксимации - почти 54 бита, хотя можно было ожидать 64.

Куда делись 10 бит?

Возникло предположение, что неточно вычисляется "эталонный" синус функцией sinl. Для проверки был проведён эксперимент, где сравнивались результаты аппроксимации двумя разными таблицами. В качестве эталона использовались несколько таблиц длиной 2¹⁵, 2¹⁶, разных степеней полинома, а сами вычисления производились 80-битными long double.

Предположение не подтвердилось, зато было замечено, что время проверки уменьшилось в несколько раз! Проверка таблицы 65536_p^1 функцией sinl заняла 279с, а аппроксимацией - 56с. Получается, что более 80% времени процессор был занят вычислением функции sinl.

В свете этого можно рекомендовать использовать на ПК аппроксимацию вместо sinl, можно получить выигрыш во времени.

А куда же делись 10 бит, я так и не понял. Есть подозрение, что точность, равная разрядности мантиссы на 32- и 64-битных float получалась за счёт запаса 80-битных вычислений, результаты которых потом огрублялись до заданных 32 и 64, а без огрубления она как раз составляла те самые 54 бита. При вычислениях на ARMах с 64-битным FPU точность вычисления упала до 50.8 бита против 53 на x86/x64.

Таблицы с коэффициентами

Их можно получить программой make_table. Там есть файл проекта на code::blocks, исходные тексты + скрипт для самостоятельной сборки с помощью gcc, а так же готовый .exe.

В качестве аргументов программа принимает:

- размер таблицы, равный 2ⁿ;
- степень полинома;
- множитель для коэффициентов для приведения их к целочисленным значениям;
- параметр AC_SHIFT, показывающий, на сколько бит множитель при коэффициенте, например, 3-й степени будет "весомее" множителя предыдущей 2-й.

Результатом работы программы является исходный код на языке С. Программа выводит его непосредственно на экран, и его можно перенаправить в файл.

Например:

.make_table.exe 64 3 0x40000000LL 3 > ..tablessine_approx_64_3_3.c

создаст файл sine_approx_64_3_3.c. Этот файл содержит таблицу коэффициентов для аппроксимации полиномами 3-й степени по 64 отрезкам, чего достаточно для получения 25-битной точности. Множитель при коэффициенте A0 = 0x40000000, а множители при коэффициентах более высоких степеней будут отличаться от соседних в 2³ = 8 раз.

В начале файла есть набор директив препроцессора, с помощью которых можно настраивать некоторые параметры таблицы:

TABLE_TYPE - тип данных таблицы, по умолчанию int32_t. Его можно "перебить", указав компилятору другой.
Q1_ONLY - указание использовать только первый квадрант, позволяет уменьшить объём таблицы в 4 раза.
DC - "добавка" к ко всем коэффициентам перед превращением их в целочисленный вид, позволяет делать арифметическое округление. По умолчанию равно 0.5 (для int32_t), для чисел с плавающей точной следует задать его равным 0.
AC0 - множитель для коэффициента 0-й степени, по умолчанию 0x40000000LL. Может быть задан другой, для таблицы типа float можно задать 1.
AC_SHIFT - на сколько бит будет сдвинут влево множитель AC0 для каждой последующей степени (механизм этого понятен при формировании множителей AC1, AC2 и так далее). Для таблицы длиной 64 строки оптимальное значение 3, для таблиц другой длины оптимальные значения лежат в районе от 1 до 5.

Каждая строка полученной таблицы выглядит примерно так:

-0.00015749713825 * AC3 + DC, -0.000000170942 * AC2 + DC, +0.098174807817 * AC1 + DC, -0.000000001187 * AC0 + DC, // 0

где
-0.00015749713825096520473455514 * AC3 - коэффициент для x³, -0.00000017094269773828251637917 * AC2 - для x², и так далее.

В конце файла, после таблицы, сформирован отчёт, какой "запас по битам" имеется у каждого коэффициента. В процессе расчёта множителей находятся таковые с максимальными абсолютными значениями. В таблице из примера для множителя 0 степени это 0.999999969786949, что почти рано 1.0, и запаса по битам нет. А максимальный коэффициент 3-го порядка равен 0.00015749713825, его можно безболезненно умножить на 2¹², и результат будет всё ещё меньше 1. С помощью запаса по битам при желании можно скорректировать сомножители AC1 ... ACn.

Если есть желание повторить проверки точности, то в папке проекта есть скрипт на lua, который может сформировать необходимые исходные файлы, и сформирует cmd-файл, запустив который получится 90 файлов с таблицами на 2² ... 2¹⁶ строк и степенями от 1 до 6. Впрочем, все они и так присутствуют

Как это использовать на практике

Для использования аппроксимации в своём проекте для начала нужно определиться, с какой точностью нужно вычислять синус. Потом нужно оценить ресурсы (объём памяти и быстродействие вычислителя), которые у вас есть.

К примеру, вы хотите, используя 32-битные целочисленные операции, сделать генератор с 24-битным ЦАП. Старший бит отведён на знак, поэтому точность аппроксимации должна быть не меньше 23 бит. На диаграмме видно, что заданная точность достигается при использовании:

- таблицы на 8192 строки при использовании аппроксимации 1-й степени,
- 512 строк при 2-й степени,
- 64 строки при 3-й степени,
- 32 при 4-й степени,
- 16 при 5-й степени, и
- 8 при 6-й степени.

При желании степень полинома можно увеличить ещё, длина таблицы при этом уменьшится.

Теперь следует определиться с ресурсами и приоритетами их использования. Если вы хотите получить максимальное быстродействие, нужно использовать аппроксимацию 1-й степени, при этом таблица на 8192 строки займёт 64к. Если вы используете какой-нибудь STM32 на Cortex-M3, то размещение таблицы в ОЗУ уменьшит время вычисления в отличии от размещения её в памяти программ, и у ST можно найти контроллеры с объёмом ОЗУ больше 64к.

Если у вас очень мало памяти и есть запас по быстродействию (FPGA например) , можно использовать аппроксимацию с высокими степенями.
Таблицы 5-й и 6-й степеней займут, соответственно, 16 * (5+1) * 4 = 384 и 8 * (6+1) * 4 = 224 байт.

Уменьшить объём ещё в 4 раза можно, используя симметричность функции синуса. Точность при этом может упасть, т.к. при аргументе функции 0x40000000, что равно 90°, будет использован 0x3fffffff, вместо0x40000001 - 0x3ffffffe, и т.д. Для использования этого фокуса, при сборке проекта в опциях компилятора нужно добавить Q1_ONLY - это уменьшит размер таблицы и добавит лишние вычисления в функции get_sine_int32() и get_sine_float().

Примеры использования

Есть 2 готовых примера - для ПК и для популярного STM32F103C8.

Что бы избавиться от лишних умножений для возведения в степень, формула
y = A2·x² + A1·x + A0
преобразована к виду
y = ((0 + A2)·x + A1)·x + A0

В упрощённом виде код выглядит так:

// angle - угол (доля полного периода в диапазоне от 0.0 до 1.0, фикс.точка 0.32)
// pN - двоичная степень количества интервалов; N = 2^pN
// ac_shift - двоичная степень отношения "весов" смежных коэффициентов
// pw - степень полинома
// table - таблица коэффициентов, одномерный массив размером 2^pN * (pw + 1)
int idx = angle >> (32 - pN); // индекс строки в таблице
uint32_t X = (angle << pN); // X = дробная часть угла * 2^32 (фикс. точка 0.32)
X >>= ac_shift; // уменьшая X, корректируем "веса" коэффициентов
const int32_t* A = table + idx * (pw + 1); // указатель на первый коэфф (при степени pw) в строке таблицы
int64_t sum = A[0];
for(int i = 1; i <= pw; i++)
{
sum *= (int32_t)X; // sum - промежуточный результат, фикс. точ. размерности 32.32
sum >>= 32; // приводим sum к размерности 64.0
sum += A[i]; // прибавляем очередной коэфф.
}
// теперь sum содержит искомое значение

ПК

Там есть 2 функции на выбор - для целочисленных вычислений get_sine_int32 и для вычисления с плавающей точкой get_sine_float. Каждая оформлена в отдельной паре файлов (.h / .c), которые нужно включить в свой проект. Так же в проект надо включить нужную таблицу с коэффициентами. Использование таблиц и функций с плавающей точкой на ПК оправдано, т.к. это позволяет ускорить вычисления.

Cortex-M3

Проект на embitz для STM32F103C8.
Вычисления синуса полиномами от 1-й до 6-й степеней оформлены в виде ассемблерных вставок в код на С. Перед каждой вставкой в комментариях есть примерно измеренное кол-во тактов на 1 цикл вычислений, когда таблица находится во flash и когда в RAM, время получается разное. Размешение таблиц в оперативной памяти позволяет немного поднять быстродействие.
Я не считаю себя знатоком ARMовского ассемблера, возможно кому-то удасться вычисления ускорить.

Здесь надо рассказать о том, почему использование коэффициента AC_SHIFT ускоряет вычисление. У Cortex-M3 есть инструкции умножения двух 32-битных чисел с получением 64-битного результата: UMULL и SMULL. Первая умножает два беззнаковых числа, и нам не подходит, потому коэффициенты есть положительные и отрицательные. Вторая, которую мы используем, умножает знаковые числа:
одно из них - это коэффициент из таблицы, знаковое число,
второе - это смещение угла внутри интервала, число беззнаковое.

Если второе число превысит 0x7fffffff, то оно блоком умножения будет интерпретировано как отрицательное и результат будет не тот, который мы ожидаем. Но инструкции, которая перемножала бы знаковое и беззнаковое число, в наборе команд нет. Мы будем использовать SMULL и позаботимся о том, что бы второе число не превышало 0x7fffffff, уменьшив его минимум в 2 раза. При этом надо будет после каждого умножения корректировать результат. Сделать это можно двумя способами:
- каждый раз удваивать результат, на что уйдут машинные такты;
- сделать в два раза больше то число, которое на него умножаем.

Более экономичным является второй способ. Если x уменьшен в 2 раза, то из формулы
y = (((0 + A3)·x + A2)·x + A1)·x + A0
понятно, что в 2 раза должен быть увеличен коэффициент A1, в 4 раза A2, и в 8 раз A3.
Если же x уменьшен в 4 раза, то множители становятся соответственно 4, 16 и 64. За это при компиляции таблиц отвечает параметр AC_SHIFT. Кроме того, его использование убивает второго зайца - коэффициенты, содержащие много нулей после запятой - а они относятся к высоким степеням полиномов, не так сильно теряют точность при округлении.

Вот и всё. Буду рад вашим замечаниям, уточнениям и предложениям.