Считаем синус быстро и точно: ч.1 - Теория

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

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


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


Кусочно-линейная аппроксимация — это, по сути, аппроксимация полиномом первой степени. Выходное значение вычисляется по формуле:
y = A1·x + A0
где x — это входной аргумент, A0, A1 — некие коэффициенты.
Для примера разобьём весь период от 0 до 2пи на 8 равных интервалов, и для каждого подберём такие коэффициенты A0 и A1, что бы максимум ошибки при любом аргументе был минимален (как это сделать, я расскажу ниже). Если мы построим график этого «синуса», он будет выглядеть примерно так (сплошная линия):

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

На нём, если приглядеться, не выполняется условие «максимум ошибки при любом аргументе был минимален», но сейчас это неважно. Максимальная ошибка составляет 0.03684497. Если найти логарифм этого числа по основанию 2, то получим -4.76239. Это можно назвать точностью вычисления в битах — она не хуже 4.76239 бит.
Теперь давайте разобьём один период (от 0 до 2π) на 64 интервала, по 16 на каждый квадрант. И рассмотрим один интервал в районе, например, 45°. График ошибки (сплошная чёрная линия) выглядит так:

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

y = A2·x² + A1·x + A0
Считаем тройку коэффициентов (об этом ниже), считаем целевой синус, и находим ошибку (сплошная чёрная линия):

Красным штрих-пунктиром поверх ошибки нарисован полином Чебышёва I-рода 3-й степени — графики снова очень похожи.
График синуса, найденный функцией 2-степени, с точностью до долей пикселя совпадает с «настоящим» синусом, и, будучи нарисованными вместе, они полностью перекрывали бы друг друга. Кстати, на картинке выше, где нарисован график с линейной интерполяцией на 8 отрезках, сиреневым пунктиром изображён синус, полученный функцией 2-й степени — он выглядит как «настоящий».

Как сильно увеличилась точность при переходе от 1-й ко 2-й степени? Для 64 интервалов на период она возросла с 10.7 до 17.63 бит. Откуда взялись эти числа, я расскажу во 2-й части.
Если мы дальше продолжим увеличивать степень полинома, то ошибки будут ещё меньше, и будут приобретать формы полиномов Чебышёва разных степеней:


На этих картинках нарисована ошибка при вычислении синуса полиномами 3-й и 4-й степеней.

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

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


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

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

Небольшое отступление, как представлять аргументы. На 32-битных процессорах полный период (круг) удобно представить как 2^32 (0x100000000), а угол, соответственно, в диапазоне от 0 до 0xffffffff. Количество интервалов аппроксимации, на которые разбит круг, равно 2^n, и угол можно интерпретировать как число с фиксированной точкой. Например, при разбиении на 64 интервала (2^6), угол будет интерпретироваться как число размерности 6.26. Здесь старшие 6 бит — это номер интервала (считая от нуля), а младшие 26 бит — смещение внутри него.
Возьмём, для примера, угол 15°. Если вычислить 15°/360° * 2^32, то получим 0000_1010_1010_1010_1010_1010_1010_1011 в двоичном виде. В формате с фикс. точкой 6.26 будет выглядеть как 000010.10101010101010101010101011. Здесь 000010 — это десятичное 2 (2-й интервал), а 10101010101010101010101011, будучи поделено на 2^26, даст 0.66666667 — это смещение внутри интервала. Или 2.66666667 = 64 * 15/360.
Смещение внутри интервала будет лежать в диапазоне 0 ≤ x < 1. К этому диапазону нужно привести корни многочленов. Корни полинома 2-й степени получатся 0.14645 и 0.85355, 3-й степени — 0.066987, 0.5 и 0.933013.

При аппроксимации полиномом 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π)
Для интервала №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^7 и 2^6 соответственно, и их значения всё равно будет лежать в пределах от -1 до +1. Интервал №15 я выбрал не случайно — на нём значение A0 максимально и близко к 1.
Вообще, в большинстве случаев коэффициенты при больших степенях можно увеличить на некий коэффициент. При целочисленных операциях это уменьшит погрешность вычислений, а на Cortex-M3 к тому же сильно сократит их время — об этом я расскажу в 3-й части.
Для вычисления таблиц других размеров в предыдущую систему уравнений вместо 64 нужно подставить нужный размер.
Для аппроксимации полиномом степени P нужно найти корни полинома Чебышёва степени P+1, и записать систему из P+1 уравнений с P+1 неизвестными, не забывая возводить корень многочлена в нужную степень 'n' при каждом An.
Если предыдущее предложение непонятно, то ничего страшного. В третьей части будет ссылка на готовый генератор таблиц и краткая инструкция к нему.

Часть 2 — Точность вычислений
Часть 3 — Практика
  • +2
  • 08 октября 2021, 14:39
  • vix

Комментарии (1)

RSS свернуть / развернуть
Очень интересно!!!
0
  • avatar
  • Aneg
  • 19 октября 2021, 08:47
Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.