Комментарии 41
Как минимум производительность C# тестировать через Stopwatch неверно.
Следует использовать https://benchmarkdotnet.org/
Также C# проверяет границы массива каждый раз, в отличии от C++.
Чтобы этого избежать можно взять указатель на массив или воспользоваться MemoryMarshal.GetArrayDataReference
Разве в простых for-циклах по диапазону оптимизатор не справляется с удалением этих проверок?
можно было сделать xf unsigned и заменить битовые операции на деление и остаток.
Вот только для x > max_int могут быть проблемы с этим методом, технически там тоже уб в этом пребразовнии. Практически если число такое большое что погрешность точности порядка пи, синус всё равно не посчитать.
Кстати implementation-defined там только для отрицательных пишут.
Интересный у вас Си, std::chrono нету, а std::cout есть )))
В этот статье написано, что в 8087 использовался CORDIC. В библиотечной функции тоже он?
Выкусите, те кто говорят, что математика не нужна в программировании! /раунд/
я последний раз про синус лет 30 назад слышал, в школе это было
Сам по себе скалярный синус не особо интересен. Интересно вычисление синуса в векторных регистрах. И тут реализаций огромное множество. Хотелось бы видеть сравнение с векторными библиотеками, например Intel SVML. Вот, например, реализация синуса под AVX-512:
https://github.com/numpy/SVML/blob/main/linux/avx512/svml_z0_sin_s_la.s
Впечатляет код?
Сравнивать с вручную написанным кодом на ассемблере — безнадёжное занятие.
Почему же? Это сравнение наглядно покажет, насколько бесперспективно написание велосипедов.
Мой личный рекорд на этом поприще
Мой личный рекорд: реализация FFT на ~10-15% быстрее, чем fftw для изображения определённого размера. Причём не на ассмеблере, а на C++ с интринсиками. Но времени убил очень много.
а вот алгоритмы типа FFT, свёртки, рекурсивных фильтров, численного интегрирования или перемножения матриц — самое то.
А там особо ничего, кроме как доступа к памяти и кэшу, оптимизировать и не получится. А сложные алгоритмы перемножения матриц с хорошей асимптотикой далеко не всегда дают хорошие результаты.
А там особо ничего, кроме как доступа к памяти и кэшу, оптимизировать и не получитсяТо был список алгоритмов, в которых важна не только скорость, но и точность, поскольку все они страдают накоплением погрешности. SSE/AVX можно использовать для повышения их точности без особого ущерба производительности — просто это менее очевидно, чем распараллеливание однотипных вычислений.
Минус полиномами Чебышева считал ещё ZX Spectrum basic, внутри своего "калькулятора"
Тут интересно порешать задачку, при сильных ограничениях в мощностях.
// начальная инициализация
cs = 1; // текущий косинус
sn = 0; // текущий синус
dx = pi / 65536; // приращение фазы
cosdx = cos(dx);
sindx = sin(dx);
// цикл генерации синусоиды
while (true)
{
tt = cs * cosdx - sn * sindx;
sn = sn * cosdx + cs * sindx;
cs = tt;
amp = sqrt(cs*cs+sn*sn);
cs /= amp;
sn /= amp;
}Правда, при такой реализации фаза неизбежно будет «плавать».
UPD: возможно, что можно увеличить устойчивость такого метода при использовании ещё и предыдущих значений.
Тут хватило бы билинейной интерполяции. Типа
F(x) = x[n] + K1*(x-x[n]) + K2*(x-x[n])^2
где x[n] ближайшее значение из таблицы
K1 можно взять или из аналогичной таблицы значений K1 посчитать как x[n+1] - x[n], линейная интерполяция
K2 взять из таблицы ранее рассчитанной, билинейная интерполяция
члены третьей и далее степени вряд ли нужны в данном контексте, при обычном ШИМ управлении двигателем, излишняя точность утонет в "шуме" общем системе.
Плюс таблица может быть без фиксированного шага по dx, в некоторых местах кривая аппроксимируется лучше, в других хуже, но это усложняет алгоритм, скорее всего необоснованно, метод интересен только теоретически.
Вот тут пример линейной интерполяции на С++ для ЦАП, подробно всё расписано, точности для практических применений хватает с запасом
1) используя многочлен низкой степени или таблицу с интерполяцией генерируем «грязную» синусоиду, с гармониками. Причём математически можно совершенно точно вычислить амплитуду всех гармоник (через разложение в ряд Фурье).
2) Пропускаем эту синусоиду через рекурсивный высокодобротный bandpass-фильтр, подавляющим первую гармонику до приемлемого уровня.
Для управления двигателем должно хватить. И его возможно даже можно улучшить, если минимизировать отклонения в частотной, а не во временной области.
Но зачем? Там же аппаратный баг у Intel, который так и не исправили. FPU FSIN и FCOS ошибаются очень сильно. Все пользуется кодом из glibc и sincos из AMD /Microsoft C library (MIT лицензия).
retrocomputing.stackexchange.com/questions/16523/intel-cpu-bug-in-the-90sfsincos
Там же аппаратный баг у Intel, который так и не исправилиБаг был у fdiv, а не fsincos, и его не только исправили, но ещё и процессоры всем поменяли.
FPU FSIN и FCOS ошибаются очень сильноА чуть более подробно можно? При каких значениях аргумента какого порядка ошибки возникают. Потому как в своих испытаниях ничего такого я не заметил.
randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion
FSIN уже не исправишь, так как это сломает ожидания целой пачки програмного обеспечения. AMD попытались, но пришлось убирать. (У них один и тот же код для это в микроинструкциях.)
В ссылке сверху все описано, но я могу дать ссылки на обсуждение в glibc.
sourceware.org/bugzilla/show_bug.cgi?id=13658
github.com/bminor/glibc/commit/b35fe25ed9b3b406754fe681b1785f330d9faf62
На x86_64 они используют sinf кстати.
Очень вероятно, что то что вы тестируете это как раз програмный код оттуда.
Я физик, поэтому бесконечно далек от таких проблем, но мне всегда казалось что для вычисления синуса логично просто решать задачу Коши x'' + x = 0, x(0) = 0, x'(0) = 1
Тем более там есть интеграл x^2 + x'^2 = 1, который позволит даже грубым методам удержать погрешность в разумных пределах.
Почему не вычисляют так?
Решать задачу коши Рунге Куттой? Это надо сделать сколько шагов, несколько десятков, сложнее чем просто подставить число в многочлен 13го порядка.
Ну а теперь попробуйте найти ответ на микроконтроллере типа Attiny 13, где 64 байта памяти, 1024 байт для программы. Для вычисления синуса можно выделить 200 байт на флеше и 8 байт оперативной памяти. Ну и уложиться желательно в 50-100 тактов процессора. Аппаратного умножения нет, деления нет.
Примерно такие же ограничения будут при формировании синусоиды на частоте 1000 МГц на ПЛИС, минимум логических элементов и высокие требования к быстродействию.
Взято из Quake III Arena?
Ну, тогда считаем корень, быстрее всех )))
float Q_rsqrt( float number )
{
const float x2 = number * 0.5F;
const float threehalfs = 1.5F;
union {
float f;
uint32_t i;
} conv = {number}; // member 'f' set to value of 'number'.
conv.i = 0x5f3759df - ( conv.i >> 1 );
conv.f *= threehalfs - x2 * conv.f * conv.f;
return conv.f;
}или вот так
float Q_rsqrt( float number )
{
long i;
float x2 ,y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y;
i = 0x5f3759df - ( i >> 1);
y = * ( float * ) &i;
y = y * (threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); //2nd iteration
return y;
}Есть и более правильная константа 0x5F375A86.
А если тронуть две другие константы, http://rrrola.wz.cz/inv_sqrt.html
При использовании многочленов Чебышева количество верных знаков ориентировочно соответствует степени многочлена. Для примера возьмём 7-ю степень — её вполне достаточно для задач визуализации:
Если под визуализацией понимать рисование на экране полной или почти полной синусоиды (с размахом вертикальной оси порядка 1), то 7 знаков - это явный перебор. Вполне хватит трех, максимум четырех. Пикселов- то у нас всего E+3. Нет смысла считать точнее.
Как посчитать синус быстрее всех на хабре