Comments 53
Функция smoothstep аппроксимирует синус с хорошей точностью.
Для своей предельной простоты и скорости - точность неплохая. Понятно, что это не для научных/технических расчетов, это для задач типа игровой анимации - точность на глаз достаточная, скорость больше некуда.
Решая уравнение 0 = cos(n*arccos(x)), находят корни полинома. У меня это сделано примерно так же (строки 149-162).
T val = cosl(M_PI * (0.5 + i) / n);Вообще, интересно, как устроено вычисление тригонометрической функции в математическом пакете или языке высокого уровня. Подозреваю, что это делается при помощи интерполяции полиномами (возможно, полиномами Чебышева).
Знатоком не являюсь, но по своему опыту могу сказать, что очень разными, зависит от платформы. На платформе x86 FPU вычисляет синус одной инструкцией fsin. На STM32 со встроенной плавающей точкой идёт много вызовов разных функций, сильно не вникал. На целочисленных CPU ещё сложнее.
На ПК, если ничего не путаю, считается первый квадрант полиномом с нечётными степенями.
Корни можно найти однажды, хранить в таблице и использовать постоянно.
Эта программа для создания таблиц на ПК, и быстродействие не важно, если вы об этом. Вычисление корней полинома нужной степени "на лету" мне кажется более элегантным и универсальным.
Знатоком не являюсь, но по своему опыту могу сказать, что очень разными, зависит от платформы. На платформе x86 FPU вычисляет синус одной инструкцией fsin.Я тоже не знаток, исхожу из здравого смысла. Даже если процессор с плавающей точкой реализует вычисление тригонометрической функции одной инструкцией, все равно где-то у него внутри должна быть микропрограмма (или что-то подобное, не знаю, как оно называется), которая сведет вычисления sin к операциям сложения и умножения чисел (и еще, возможно, поиска в предопределенной таблице). Поэтому, мне кажется, значения тригонометрических функций должны вычисляться не так уж и медленно, и порядок именно такой — sin вычисляется с помощью полинома, а не наоборот, и схема должна быть примерно такой, как вы описываете или похожей (повторюсь, это мое личное мнение)
Во-первых, в x86 ЦПУ есть просто инструкции синуса и косинуса.
Если встроенные инструкции не использовать, (например, для векторизации, для воспроизводимости, или когда инструкции недоступны), то используется некий полином. Например, вот что делает библиотека SLEEF: https://github.com/shibatch/sleef/blob/85440a5e87dae36ca1b891de14bc83b441ae7c43/src/libm/sleefdp.c#L789-L840
(сначала аргумент приводится к определенному диапазону, и затем используется полином 8 порядка)
Update: вот реализация из musl libc (в свою очередь, основана на реализации из freebsd):
static const double
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
double __sin(double x, double y, int iy)
{
double_t z,r,v,w;
z = x*x;
w = z*z;
r = S2 + z*(S3 + z*S4) + z*w*(S5 + z*S6);
v = z*x;
if (iy == 0)
return x + v*(S1 + z*r);
else
return x - ((z*(0.5*y - v*r) - y) - v*S1);
}
Ну, это потому, что здесь у него миллион корней
Но вот для небольших порядков, вроде единиц-десятков, считать их значения через косинусы мне кажется нецелесообразным. Наоборот (наоборот!), есть смысл вычислять косинусы с помощью полиномов Чебышева.
Не имея большого опыта работы с полиномами Чебышева, хотел тем не менее поинтересоваться — часто ли в реальной жизни нужно вычислять значения этих полиномов порядка миллион, да еще вне их естественной области определения? И какого порядка ошибки при этом возникают? В вашей последней публикации таки появились числа типа 0.9 + 10^-16i, которые показывают, что комплексные числа тоже могут вызывать проблемы (возможно, в этом конкретном случае без них не обойтись, но и с ними тоже приходится задумываться)
есть смысл вычислять косинусы с помощью полиномов Чебышева.В некоторых случаях — да, есть, и в моей последней статье такой случай тоже упоминается при вычислении окна Нуталла (причём результат получен аналитически, а не численно, и такой фокус тоже не всегда срабатывает). Так как вычисляет автор статьи — не факт, потому что есть методы поэффективнее (но спорить не буду).
Не имея большого опыта работы с полиномами Чебышева, хотел тем не менее поинтересоваться — часто ли в реальной жизни нужно вычислять значения этих полиномов порядка миллион, да еще вне их естественной области определения?Есть такая оконная функция — Дольфа-Чебышева, которая как раз считается через многочлены Чебышева, в том числе и за пределами диапазона (-1;1). А степень многочлена Чебышева в ней линейно зависит от количества отсчётов оконной функции, и окно длиной всего лишь в секунду может занимать 192000 отсчётов. До миллиона, как видите, совсем недалеко (а мой спектроанализатор поддерживает семплы до 8 миллионов отсчётов — что тоже не так уж и много при пересчёте на минуты).
Первое, что приходит в голову — записать рекуррентное соотношение для многочленов Чебышева в матричном виде и воспользоваться алгоритмом быстрого возведения в степень. Для миллионной степени будет примерно 20 итераций, что не так уж и много.
Ну это примерно как быстрый расчет чисел Фибоначчи (кажется, даже на Хабре когда-то был разбор этого метода).
Кстати, как следствие матричной записи, можно привести матрицу к диагональному виду и получить ответ в известном виде (правда, с использованием комплексных чисел):
Tn(x) = { [x+i√(1−x2)]n + [x−i√(1−x2)]n } / 2.
ChebyshevT[27, x] = ChebyshevT[3, ChebyshevT[3, ChebyshevT[3, x]]]
Так вычислять полиномы Чебышева имеет смысл для высоких степеней, что бы избежать ошибок вычисления высоких степеней и суммирования близких величин.
Ссылки на гит сбились, выдает 404
Разница в точности между ARM и x86 не может ли быть связана с разными ключами компиляции? Разные там strict / fast math? А свежие MSVC компиляторы, к примеру, вообще не умеют в 80-битный long double, у них всегда 64 бита. Какие компиляторы и настройки использовались?
Какие компиляторы и настройки использовались?
Везде gcc, настроек минимум. Я не являюсь знатоком компиляторов - меняю настройки, только если совсем собираться не хочет.
В разных проектах встречаются скрипты (*.cmd, *.sh) для компиляции, можете взглянуть.
А свежие MSVC компиляторы, к примеру, вообще не умеют в 80-битный long double, у них всегда 64 бита.
Хм. Не зря я их не люблю...
Для расчета референсных значений можно использовать boost::multiprecision
Спасибо!
Помимо максимальной ошибки, для тригонометрических функций часто важны и другие свойства:
Непрерывность
Дифференцируемость
Правильные значения при "специальных" значениях аргумента (0, 30, 45, 60, 90, ...)
Непревышение 1 модуля значения (чтобы не получить сюрпризы при `sqrt(1 - sin^2)`)
А чем не устроило разложение в ряд Тейлора? Вычисляется просто, итеративно, остаточный член O(x^2n)
Присоединяюсь к вопросу. В любом калькуляторе все функции через эти ряды считаются.
разложение в ряд Тейлора
А вот кстати...
Здесь автор этим и занимается, делая основной упор на точности вычислений с плавающей точкой (формула оттуда же):
Было бы интересно сравнить оба метода по точности и быстродействию на фиксированной точке. Возможно, займусь в обозримом будущем
чем не устроило разложение в ряд Тейлора
Подумал немного, и готов дать ответ.
Рядом Тейлора синус считают, емнип, только для первого квадранта, и там целевая функция непрерывна. В методе, описанном в статье, есть похожий случай, когда таблица состоит из 4 строк, и в этом случае аппроксимация полиномами Чебышёва может дать худший результат при том же количестве операций умножения (я исхожу из того, что степени Х в рядах Тейлора растут в 2 раза быстрее. Поправьте меня, если ошибаюсь).
Но при использовании рядов Тейлора я не вижу способа аппроксимации небольшого интервала, например от 110/256 до 111/256 части полного оборота. А полиномами Чебышёва это можно сделать, найдя корни именно для этого промежутка. Такое сужение диапазона аргумента даст прирост точности, что можно конвертировать в уменьшение кол-ва математических операций при снижении точности до нужного значения.
Но сравнение обоих методов на 1-м квадранте, целыми числами, я всё-таки сделаю. Наверное...
<Рядом Тейлора синус считают, емнип, только для первого квадранта, и там целевая функция непрерывн...>
Вы неправильно помните.
Разложение точки в окрестности 0 это только частный случай ряда Тейлора - ряд Маклорена. Ряд Тейлора применим для окрестности любой точки, а не только 0.
А почему нет сравнения с другими методами аппаратного вычисления приближения синуса, например, CORDIC (кстати там выше ссылка была в одном из комментариев)?
У последнего, кстати, очень удобная интерпретация вычисления в битах, т.к., собственно, каждый шаг уточнения является сложением и делением на два (сдвиг на бит).
Даже если тот же CORDIC не может дать нужную точность или, например, на МК считается плохо (на ПЛИС, кстати, очень даже хорошо работает, т.к. там очень хорошее целочисленное представление получается), то всё равно можно было бы уделить этому хотя бы пару абзацев.
Кстати, по поводу МК (пусть и на ARM-ядре) - уже сейчас есть STM32G4 с CORDIC-сопроцессором на борту (https://www.compel.ru/lib/140072).
А куда же делись 10 бит, я так и не понял. Есть подозрение, что точность, равная разрядности мантиссы на 32- и 64-битных float получалась за счёт запаса 80-битных вычислений, результаты которых потом огрублялись до заданных 32 и 64, а без огрубления она как раз составляла те самые 54 бита. При вычислениях на ARMах с 64-битным FPU точность вычисления упала до 50.8 бита против 53 на x86/x64.
правильно ли я понимаю, что у вас все исходные данные были с 80-бит точностью и результат был доступен в формате 80-бит, но точность получилась все равно duble? Была на хабре очень интересная статья про синус-косинус и там приводился пример пакета для вычислений с double-double и quad-double точностью. Считает медленно, но если нужен более точный результат, то очень хорошая штука.
правильно ли я понимаю, что у вас все исходные данные были с 80-бит точностью и результат был доступен в формате 80-бит
Нет. Исходные данные, как и промежуточные результаты заданы либо 32, либо 64, либо 80 бит float, зависит от значения TABLE_TYPE.
Но внутри FPU все вычисления производятся в формате 80 бит. И я не уверен, что компилятор (точнее код, им сгенерённый) гоняет промежуточные результаты между FPU и памятью, огрубляя их. То есть, независимо от типа TABLE_TYPE, все вычисления скорее всего производятся как 80-битные, и только в самом конце происходит огрубление (приведение к типу).
class Program
{
static double ChebyshevSin(double x)
{
double xx = x * x;
return
x*(1.570796326794895 +
xx*(-0.6459640975061731 +
xx*(0.07969262624514188 +
xx*(-0.004681754128869211 +
xx*(0.0001604411632629485 +
xx*(-3.598802460206541e-6 +
xx*(5.687768400923302e-8 -
6.434815323635176e-10*xx)))))));
}
static void Main(string[] args)
{
double max = 0.0;
int count = 100000000;
for (int i = 0; i < count; i++)
{
double x = i / (double)count;
double y = Math.Sin(x * Math.PI / 2.0) - ChebyshevSin(x);
max = Math.Max(max, Math.Abs(y));
}
Console.WriteLine("max error: {0}",max);
Console.ReadKey();
}
}Максимальное отклонение составило 5.55*10-16, то есть как минимум 15 знаков после запятой совпадают. Также на SSE код можно ускорить почти в два раза, если разбить многочлен на два и считать параллельно, а результаты сложить.
static double ChebyshevSin(double x)
{
double xx = x * x;
return
x * (0.9999999999999999980 +
xx * (-0.1666666666666666190 +
xx * (0.00833333333333299283 +
xx * (-0.0001984126984115937723 +
xx * (2.75573192045707547e-6 +
xx * (-2.50521063803320390e-8 +
xx * (1.605891861118114436e-10 +
xx * (-7.64251091766084203e-13 +
2.71665816482462532e-15 * xx))))))));
}библиотечная: 50.6 бит,
Чебышевым: 50.45 бит.
Разницу считал тоже в Wolfram-е.
А есть ли у кого вычисленная таблица синусов и косинусов в формате C++ include file? В градусах, от 0 до 360
Как посчитать синус быстро