Комментарии 58
Круто конечно, но где нужно на практике знать синус с такой потрясающей точностью?
Хороший вопрос.
Во-первых хочу заметить, что одна из целей статьи показать, что от «перемены мест слагаемых сумма меняется» и показать методы улучшения точности. Если можно вычислить точнее в одно и тоже время, почему так не сделать?
Ну а если говорить по синус, то частый пример — численное дифференцирование функции. Оно требуется в решении дифференциальных уравнений и оптимизации функции. Самый простой пример. cos(x) = (sin(x + dx) — sin(x)/dx. Ясно, что при уменьшении dx, ошибка в последнем знаке будет давать гораздо большую ошибку функции cos. Вычислении численного градиента, например, используется при минимизации функций. Минимизация с неточным градиентом может потребовать больше шагов и/или привести к неточному решению. Вычисление производных более высокого порядка, которое может быть нужно для численного решения дифференциальных уравнений, даст ещё большую ошибку. К тому же для многих уравнений ошибка накапливается и спустя какое-то время становится существенной.
Из самых простых примеров дифференциальных уравнений — это обычный маятник. Только в первом приближении уравнение его движения линейно (Маятник).
Если Вас не интересует точность double, то, конечно, можно использовать числа float. Для них все перечисленные методы работают абсолютно так же. А последний знак там всего лишь 0.000 000 119
Во-первых хочу заметить, что одна из целей статьи показать, что от «перемены мест слагаемых сумма меняется» и показать методы улучшения точности. Если можно вычислить точнее в одно и тоже время, почему так не сделать?
Ну а если говорить по синус, то частый пример — численное дифференцирование функции. Оно требуется в решении дифференциальных уравнений и оптимизации функции. Самый простой пример. cos(x) = (sin(x + dx) — sin(x)/dx. Ясно, что при уменьшении dx, ошибка в последнем знаке будет давать гораздо большую ошибку функции cos. Вычислении численного градиента, например, используется при минимизации функций. Минимизация с неточным градиентом может потребовать больше шагов и/или привести к неточному решению. Вычисление производных более высокого порядка, которое может быть нужно для численного решения дифференциальных уравнений, даст ещё большую ошибку. К тому же для многих уравнений ошибка накапливается и спустя какое-то время становится существенной.
Из самых простых примеров дифференциальных уравнений — это обычный маятник. Только в первом приближении уравнение его движения линейно (Маятник).
Если Вас не интересует точность double, то, конечно, можно использовать числа float. Для них все перечисленные методы работают абсолютно так же. А последний знак там всего лишь 0.000 000 119
P.S. Большая точность в расчёте sin далеко не всегда. Тему быстрых неточных вычислений я раскрою в одной из следующих статей.
Спасибо, буду ждать :)
А насчет точности, просто со времен курсового проекта по термодинамике приучили, что хватает четырех значащих цифр после запятой, вот и проецирую на всё.
А насчет точности, просто со времен курсового проекта по термодинамике приучили, что хватает четырех значащих цифр после запятой, вот и проецирую на всё.
нужна не всегда?
Забегу вперёд. В компиляторах GNU есть опция fast-math. При её включении используется код, который считает синус очень быстро, но с точностью 4 ULP. glibc source
Я подробно буду обсуждать это в одной из следующих статей.
Я подробно буду обсуждать это в одной из следующих статей.
Там где может происходить накопление ошибки, ну, например, суммирование таких вот синусов. Там нужно контролировать накопленную ошибку и минимизировать ее.
А набор данных на котором среднее отклонение и СКО считали он что из себя представляет?
Спасибо за вопрос. Действительно надо было пару слов про это сказать. Под катом в начале статьи есть полная программа. Набор данных — линейное заполнение 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());
}
deleted
Автор, мне кажется, что зря вы взяли именно синус. Надо было взять какую-нибудь экзотическую функцию, для которой нет или мало известных оптимальных решений.
Почему? Потому, что его вычисление через ряд Тейлора — далеко не оптимальный метод. Есть и другие методы, такие, как:
1) Разложение в ряд по полиномам Чебышева
2) Приближение минимакс-полиномом
3) CORDIC
4) Различные варианты, основанные на табличной интерполяции
Также, если нам нужны синусы не произвольных углов, а, допустим, надо рассчитать таблицу синуса по равноотстоящим углам — то для этого есть специальные методы, основанные на рекуррентных соотношениях. Например:
y[n] = a*y[n-1] — y[n-2]
Результатом применения этого рекуррентного соотношения являются значения функции синус (или косинус) на равноотстоящих углах x[i] = x[i-1]+d.
Разумеется, в этой формуле накапливаются ошибки округления, поэтому придуманы более продвинутые рекуррентные формулы, гарантирующие стабильность, имеющие меньшие погрешности и т.д.
В общем, в области синуса уже много что копано и перекопано. Чтобы пытаться что-то там улучшить, надо сначала изучить то, что уже сделано.
Почему? Потому, что его вычисление через ряд Тейлора — далеко не оптимальный метод. Есть и другие методы, такие, как:
1) Разложение в ряд по полиномам Чебышева
2) Приближение минимакс-полиномом
3) CORDIC
4) Различные варианты, основанные на табличной интерполяции
Также, если нам нужны синусы не произвольных углов, а, допустим, надо рассчитать таблицу синуса по равноотстоящим углам — то для этого есть специальные методы, основанные на рекуррентных соотношениях. Например:
y[n] = a*y[n-1] — y[n-2]
Результатом применения этого рекуррентного соотношения являются значения функции синус (или косинус) на равноотстоящих углах x[i] = x[i-1]+d.
Разумеется, в этой формуле накапливаются ошибки округления, поэтому придуманы более продвинутые рекуррентные формулы, гарантирующие стабильность, имеющие меньшие погрешности и т.д.
В общем, в области синуса уже много что копано и перекопано. Чтобы пытаться что-то там улучшить, надо сначала изучить то, что уже сделано.
Спасибо за развёрнутый комментарий. Действительно интересно прочитать другую точку зрения. Я на самом деле думаю, что если вы добавите ссылки и примеры кода для разных методов и т.п., то это может быть очень интересная статья на хабре. По крайней мере я бы её с удовольствием бы сам прочитал. Пока, что некоторые из ваших утверждений кажутся мне не совсем очевидными.
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».
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».
1) Интерполяция бывает не только линейной. На практике сплайн 3–4 порядка на таблице из пары десятков точек дает неплохие результаты.
2) Рекуррентные формулы широко используются в разного рода синтезаторах, где надо считать и комбинировать большое количество синусов.
3) Раз вам важно только то, что это полиномы, то не стоит говорить про ряд Тейлора, ибо он почти никогда не применяется для целей аппроксимации функций.
4) CORDIC предназначен для вычислений без использования операции умножения и на соответствующем оборудовании он обойдет вычисления по полиномам (не говоря уже о никуда не годном ряде Тейлора). Но да, на современном железе CORDIC безнадежно устарел.
На мой взгляд, задача на практическое вычисление функций обязана начинаться с четкого определения требуемой точности. Например, если нам нужен синус для рисования кругов на экране, то считать его до последнего ulp в double полностью бессмысленно, аппроксимации с точностью ±0.0001 будет более чем достаточно. Соответственно, постановка в виде «считать как можно быстрее с как можно большей точностью» является некорректной. Более того, расчет с точностью ±1ulp обычно нужен только разработчикам библиотек, никакие практические цели такой точности не требуют. Таким образом, если возникла проблема «что-то это считает слишком медленно, надо ускорить», то первое, что надо сделать — это уменьшить точность.
В общем, пока статья не тянет на «хорошие практики работы с плавающей точкой», скорее на набор частых заблуждений. Ни слова не сказано про выбор точности, хотя это является определяющим фактором финальной производительности. Использован ряд Тейлора без надписи жирными буквами, что это только для примера, для продакшена его использовать ни в коем случае нельзя.
2) Рекуррентные формулы широко используются в разного рода синтезаторах, где надо считать и комбинировать большое количество синусов.
3) Раз вам важно только то, что это полиномы, то не стоит говорить про ряд Тейлора, ибо он почти никогда не применяется для целей аппроксимации функций.
4) CORDIC предназначен для вычислений без использования операции умножения и на соответствующем оборудовании он обойдет вычисления по полиномам (не говоря уже о никуда не годном ряде Тейлора). Но да, на современном железе CORDIC безнадежно устарел.
На мой взгляд, задача на практическое вычисление функций обязана начинаться с четкого определения требуемой точности. Например, если нам нужен синус для рисования кругов на экране, то считать его до последнего ulp в double полностью бессмысленно, аппроксимации с точностью ±0.0001 будет более чем достаточно. Соответственно, постановка в виде «считать как можно быстрее с как можно большей точностью» является некорректной. Более того, расчет с точностью ±1ulp обычно нужен только разработчикам библиотек, никакие практические цели такой точности не требуют. Таким образом, если возникла проблема «что-то это считает слишком медленно, надо ускорить», то первое, что надо сделать — это уменьшить точность.
В общем, пока статья не тянет на «хорошие практики работы с плавающей точкой», скорее на набор частых заблуждений. Ни слова не сказано про выбор точности, хотя это является определяющим фактором финальной производительности. Использован ряд Тейлора без надписи жирными буквами, что это только для примера, для продакшена его использовать ни в коем случае нельзя.
Вы во многом правы в этой фразе:
Более того, расчет с точностью ±1ulp обычно нужен только разработчикам библиотек, никакие практические цели такой точности не требуют
но не во всём. Я как раз разработчик библиотек. И мне бывает нужна точность, которая на десятки (иногда сотни) порядков выше 1ulp по причине, указанной в этой статье. Один из практических смыслов такой точности — чтобы округление результата было корректным до последнего бита, чтобы обеспечить полную повторяемость результата не зависимо ни от компилятора, ни от вычислительной машины. А в обыденной жизни да, я согласен, что это не нужно, из-за последнего бита даже ракета не упадёт… хотя, не знаю. Знаю только, что из-за разницы даже в последнем бите бывает, что некоторые программы работают слишком по-разному, но это из-за наличия в них кучи других ошибок. Ошибка в округлении последнего бита является лишь катализатором для их проявления.
CORDIC предназначен для вычислений без использования операции умножения и на соответствующем оборудовании он обойдет вычисления по полиномам (не говоря уже о никуда не годном ряде Тейлора). Но да, на современном железе CORDIC безнадежно устарел.
CORDIC может быть и устарел, а сам метод Shift-and-Add, лежащий в его основе — нет. Сейчас вообще всё идёт по направлению комбинации всех существующих методов. В зависимости от ситуации, каждый метод может оказаться хорошим, поэтому однозначно говорить, что уже можно сбросить со счетов, а что нельзя — не следует. Даже ряд Тейлора, хоть и плохая идея, но, например, для функции exp(x) около нуля, если нужно быстро на коленках собрать хоть что-то, вполне себе рабочий вариант. А ещё лучше он подходит для стартовой точки в методе Ньютона, потому что благодаря простым коэффициентам легко переводится на арифметику с фиксированной запятой без лишней возни. Каждому методу — своя ситуация, где он хорош. Хотя в нашей общей практике, да, не применяется.
Я думаю, что вы имели ввиду не интерполяцию, а разложение в ряд Тейлора вокруг табличных значений.Не знаю, что имел в виду автор исходного комментария, но есть и другой способ вычисления синуса/косинуса с табличными значениями — через комплексное умножение, как это используется в быстром преобразовании Фурье. Например, зная значения для sin|cos(2) из таблицы, для вычисления sin|cos(2.1) достаточно посчитать sin|cos(0.1) через ряд Тейлора, а затем комплексно перемножить с sin|cos(2).
пример в Wolfram Mathematica
In[1]:= Cos[2.1] + I Sin[2.1]
Out[1]= -0.504846 + 0.863209 I
In[2]:= (Cos[2] + I Sin[2]) (Cos[0.1] + I Sin[0.1])
Out[2]= -0.504846 + 0.863209 I
In[3]:= (Cos[2] + I Sin[2]) (Cos[0.05] + I Sin[0.05]) (Cos[0.05] + I Sin[0.05])
Out[3]= -0.504846 + 0.863209 IРяд тейлора или полиномы Чебышева/минимакс-полиномом это всё равно полиномы. Принципиальной разницы быть не должно в самом коде. Если вы думаете по-другому — пишите.
Я думаю очень по-другому. Всё дело в том, что степени этих полиномов будут различными. Минимакс-полином может оказаться вдвое короче ряда Тейлора (в зависимости от выбранного отрезка). Например, для функции exp2(x) на отрезке [0,1] если нужна точность 54 бита, то хватает минимакс-полинома 10-й степени. А ряд Тейлора будет намного длиннее, 15-я степень. При этом его реализация будет намного сложнее. В результате моих экспериментов, Тейлор проиграл Минимаксу по скорости в 8 раз. Вы только вдумайтесь. Да, код можно ускорить, ну, раза в 3. А куда деть 15-ю степень? Поэтому разница всё-таки есть.
Не могли бы вы уточнить смысл столбцов таблицы и чисел в них?
Конечно, напишу немного больше. Спасибо за комментарий.
Пусть 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
Я буду очень благодарен, если вы напишите своё мнение о понятности этого комментария. Если всё ОК, то я добавлю это в исходную статью.
Пусть 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
Я буду очень благодарен, если вы напишите своё мнение о понятности этого комментария. Если всё ОК, то я добавлю это в исходную статью.
Спасибо за подробное объяснение
Не очень понятно, почему побитово сравнивается результат нашей и libc, а не sin(float128). И какой смысл при этом несут проценты, если при этом возможно что и при совпадающих побитово результатах они далеки от sin(float128)?
Спасибо за вопрос. Ответ достаточно простой. Потому для double(float64) просто недостижима точность quad (float128) потому, что в переменной double меньше значащих знаков. Всё, что мы можем сделать это попытаться как можно ближе подобраться к теоретическому пределу в 0.5ULP. И данное сравнение показывает кто всё-таки ближе. Написанная функция, или встроенная.
Благодарю вас за упоминания моих статей и за вашу работу. Действительно тема очень актуальна (в узких кругах) и хорошо, что ею интересуются. Будь у меня больше времени, я обязательно дал бы несколько конкретных советов о том, как было бы лучше сделать то, к чему вы пришли в этой статье, но этого времени нет, кроме того, вы обещали продолжение и возможно там будет как раз то, что нужно.
У меня есть желание поделиться ссылкой: существует книга Jean-Michel Muller, “Elementary Functions: Algorithms and Implementation”, 2016, в этой книге огромное количество знаний о вычислении элементарных функций. В принципе, её достаточно чтобы знать о них почти всё, что сейчас применяют в этой области, за исключением каких-то продвинутых алгоритмов, изобретённых недавно. В том числе и общую теорию. Например, там сказано, что не следует вычислять синус на отрезке от 0 до Pi/2, потому что это невыгодно, отрезок нужно уменьшать минимум до [0, Pi/256] и уже дальше особыми формулами распространять на весь период. Также есть целая глава, где сказано, что синус (и любую другую элементарную функцию на монотонном интервале) можно вычислить используя только сдвиги и сложения, а также небольшую таблицу предподсчитанных заранее констант (Shift-and-Add Algorithms).
Ещё рекомендую для ознакомления библиотеку libquadmath где почти все элементарные функции реализованы для float128 и весь код хорошо виден. Посмотрите как там реализован синус (файл sincosq_kernel.c), там сразу одновременно вычисляется синус и косинус. Метод не очень хитрый, но не очевидный: там и редукция аргумента, и minimax-полином (если я верно понял, это именно он), и финальная корректировка с помощью таблиц. В общем, я так понял, вы хотели рассказать об этом в следующей статье? Мне просто показалось, что это очень большой труд — всё это рассказать, и однажды это хотел сделать я, но прикинул, что нужно для этого ещё пару лет потренироваться на тех статьях, что я пишу сейчас. Поэтому элементарные функции решил пока оставить. Очень сложная для популярного изложения (для меня) тема. Само только вдумчивое чтение упомянутой книги далось не с одного месяца.
Я буду очень рад, если вы сможете сделать эту работу, на русском языке этого и правда не хватает катастрофически. Я согласен, что действительно нужно обязательно показать сначала, что разложением в ряд Тейлора считать не вполне хорошая идея, и есть методы получше, потому что большинство как раз об этом и не знает. Я догадываюсь, какой следующий метод вы опишите, но не буду делать спойлеры :)
PS. А ещё я начал бы с экспоненты или 2 в степени x, эти функции кажутся намного проще. Но это уже как вам угодно.
Спасибо, Артём. Я уважительно отношусь к вашему времени, поэтому буду краток. Действительно в следующей статье я подробно расскажу метод, который вы обрисовали сейчас. Вот ссылка на код, который будет обсуждаться в следующей статье libm sin.
Так же я Вам очень благодарен в том, что вы единственный, кто предоставил ссылку на код. Давайте же его посмотрим вместе. А именно на 126 строчку вот этого файла sinq_kernel.c. Я понимаю, что у вас мало времени, но, пожалуйста, потратьте ещё пару минут и посмотрите на него и одновременно посмотрите на 2-ую формулу моей статьи. Они одинаковы! Просто точь-в-точь!
Так же я Вам очень благодарен в том, что вы единственный, кто предоставил ссылку на код. Давайте же его посмотрим вместе. А именно на 126 строчку вот этого файла sinq_kernel.c. Я понимаю, что у вас мало времени, но, пожалуйста, потратьте ещё пару минут и посмотрите на него и одновременно посмотрите на 2-ую формулу моей статьи. Они одинаковы! Просто точь-в-точь!
Конечно, я тут спорить не готов в силу того, что в своей работе до тригонометрических функций ещё не добрался. Превосходно, если они одинаковые (хотя сами реализации и разбивка отрезка точно разная), но я вам в одном из комментариев выше привел пример, что ряд Тейлора может проигрывать minimax-полиному. Далее я всецело доверяюсь вашему мнению, потому что как уже сказано ранее, не чувствую в себе уверенность в том, чтобы спорить, когда сам ещё до сложной реализации элементарных тригонометрических функций не добрался. Повторюсь, что я только ЗА то, чтобы вы смогли завершить свой труд, потому что я в себе таких сил сейчас не вижу. А то я думал, что мне придётся это делать через год-два, а вот глядишь, может и не придётся :)
Спасибо Артём за ваш ответ В свою очередь я тоже, конечно, не спорю, что для других функций другие методы могут быть лучше. Рад что мы нашли общий язык. Сами понимаете, что не все полезут в код этой библиотеки по вашей ссылке и не все разберутся, что там на самом деле не минмакс полином. Поэтому мне нужно было отдельно это отметить. В остальном я прочитал ваши комментарии с большим интересом и принял к сведению. Замечание, что ряд Тейлора это не универсальная вещь, конечно, полезно.
Текущим моим комментарием я не хочу отвечать на негативные отзывы. Я хочу обратится к людям, которым статья понравилась, и, прочитав огромное количество комментариев (которые уже превышают размер статьи) могут усомниться в том, что я написал. Если вы начали сомневаться подумайте о следующих вещах.
1) «Такая точность не нужна» в различных вариациях. Безусловно она нужна не всегда. Но, если бы она никогда не была нужна, зачем вообще тогда придумали числа double? Зачем в С они по умолчанию (напр константа 1.0 это double, а 1.0f float (32-bit)). sin(x) это double, а sinf float (32-bit). Зачем nvidia добавляет в свои топовые карты поддержку double, если это ни кому не нужно?
2) «Ряд Тейлора безнадёжно устарел» в различных вариациях. Утверждение, конечно, сильное но совершенно не обоснованное. Никто (кроме ArtemKaravaev) не написал код лучше, либо не дал ссылку. А ссылка от Артёма приводит к полиному, который выглядит как ряд Тейлора и считается как ряд Тейлора. От того, что его назвать по-другому не изменится ничего.
3) «Статья ничего нового не рассказывает». Может быть, но если всё до нельзя банально, то почему никто из комментаторов не ответил на вопрос в конце? Почему никто не объяснил в двух словах почему у sin_e4 и sin_e5 результаты разные? И так, чтобы всем стало понятно.
Хочу сказать спасибо и обратиться к тем, кому понравилась эта статья. Если у Вас есть желание меня поддержать, 10 минут времени и базовые знания программирования, то я буду очень благодарен, если вы перейдёте по ссылке от Артёма (sinq_kernel.c) и напишите в комментариях является ли формула в строке 126 одинаковой с второй формулой в моей статье или нет.
1) «Такая точность не нужна» в различных вариациях. Безусловно она нужна не всегда. Но, если бы она никогда не была нужна, зачем вообще тогда придумали числа double? Зачем в С они по умолчанию (напр константа 1.0 это double, а 1.0f float (32-bit)). sin(x) это double, а sinf float (32-bit). Зачем nvidia добавляет в свои топовые карты поддержку double, если это ни кому не нужно?
2) «Ряд Тейлора безнадёжно устарел» в различных вариациях. Утверждение, конечно, сильное но совершенно не обоснованное. Никто (кроме ArtemKaravaev) не написал код лучше, либо не дал ссылку. А ссылка от Артёма приводит к полиному, который выглядит как ряд Тейлора и считается как ряд Тейлора. От того, что его назвать по-другому не изменится ничего.
3) «Статья ничего нового не рассказывает». Может быть, но если всё до нельзя банально, то почему никто из комментаторов не ответил на вопрос в конце? Почему никто не объяснил в двух словах почему у sin_e4 и sin_e5 результаты разные? И так, чтобы всем стало понятно.
Хочу сказать спасибо и обратиться к тем, кому понравилась эта статья. Если у Вас есть желание меня поддержать, 10 минут времени и базовые знания программирования, то я буду очень благодарен, если вы перейдёте по ссылке от Артёма (sinq_kernel.c) и напишите в комментариях является ли формула в строке 126 одинаковой с второй формулой в моей статье или нет.
Строка 126 называется "схема Горнера" и является стандартным способом вычисления полиномов (мне его даже в школе рассказывали). По сути, вы заявляете, что любой полином «выглядит как ряд Тейлора», что является в корне неверным.
Во-первых, вы путаете понятия математического разложения функции, и схемы её вычисления. Если вам так будет проще, то в строке 126 ряд Тейлора вычисляется по схеме Горнера. Абсолютно так же как мой пример sin_e4. Если сомневаетесь, то раскройте скобки в выражении и приведите получившийся полином к каноническому виду. Вы увидите, что это не что иное как запись разложения синуса в ряд Тейлора.
Во-вторых, я не заявляю, что любой полином выглядит как ряд Тейлора. Я говорю, что какие бы вы полиномы не использовали, Вам не удастся достичь вычислительной сложности (меньше сложений/умножений), лучше чем у ряда Тейлора при такой же «математической» (все операции проводятся абсолютно точно) точности при приведении конечного выражения к каноническому виду. Считаете по-другому — так приведите контрпример.
P.S. Есть один способ улучшения точности последнего члена разложения. О нём я поговорю позже.
Во-вторых, я не заявляю, что любой полином выглядит как ряд Тейлора. Я говорю, что какие бы вы полиномы не использовали, Вам не удастся достичь вычислительной сложности (меньше сложений/умножений), лучше чем у ряда Тейлора при такой же «математической» (все операции проводятся абсолютно точно) точности при приведении конечного выражения к каноническому виду. Считаете по-другому — так приведите контрпример.
P.S. Есть один способ улучшения точности последнего члена разложения. О нём я поговорю позже.
в строке 126 ряд Тейлора вычисляется по схеме ГорнераОчевидно, что нет. Видно, что, например, SIN2 отличается от 1/5! в 8-ми последних десятичных знаках (или 6-ти шестнадцатеричных) или SIN8 ≈ 2.81068755e-15 заметно отличается от 1/17! ≈ 2.81145725e-15. Любая аппроксимация синуса полиномом будет похожа на ряд Тейлора, особенно в первых членах, вот только в этих небольших отличиях и скрывается принципиальная разница в точности аппроксимации.
Я говорю, что какие бы вы полиномы не использовали, Вам не удастся достичь вычислительной сложности (меньше сложений/умножений), лучше чем у ряда Тейлора при такой же «математической» (все операции проводятся абсолютно точно) точности при приведении конечного выражения к каноническому виду.Это тоже не соответствует действительности. Членов ряда Тейлора для достижения требуемой точности надо гораздо больше. Про это уже выше ArtemKaravaev написал.
Начну со конца. Ну напишите, наконец-то полином Чебышева, который будет сходится быстрее, а расчётная сложность будет меньше. Уже написана программа для тестирования. Вам реально там остаётся добавить одну (!) функцию. Ну покажите нам всем это! Это же, я надеюсь, не секрет какой-то?
По первому вопросу DrSmile вы реально думаете, что читатели (или я) не понимают, что небольшая вариация коэффициентов в разложении может уменьшить общую ошибку? Это абсолютно банальное утверждение по вашему в статье поменяет что? Она поменяет схему оптимизации скорости? Она поменяет вывод что, сумму надо вести с мелких элементов? С другой стороны, если вы всё прекрасно знаете, то вы не можете не понимать, что количество и значение этих коэффициентов зависит от пределов вычисления синуса. Для [0, pi/2] одна, для [0, pi/256] совершенно другая. И что в самом начале это давать просто нерационально. Откройте код и посмотрите следующий метод, который я буду описывать.
Я поставил цель при написании этой статьи не самопиар, а планомерное объяснение шаблонов работы с числами с плавающей точкой и не только. Я правда призываю Вас не тратить ни свои ни мои силы на написание этих бесполезных комментариев, а взять и написать статью лучше. Получить много кармы от радостных подписчиков. Я, правда, считаю, что подбор этих коэффициентов выходит за рамки данного цикла. Как я написал в PS я только покажу подбор последнего члена. Это реально можно сделать аналитически красиво.
К сожалению я вижу в Ваших комментариях только желание принизить данную работу. У каждого, конечно, есть право на своё мнение и желания. У меня тоже есть право больше не отвечать на ваши комментарии, чем я и воспользуюсь.
По первому вопросу DrSmile вы реально думаете, что читатели (или я) не понимают, что небольшая вариация коэффициентов в разложении может уменьшить общую ошибку? Это абсолютно банальное утверждение по вашему в статье поменяет что? Она поменяет схему оптимизации скорости? Она поменяет вывод что, сумму надо вести с мелких элементов? С другой стороны, если вы всё прекрасно знаете, то вы не можете не понимать, что количество и значение этих коэффициентов зависит от пределов вычисления синуса. Для [0, pi/2] одна, для [0, pi/256] совершенно другая. И что в самом начале это давать просто нерационально. Откройте код и посмотрите следующий метод, который я буду описывать.
Я поставил цель при написании этой статьи не самопиар, а планомерное объяснение шаблонов работы с числами с плавающей точкой и не только. Я правда призываю Вас не тратить ни свои ни мои силы на написание этих бесполезных комментариев, а взять и написать статью лучше. Получить много кармы от радостных подписчиков. Я, правда, считаю, что подбор этих коэффициентов выходит за рамки данного цикла. Как я написал в PS я только покажу подбор последнего члена. Это реально можно сделать аналитически красиво.
К сожалению я вижу в Ваших комментариях только желание принизить данную работу. У каждого, конечно, есть право на своё мнение и желания. У меня тоже есть право больше не отвечать на ваши комментарии, чем я и воспользуюсь.
Точные… вычисления ...
double sin...Если честно, ждал более точного, нежели double. Например чего-нибудь, где требуемая точность задается.
Жаль Вас разочаровывать, но таких исследований в данном цикле статей не будет. Впрочем, метод, который я буду обсуждать в следующей статье может быть адаптирован для любой заданной точности.
Да метод и этой статьи можно адаптировать, если правильно вычислить n (не 25) и не для IEEE double. Точность достижима любая.
Просто вопрос: если в прикладной задаче хватает точности IEEE double, то целесообразно ли переписывать реализацию стандартных функций?
Просто вопрос: если в прикладной задаче хватает точности IEEE double, то целесообразно ли переписывать реализацию стандартных функций?
Ответ на Ваш вопрос не нельзя назвать простым. Ни одна функция в стандартной библиотеке не гарантирует абсолютной точности 0.5ULP. Это, как абсолютный ноль, величина не достижимая. Ну по крайней мере для сложных функций. Поэтому точность функции это всегда компромисс между её скоростью и самой точностью. Главный мой интерес, как программиста — разобраться, как устроена библиотека, повторить, а может быть в чём-то превзойти. Например, при такой же точности сделать быстрее.
Я могу ошибаться, но разве когда вы делаете что-то вроде double\float x = sin(double\float any value); исполняется не ассемблерная инструкция? Тогда же просто нереально сделать быстрее при сохранении точности.
При вызове функции из библиотеки, она, конечно, исполняется в машинных кодах. Но этот год генерируется из С кода s_sin.c, так и из ассемблера svml_d_sin4_core_avx2.S. Безусловно это обсудим в следующих статьях.
Спасибо за эту ссылку, реализацию в машинном коде мне реально было интересно увидеть. Извините, если ещё один вопрос покажется глупым: вы не знаете, почему в реализации синуса библиотеки libc для x86_64 так много кода (т.е. они «вручную» нормируют аргумент sin(arg) = sin(N*Pi + R) = (-1)^N * sin® и вычисляют sin® is approximated by corresponding polynomial), разве в x86_64 нет просто инструкции типа sin/fsin?
Спасибо за замечательный вопрос. Конечно инструкция fsin есть и для аргумента double, но
1) У неё есть одна неприятная особенность, связанная с точностью intel fsin. Вызывать её просто так нельзя (см. график и «Conclusion» в статье). Те, кто считают CORDIS новым, современным методом, прочитайте, пожалуйста первый параграф документа по ссылке.
2) Она медленнее. На моём железе (i5 mobile) ~ в 4 раза. Через некоторое время я отредактирую исходную статью (полный текст программы) и добавлю туда ассемблерный вызов синуса. Вы сможете у себя проверит разницу в производительности.
1) У неё есть одна неприятная особенность, связанная с точностью intel fsin. Вызывать её просто так нельзя (см. график и «Conclusion» в статье). Те, кто считают CORDIS новым, современным методом, прочитайте, пожалуйста первый параграф документа по ссылке.
2) Она медленнее. На моём железе (i5 mobile) ~ в 4 раза. Через некоторое время я отредактирую исходную статью (полный текст программы) и добавлю туда ассемблерный вызов синуса. Вы сможете у себя проверит разницу в производительности.
Т.е. простыми словами, FSIN архитектуры x86_64 (2000 год) такой же «точный», как и в реализации >10 лет назад (до 1990) и для практического применения куда лучше (читать: обязательно) использовать код libc, т.к. он точнее и даже внезапно быстрее.
Простыми словами, разработчики glibc очень умные люди и не стали бы они делать плохо.
Чуть более сложно: Функция fsin некоторых случаях даже точнее за счёт использования расширенных регистров. Но в некоторых случаях она даёт неприемлемый результат, и в добавок ко всему, сильно медленнее. Я очень подозреваю, что данная функция оставлена только для совместимости с ранними программами. Косвенно на это указано во втором параграфе документа.
P.S. Статья обновлена. Можно сравнивать.
Чуть более сложно: Функция fsin некоторых случаях даже точнее за счёт использования расширенных регистров. Но в некоторых случаях она даёт неприемлемый результат, и в добавок ко всему, сильно медленнее. Я очень подозреваю, что данная функция оставлена только для совместимости с ранними программами. Косвенно на это указано во втором параграфе документа.
P.S. Статья обновлена. Можно сравнивать.
А можете объяснить, почему оптимизированная версия программы быстрее аж в 10000 раз?
Вот это постоянно не пересчитывается
pow(x, i) / bm::factorial(i)
Если вкратце, то возведение в произвольную степень очень затратная функция. Вот её код. e_pow.c Просто посмотрите на сколько он сложнее, чем x * x.
По поводу точности sin_e4, там для малых x суммирование с 1 будет убивать точность, а в sin_e5 происходит суммирование более-менее сопоставимых величин.
Ход мыслей у Вас очень правильный. Но утверждение 2: «происходит суммирование более-менее сопоставимых величин.» не верно. Величины 1.0 и y, абсолютно так же несапоставимы как и x и x * y, если x не 0 разумеется. Попробуете докрутить?
Могу еще предположить, что при преобразовании выражения x+x*y в x*(1-y) нам надо прозвести деление, что-то типа x*(1 + x*y/x), а в вычислениях с плавающей точкой это не всегда можно выполнить без потери точности. Ну и это начинает оказывать сильное влияние, т.к. потеря точности происходит на последней итерации, т.е. в самом значимом слагаемом ряда.
Да. Вы правы. То, что это происходит на последней итерации это реально важно. Тут проще будет понять на примере. Представьте, что у Вас десятичное плавающее число с 4-мя знаками точности. Например x*(1+y). y=1.234E-2, x=5.678E-2, 1=1.000E0. Попробуйте провести операции с данными числами, не забывая округлять.
Вопрос, где же ссылки на стандартную теорию численных методов? Это ничто иное, как решение простых задач из данной дисциплины.
Это практическая работа. Научности тут не будет. И не потому, что я не умею так писать. Умею. Я учёный, у меня больше 20 работ в иностранных журналах. Этим циклом я бы хотел дать читателям возможность руками пощупать все подводные камни. Если вам кажется, что задачи слишком простые, может в следующих статьях вы найдёте для себя что-то более сложное.
Возможно, это очень глупый вопрос, но я не нашёл на него ответа: почему числа в стобцах «лучше наша» и «лучше libm» не дают в сумме 100%?
Простите, не смог пройти мимо. Избавился от sign.
double sin_e2(double x) {
double result = 0;
double xx = -x * x;
double pw = x;
double fti = 1.0;
for(int i = 1; i < 25; i += 2) {
fti /= i;
result += pw * fti;
fti /= ( i + 1 );
pw *= xx;
}
return result;
}
array<double, 25> fc;
…
double res = fc[25];
Точно не double res = fc[24];?
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Точные и быстрые вычисления для чисел с плавающей точкой на примере функции синуса. Введение и часть 1