Пишем numpy-модуль для ускорения математических функций с помощью SIMD-инструкций

[Imposeren]

странно, что в numpy этого нет «из коробки».

Насколько я знаю, то обвязку для C можно еще делать на Cython и синтаксис там будет привычнее для тех кто привык к python (как-то так наверное). Не думали глянуть такой вариант привязки? Ну а может и ваш и этот вариант, да сравнить потом.

[nikolaynag]

странно, что в numpy этого нет «из коробки».

Согласен, но тут видимо проблема архитектурной зависимости. Все-таки не на всех процессорах поддерживается AVX, а реализовать детектирование поддерживаемых SIMD-интрукций и выбор оптимальных функций для каждого типа данных при импорте модуля достаточно сложно.

По поводу Cython — его наверное можно было бы использовать, но тут не просто вызов С-функции из python, а еще и работа по упаковыванию элементов numpy-массива в группы по 4 элемента и последующая распаковка обратно. В любом случае на Cython вряд ли получится быстрее, чем на C :)

[Lockal]

Точность операций над числами с плавающей точкой нельзя оценить одним числом. На то и называется точка плавающей, что посчитать sin(2e20 * pi) возможно, но абсолютно бессмысленно.

[nikolaynag]

Конечно, абсолютная точность у sin(2e20*pi) и sin(2*pi) будет разной (как впрочем она будет различаться для самих чисел double, содержащих 2*pi или 2e20*pi):

xsin(1e+00*M_PI) = 1.22465e-16, sin(1e+00*M_PI) = 1.22465e-16
xsin(1e+01*M_PI) = -1.22465e-15, sin(1e+01*M_PI) = -1.22465e-15
xsin(1e+02*M_PI) = 1.96439e-15, sin(1e+02*M_PI) = 1.96439e-15
xsin(1e+03*M_PI) = -3.21417e-13, sin(1e+03*M_PI) = -3.21417e-13
xsin(1e+04*M_PI) = -4.85682e-13, sin(1e+04*M_PI) = -4.85682e-13
xsin(1e+05*M_PI) = -3.39607e-11, sin(1e+05*M_PI) = -3.39607e-11
xsin(1e+06*M_PI) = -2.23191e-10, sin(1e+06*M_PI) = -2.23191e-10xsin(1e+10*M_PI) = -8.96779e+183, sin(1e+10*M_PI) = -2.23936e-06
xsin(1e+07*M_PI) = 5.62056e-10, sin(1e+07*M_PI) = 5.62056e-10
xsin(1e+08*M_PI) = -3.90829e-08, sin(1e+08*M_PI) = -3.90829e-08
xsin(1e+09*M_PI) = -3.32014e-08, sin(1e+09*M_PI) = -3.32014e-08

Именно потому в статье я и привожу относительную погрешность в сравнении с обычным синусом (т.е. (xsin(x)-sin(x))/sin(x)). Показатель конечно не идеальный (а при нулевом знаменателе и вовсе неопределенный), но представление о погрешности дает. С большими числами есть другая пробелема — недавно я обнаружил, что синус в SLEEF для чисел порядка 1e10 и больше, дает расходимость:

xsin(1e+10*M_PI) = -8.96779e+183, sin(1e+10*M_PI) = -2.23936e-06

Так что Вы правы насчет того, что с большими числами нужно быть осторожнее.