Считаем синус быстро и точно: ч.2 - Точность вычислений

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



Разные комбинации размера раблицы и степени полинома дают разную погрешность. Не имея представления, как определить её аналитически, я решил считать «в лоб», перебирая все возможные комбинации, и прогоняя через каждую все углы от 0 до 0xffffffff. Для этого был создан проект accuracy_test. Суть его в следующем:
  • в него включены несколько (90) таблиц, степенью полинома от 1 до 6 и размером от 4 до 65536. В аргументах при запуске можно указать набор или диапазон проверяемых таблиц;
  • с помощью заданной таблицы производится аппроксимация синуса по всем возможным 4294967296 аргументам. Результат аппроксимации сравнивается со значением, полученным стандарной функцией sinl(), получается разница-ошибка;
  • находится самая максимальная в абсолютном значении ошибка — это самый худший случай аппроксимации.


Количество аргументов от 0 до 2^32 велико, но не бесконечно, и их перебор занимает какое-то конечное время. Мой ноутбук с Intel Core-I7 перебирает их примерно за 6 минут, а одноплатному компьютеру на Allwinner A20 на это требуется 25 минут. Причём на ПК, как показал один эксперимент, львинную часть времени занимает вычисление «эталонного» значения функцией sinl.
Результаты проверок всех 90 таблиц в конечном итоге сводятся в одну небольшую табличку, на основании которой и построена диаграмма в начале статьи.
На ней видно, что с увеличением размера таблицы точность растёт, и тем быстрее, чем выше степень полинома, при этом не поднимаясь выше 54 бит. Это самая большая точность, которую мне удалось достичь, при этом вычисления производились на ПК (что немаловажно), с помощью 80-битных чисел с плавающей точкой (long double).

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

Здесь точность не превышает 30.36 бит. Но с учётом того, что старший бит отведён под знак числа, ошибка получается менее 1 бита.
На этой диаграмме присутствуют не все комбинации размера таблицы и степени полинома. Это связано с тем, что некоторые коэффициенты при малом количестве интервалов превышают 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 лучше не использовать.
Результаты всех проверок сведены в несколько файлов:
bit_accuracy_FLOAT_ALL.ods — для чисел с плавающей точкой;
bit_accuracy_INT32_ALL.ods — для целых 32-битных чисел;
bit_accuracy_FLOAT_E8_ALL.ods — для чисел с плавающей точкой, когда в качестве эталона использована не sinl, а аппроксимация теми же таблицами.
О последнем надо сказать особо.
Сравним числа с плавающей точкой и какие точности на них достигнуты:
  • float 32-бит, мантисса 23+1 бит, точность аппроксимации — 24 бита.
  • float 64-бит, мантисса 52+1 бит, точность аппроксимации — 53 бита.
  • float 80-бит, мантисса 63+1 бит, точность аппроксимации — почти 54 бита, хотя можно было ожидать на 10 больше.
Куда делись 10 бит?
Возникло предположение, что неточно вычисляется «эталонный» синус функцией sinl. Для проверки был проведён эксперимент, где сравнивались результаты аппроксимации двумя разными таблицами. В качестве эталона использовались несколько таблиц длиной 2^15 и 2^16, разных степеней полинома, а вычисления производились числами 80-битными long double.
Предположение не подтвердилась, зато было замечено, что время проверки уменьшилось в несколько раз. Напрмер, проверка таблицы 65536_p^1 функцией sinl заняла 279с, а той же функцией get_sine_float — 56с. Получается, более 80% времени процессор был занят вычислением функции sinl. В свете этого можно рекомендовать использовать на ПК аппроксимацию вместо sinl, когда время выполнения критично.
А куда делись 10 бит, я так и не понял. Есть подозрение, что точность, равная разрядности мантиссы на 32- и 64-битных float получалась за счёт запаса 80-битных чисел, в которых производилось вычисление, а их результат потом огрублялся до заданных. При вычислениях на ARMах с 64-битным FPU точность вычисления упала до 50.8 бита.

Альтернативные методы аппроксимации


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

Теперь более подробно, если кому интересно.

С использованием рядов Маклорена синус считается по следующей формуле:



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

Вынесем повторяющиеся фрагменты за скобки, например так:



или так:



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

Далее найдём значения синуса на интервале от о до 2пи обоими методами, и построим их графики вместе с графиком «настоящего» синуса:


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


Видно, что ошибка приближения полиномом Чебышёва в разы больше таковой у ряда Маклорена.

Теперь давайте посмотрим, как ошибки аппроксимации зависят от количества членов аппроксимирующей функции.

Графики ошибок аппроксимации рядами Маклорена разной длины выглядят так:



Графики ошибок аппроксимации полиномами Чебышёва имеют следующий вид:



Графики ошибок выглядят по разному, но в обоих случаях величина ошибки падает с увеличением количества членов ряда.

На следующей диаграмме показана сравнительная точность аппроксимации обоими методами в зависимости от количества членов ряда (полинома):



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

Однако для аппроксимации полиномами Чебышёва полный период (от 0 до 2пи) можно разбить не на 4, как это делалось во всех предыдущих вычислениях, а, например, на 512 равных интервалов. На диаграмме этот случай представлен серой линией. Точность в 53.94 бита при 6 членах (сответствует полиному 5-й степени) превышает таковую у ряда Маклорена в 51.62 бита при 10 членах ряда. Справедливости ради надо отметить, что при этом страдает такое свойство функции, как непрерывность, т.к. на границах интервалов будут разрывы, что в некоторых случаях может быть критичным.

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

Часть 1 — Теория
Часть 3 — Практика
  • +8
  • 08 октября 2021, 14:43
  • vix

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

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