Спасибо, Артём. Я уважительно отношусь к вашему времени, поэтому буду краток. Действительно в следующей статье я подробно расскажу метод, который вы обрисовали сейчас. Вот ссылка на код, который будет обсуждаться в следующей статье libm sin.
Так же я Вам очень благодарен в том, что вы единственный, кто предоставил ссылку на код. Давайте же его посмотрим вместе. А именно на 126 строчку вот этого файла sinq_kernel.c. Я понимаю, что у вас мало времени, но, пожалуйста, потратьте ещё пару минут и посмотрите на него и одновременно посмотрите на 2-ую формулу моей статьи. Они одинаковы! Просто точь-в-точь!
Конечно, напишу немного больше. Спасибо за комментарий.
Пусть ypi истинное значение синуса, вычисленное в точке xi. Истинным назовём значение вычисленное с помощью функции sin(float128). yni — значение полученное с помощью представленной функции sin_e1-5. N — общее количество точек.
1) Среднее ошибки SUM(yni — ypi)/N
2) STD Среднеквадратичное отклонение величины (wiki, вторая сверху формула)
Введём ещё переменную yoi. Это значение синуса, вычисленное с помощью встроенной функции библиотеки glibc sin(double).
Представьте, если значение yoi и yni не совпадают. Тогда (цитирую код функции).
if( <результат встроенной и нашей функции не совпадает побитово >) {
if( fabsq(q - dold.x) <= fabsq(q - dnew.x) )
co++;
else
cn++;
}
увеличиваем счётчик правильности функции для встроенной или нашей функции. Где q это ypi
3) Лучше наша. cn / N
4) Лучше libm. co / N
Я буду очень благодарен, если вы напишите своё мнение о понятности этого комментария. Если всё ОК, то я добавлю это в исходную статью.
Забегу вперёд. В компиляторах GNU есть опция fast-math. При её включении используется код, который считает синус очень быстро, но с точностью 4 ULP. glibc source
Я подробно буду обсуждать это в одной из следующих статей.
Спасибо за развёрнутый комментарий. Действительно интересно прочитать другую точку зрения. Я на самом деле думаю, что если вы добавите ссылки и примеры кода для разных методов и т.п., то это может быть очень интересная статья на хабре. По крайней мере я бы её с удовольствием бы сам прочитал. Пока, что некоторые из ваших утверждений кажутся мне не совсем очевидными.
1) Если дословно понимать утверждение 4 то это очень плохой метод. Интерполяция между табличными значениями даёт очень плохую точность при больших затратах времени расчёта. Я думаю, что вы имели ввиду не интерполяцию, а разложение в ряд Тейлора вокруг табличных значений. Этот метод требует существенно меньше памяти, чем кубическая интерполяция, а точность даёт несравненно выше. О нём я расскажу в следующей статье.
2) Рекуррентные формулы не имеют практической ценности. Даже если абсолютно убрать накопление ошибки, я правда не знаю таких методов, то их использование целесообразно только в таком сценарии: сначала рассчитывается таблица, а потом по этой таблице берутся значения sin. Это таблицу всегда можно/нужно предрасчитать перед началом работы. Если у вас в программе синус считается от силы раз 10, то смысл предрасчитывать таблицу. А если триллионы раз, то ограничение по производительность будет как раз доступ к этой таблице (извлечение из памяти) и т.п., а не её расчёт в начале программы.
3) Ряд тейлора или полиномы Чебышева/минимакс-полиномом это всё равно полиномы. Принципиальной разницы быть не должно в самом коде. Если вы думаете по-другому — пишите.
4) CORDIC. Просто по определению метода угол вращения должен быть arctg(1/(2^(i−1))) Т.е. предполагая, что для малых x arctan(x) ~= x мы получим что нужно 54 шага, чтобы выйти из точности double. Очевидно, что ряд Тейлора сходится быстрее.
При написании статьи я внимательно просмотрел исходный код нескольких математических библиотек. Выбрал glibc и пишу свои статьи на основе практик оттуда. В следующей статье я буду давать конкретные ссылки на исходный код. Если вы считаете что libm от GNU это плохая практика и вы можете сделать лучше, то я призываю Вас не ограничиваться комментарием к моей статье а написать свою. Это всем будет интересно.
Ну и напоследок, повторю, что написал в начале статьи. Она предназначена для того, чтобы показать «хорошие» не очевидные практики работы с числами с плавающей точкой. Я считаю, что «I did my best».
Спасибо за вопрос. Действительно надо было пару слов про это сказать. Под катом в начале статьи есть полная программа. Набор данных — линейное заполнение 16 млн числами от 0 до pi/2.
vector<double> xv;
xv.resize(8 * 2000000);
// Линейное заполнение массива значениями от 0 до M_PI/2
for (int i = 0; i < xv.size(); i++) {
// Максимальное значение в массиве изменять здесь.
xv[i] = (M_PI / 2) * i / double(xv.size());
}
Хороший вопрос.
Во-первых хочу заметить, что одна из целей статьи показать, что от «перемены мест слагаемых сумма меняется» и показать методы улучшения точности. Если можно вычислить точнее в одно и тоже время, почему так не сделать?
Ну а если говорить по синус, то частый пример — численное дифференцирование функции. Оно требуется в решении дифференциальных уравнений и оптимизации функции. Самый простой пример. cos(x) = (sin(x + dx) — sin(x)/dx. Ясно, что при уменьшении dx, ошибка в последнем знаке будет давать гораздо большую ошибку функции cos. Вычислении численного градиента, например, используется при минимизации функций. Минимизация с неточным градиентом может потребовать больше шагов и/или привести к неточному решению. Вычисление производных более высокого порядка, которое может быть нужно для численного решения дифференциальных уравнений, даст ещё большую ошибку. К тому же для многих уравнений ошибка накапливается и спустя какое-то время становится существенной.
Из самых простых примеров дифференциальных уравнений — это обычный маятник. Только в первом приближении уравнение его движения линейно (Маятник).
Если Вас не интересует точность double, то, конечно, можно использовать числа float. Для них все перечисленные методы работают абсолютно так же. А последний знак там всего лишь 0.000 000 119
Так же я Вам очень благодарен в том, что вы единственный, кто предоставил ссылку на код. Давайте же его посмотрим вместе. А именно на 126 строчку вот этого файла sinq_kernel.c. Я понимаю, что у вас мало времени, но, пожалуйста, потратьте ещё пару минут и посмотрите на него и одновременно посмотрите на 2-ую формулу моей статьи. Они одинаковы! Просто точь-в-точь!
Пусть ypi истинное значение синуса, вычисленное в точке xi. Истинным назовём значение вычисленное с помощью функции sin(float128). yni — значение полученное с помощью представленной функции sin_e1-5. N — общее количество точек.
1) Среднее ошибки SUM(yni — ypi)/N
2) STD Среднеквадратичное отклонение величины (wiki, вторая сверху формула)
Введём ещё переменную yoi. Это значение синуса, вычисленное с помощью встроенной функции библиотеки glibc sin(double).
Представьте, если значение yoi и yni не совпадают. Тогда (цитирую код функции).
увеличиваем счётчик правильности функции для встроенной или нашей функции. Где q это ypi
3) Лучше наша. cn / N
4) Лучше libm. co / N
Я буду очень благодарен, если вы напишите своё мнение о понятности этого комментария. Если всё ОК, то я добавлю это в исходную статью.
Я подробно буду обсуждать это в одной из следующих статей.
1) Если дословно понимать утверждение 4 то это очень плохой метод. Интерполяция между табличными значениями даёт очень плохую точность при больших затратах времени расчёта. Я думаю, что вы имели ввиду не интерполяцию, а разложение в ряд Тейлора вокруг табличных значений. Этот метод требует существенно меньше памяти, чем кубическая интерполяция, а точность даёт несравненно выше. О нём я расскажу в следующей статье.
2) Рекуррентные формулы не имеют практической ценности. Даже если абсолютно убрать накопление ошибки, я правда не знаю таких методов, то их использование целесообразно только в таком сценарии: сначала рассчитывается таблица, а потом по этой таблице берутся значения sin. Это таблицу всегда можно/нужно предрасчитать перед началом работы. Если у вас в программе синус считается от силы раз 10, то смысл предрасчитывать таблицу. А если триллионы раз, то ограничение по производительность будет как раз доступ к этой таблице (извлечение из памяти) и т.п., а не её расчёт в начале программы.
3) Ряд тейлора или полиномы Чебышева/минимакс-полиномом это всё равно полиномы. Принципиальной разницы быть не должно в самом коде. Если вы думаете по-другому — пишите.
4) CORDIC. Просто по определению метода угол вращения должен быть arctg(1/(2^(i−1))) Т.е. предполагая, что для малых x arctan(x) ~= x мы получим что нужно 54 шага, чтобы выйти из точности double. Очевидно, что ряд Тейлора сходится быстрее.
При написании статьи я внимательно просмотрел исходный код нескольких математических библиотек. Выбрал glibc и пишу свои статьи на основе практик оттуда. В следующей статье я буду давать конкретные ссылки на исходный код. Если вы считаете что libm от GNU это плохая практика и вы можете сделать лучше, то я призываю Вас не ограничиваться комментарием к моей статье а написать свою. Это всем будет интересно.
Ну и напоследок, повторю, что написал в начале статьи. Она предназначена для того, чтобы показать «хорошие» не очевидные практики работы с числами с плавающей точкой. Я считаю, что «I did my best».
Во-первых хочу заметить, что одна из целей статьи показать, что от «перемены мест слагаемых сумма меняется» и показать методы улучшения точности. Если можно вычислить точнее в одно и тоже время, почему так не сделать?
Ну а если говорить по синус, то частый пример — численное дифференцирование функции. Оно требуется в решении дифференциальных уравнений и оптимизации функции. Самый простой пример. cos(x) = (sin(x + dx) — sin(x)/dx. Ясно, что при уменьшении dx, ошибка в последнем знаке будет давать гораздо большую ошибку функции cos. Вычислении численного градиента, например, используется при минимизации функций. Минимизация с неточным градиентом может потребовать больше шагов и/или привести к неточному решению. Вычисление производных более высокого порядка, которое может быть нужно для численного решения дифференциальных уравнений, даст ещё большую ошибку. К тому же для многих уравнений ошибка накапливается и спустя какое-то время становится существенной.
Из самых простых примеров дифференциальных уравнений — это обычный маятник. Только в первом приближении уравнение его движения линейно (Маятник).
Если Вас не интересует точность double, то, конечно, можно использовать числа float. Для них все перечисленные методы работают абсолютно так же. А последний знак там всего лишь 0.000 000 119