Продолжаю цикл статей по работе с плавающей точкой. В первой статье я дал небольшое математическое введение и показал самый простой и очевидный способ вычисления синуса с примерами программ с разными «подводными камнями». Сегодня статья будет немного другая по стилю. Здесь не будет практики, зато мы копнём глубже математику и залезем в святая-святых — код стандартной библиотеки. Так же я дам ответ на вопрос в конце первой статьи. Итак, поехали.
Очевидно, чем быстрее убывает по модулю ряд Тейлора, тем меньше слагаемых нужно, чтобы добиться необходимой точности. И тем, кажется, точнее будет результат (ниже это будет обсуждено подробнее). Для сравнения, например, возьмём член седьмой степени ряда Тейлора () при и . Значения выражения будут и соответственно. Огромная разница, не правда ли? Так давайте попробуем найти способ уменьшить верхний предел интервала расчёта функции синуса.
Для того, чтобы понять этот способ нам нужно вернуться на первый курс института и вспомнить определения ряда Тейлора (wiki). В двух словах: зная функцию и её производные в какой-то точке можно найти значения функции в окрестностях этой точки разложением в ряд Тейлора. Для функции синуса это означает следующее
Что даёт нам этот подход с практической точки зрения? Представьте себе, что у нас есть интервал от до . Выберем на этом интервале 10 линейно распределённых точек (выбор не оптимальный): , , , . Для каждой точки посчитаем табличку с синусом и его производными в данной точке. Теперь можно модифицировать функцию так, что при получении значения функция берёт ближайшее значение и раскладывает в ряд вокруг точки , а не вокруг нуля ().
Если вернуться ещё дальше, в старшие классы школы, то можно вспомнить одну очень важную формулу:
А дальше всё так же как и в предыдущем параграфе. Выбираем точки внутри интервала, считаем для них синус и косинус, а при вызове функции синуса ищем ближайшую и, используя формулу выше, вычисляем синус используя малое значение .
Подумайте, какой из этих двух способов лучше выбрать, а пока мы перейдём от математики к практическим вычислениям.
Пришлось спросить совета у интернета, как же всё-таки называется . Оказывается распределительное свойство. Вернёмся к вопросу, который я задал в конце первой части. А именно, почему математически эквивалентные выражения и могут давать разные результаты в вычислениях с плавающей точкой? Проще всего это проиллюстрировать на примере. Возьмём гипотетическую систему, которая работает с числами с плавающей точкой в десятичном формате с 4-мя знаками точности. Предположим, что , , а . Для начала возьмём выражение со скобками и посчитаем его шаг за шагом, не забывая округлять на каждом шаге:
1)
2)
Полученный ответ
Теперь давайте так же шаг за шагом вычислим второе выражение:
Полученный ответ
Истинный ответ 0.0574806652.
Как видите, полученный ответ во втором случае гораздо ближе к истинному, чем в первом. Если объяснять это на пальцах, то представьте себе, что когда в первом случае мы добавляем к 1.0 число мы просто выкидываем последние две цифры. Их больше нет. Во-втором же случае, выкидывание происходит в самом конце, уже после умножения. Т.е. во втором случае операция(-ии) умножения точнее.
Вроде на этом можно и закончить, но присмотритесь ещё раз к первому способу и скажите какой будет результат вычисления . И… У нас есть способ округления чисел с плавающей точкой! Не проходите мимо этого примера. Дайте себе время в нём разобраться. Округление чисел будет нами очень интенсивно использоваться далее в этой и следующих статьях.
Подметим ещё одну особенность данного выражения. Представьте себе, что нам точности в 4 цифры у переменной недостаточно. Что делать? И тут ответ у нас уже есть — представить число в виде и хранить её в памяти в виде суммы двух цифр. И, соответственно, проводить операции (например умножения) отдельно для обоих слагаемых. Этот приём более подробно описан в статье Сложение двух чисел с плавающей запятой без потери точности.
В предыдущей статье я так же написал, что у способа есть одна неприятная особенность. И она заключается вот в чём. Число всегда отрезается по последней значащей цифре числа . Это значит, что в независимости от числа , если , то ошибка в последнем знаке возможна всегда даже при малых . А это недопустимо в подходе в следующей главе.
Ну как? Выбрали какой из двух способов описанных в начале статьи Вы выбрали для точного расчёта синуса? Тут могу можно поздравить всех. Оба способа правильные. Более того — абсолютно идентичные. Хоть поверьте, хоть проверьте. Ниже я буду использовать школьные формулы. Они проще для объяснения.
Вооружившись знаниями полученными в предыдущей и в этой статье можно легко понять код стандартной библиотеки. Откроем файл s_sin.c и найдём там функцию __sin:
Код её достаточно простой. Легко понять, что она вызывает разный набор функций в зависимости от пределов входной переменной. В этой статье обсудим участок кода 218-224 для углов 2^-26<|x|< 0.855469. Видно что в данном участке кода вызывается функция do_sin (x, 0). На данной функции мы остановимся подробнее:
1) Примем, что dx=0 для упрощения объяснений.
2) 129-130 Сразу вначале понимаем, что когда abs(x)<0.126, т.е. x величина достаточно малая, вычисляем синус напрямую рядом Тейлора. Поскольку в данном макросе нет ничего существенно нового, по сравнению с кодом, описанным в предыдущей части, то идём сразу дальше.
3) 136-137. Округление числа, о котором мы говорили выше. По сути x раскладывается на 2 части. Большая часть u и остаток x. Для гипотетического примера, предположим что у нас есть число 0.345678. После данной операции данное число разложится на два u=0.34, а х станет 0.005678.
4) 140-142. Вычисление синуса (s) и косинуса © от x из предыдущего пункта рядом Тейлора. Прошу заметить, что cos(x)=1-c, потому что в разложении нет первого члена 1.0, а остальные члены идут с противоположными знаками (см. исходный файл), чем в ряду Тейлора для косинуса.
5) 143. Извлекаются табличные значения для переменной u. Если использовать наш гипотетический пример сверху, при u=0.34 берётся элемент таблицы под номером 34. sin(u)=sn+ssn, cos(u)=cs+ccs. sn и cs — «большие» части значения синуса и косинуса в точке u, а snn и ccs — малые.
6) 144-145. Используя вторую формулу из этой статьи получим sin(u+x)=(sn+ssn)*(1-c)+(cs+ccs)*s. Уже зная, как правильно складывать и умножать числа с плавающей точкой, раскройте скобки и полученное выражение сравните с 144-145. Кто проведёт вычисления — получит в конце небольшой сюрприз.)))
На самом деле я описал лишь самую простую часть вычисления синуса данным способом. За «бортом» осталось очень много математики. Например, как вычислить размер таблицы и сами элементы в ней? Откуда взялись магические числа 0.126 и 0.855469? Когда обрубать вычисление рядом Тейлора? Поправки в коэффициенты ряда Тейлора, для уточнения результата.
Всё это, конечно, интересно, но, объективно, представленный способ имеет много недостатков: необходимо вычисление синуса (s) и косинуса © одновременно, что требует в два раза больше вычислений ряда Тейлора1. Умножение на табличные значения, как мы видим, тоже не бесплатно. Так же хранение таблицы в 3520 байт в оперативной памяти, конечно, не проблема, а вот обращаться к ней (даже в кеше) может быть накладно.
Поэтому в следующей части мы попытаемся избавится от таблички и посчитать синус в интервале [0.126, 0.855469] напрямую, но более точно, чем в первой главе.
Прежде чем закончить — вопрос на сообразительность. Число big в данном примере 52776558133248=3*2^44. Откуда взялось такое число, на не, например, 2^45? Сформулирую вопрос более точно. Почему при округлении чисел оптимально число 3*2^N, а не, например, 2^(N+1)? Другой вопрос, какое N нужно выбрать, чтобы округлить число до целого?
1Стоит отметить, что существенное преимущество у такого подхода может появится при одновременном вычислении синуса и косинуса от одного и того же угла. Вторую функцию можно вычислить почти бесплатно.
Снова немного математики
Очевидно, чем быстрее убывает по модулю ряд Тейлора, тем меньше слагаемых нужно, чтобы добиться необходимой точности. И тем, кажется, точнее будет результат (ниже это будет обсуждено подробнее). Для сравнения, например, возьмём член седьмой степени ряда Тейлора () при и . Значения выражения будут и соответственно. Огромная разница, не правда ли? Так давайте попробуем найти способ уменьшить верхний предел интервала расчёта функции синуса.
Разложение в ряд вокруг заданных значений
Для того, чтобы понять этот способ нам нужно вернуться на первый курс института и вспомнить определения ряда Тейлора (wiki). В двух словах: зная функцию и её производные в какой-то точке можно найти значения функции в окрестностях этой точки разложением в ряд Тейлора. Для функции синуса это означает следующее
Что даёт нам этот подход с практической точки зрения? Представьте себе, что у нас есть интервал от до . Выберем на этом интервале 10 линейно распределённых точек (выбор не оптимальный): , , , . Для каждой точки посчитаем табличку с синусом и его производными в данной точке. Теперь можно модифицировать функцию так, что при получении значения функция берёт ближайшее значение и раскладывает в ряд вокруг точки , а не вокруг нуля ().
Использование тригонометрических преобразований
Если вернуться ещё дальше, в старшие классы школы, то можно вспомнить одну очень важную формулу:
А дальше всё так же как и в предыдущем параграфе. Выбираем точки внутри интервала, считаем для них синус и косинус, а при вызове функции синуса ищем ближайшую и, используя формулу выше, вычисляем синус используя малое значение .
Подумайте, какой из этих двух способов лучше выбрать, а пока мы перейдём от математики к практическим вычислениям.
Распределительное свойство умножения в мире чисел с плавающей точкой
Пришлось спросить совета у интернета, как же всё-таки называется . Оказывается распределительное свойство. Вернёмся к вопросу, который я задал в конце первой части. А именно, почему математически эквивалентные выражения и могут давать разные результаты в вычислениях с плавающей точкой? Проще всего это проиллюстрировать на примере. Возьмём гипотетическую систему, которая работает с числами с плавающей точкой в десятичном формате с 4-мя знаками точности. Предположим, что , , а . Для начала возьмём выражение со скобками и посчитаем его шаг за шагом, не забывая округлять на каждом шаге:
1)
2)
Полученный ответ
Теперь давайте так же шаг за шагом вычислим второе выражение:
Полученный ответ
Истинный ответ 0.0574806652.
Как видите, полученный ответ во втором случае гораздо ближе к истинному, чем в первом. Если объяснять это на пальцах, то представьте себе, что когда в первом случае мы добавляем к 1.0 число мы просто выкидываем последние две цифры. Их больше нет. Во-втором же случае, выкидывание происходит в самом конце, уже после умножения. Т.е. во втором случае операция(-ии) умножения точнее.
Вроде на этом можно и закончить, но присмотритесь ещё раз к первому способу и скажите какой будет результат вычисления . И… У нас есть способ округления чисел с плавающей точкой! Не проходите мимо этого примера. Дайте себе время в нём разобраться. Округление чисел будет нами очень интенсивно использоваться далее в этой и следующих статьях.
Подметим ещё одну особенность данного выражения. Представьте себе, что нам точности в 4 цифры у переменной недостаточно. Что делать? И тут ответ у нас уже есть — представить число в виде и хранить её в памяти в виде суммы двух цифр. И, соответственно, проводить операции (например умножения) отдельно для обоих слагаемых. Этот приём более подробно описан в статье Сложение двух чисел с плавающей запятой без потери точности.
В предыдущей статье я так же написал, что у способа есть одна неприятная особенность. И она заключается вот в чём. Число всегда отрезается по последней значащей цифре числа . Это значит, что в независимости от числа , если , то ошибка в последнем знаке возможна всегда даже при малых . А это недопустимо в подходе в следующей главе.
Как это работает на примере библиотеки GNU
Ну как? Выбрали какой из двух способов описанных в начале статьи Вы выбрали для точного расчёта синуса? Тут могу можно поздравить всех. Оба способа правильные. Более того — абсолютно идентичные. Хоть поверьте, хоть проверьте. Ниже я буду использовать школьные формулы. Они проще для объяснения.
Вооружившись знаниями полученными в предыдущей и в этой статье можно легко понять код стандартной библиотеки. Откроем файл s_sin.c и найдём там функцию __sin:
Код её достаточно простой. Легко понять, что она вызывает разный набор функций в зависимости от пределов входной переменной. В этой статье обсудим участок кода 218-224 для углов 2^-26<|x|< 0.855469. Видно что в данном участке кода вызывается функция do_sin (x, 0). На данной функции мы остановимся подробнее:
1) Примем, что dx=0 для упрощения объяснений.
2) 129-130 Сразу вначале понимаем, что когда abs(x)<0.126, т.е. x величина достаточно малая, вычисляем синус напрямую рядом Тейлора. Поскольку в данном макросе нет ничего существенно нового, по сравнению с кодом, описанным в предыдущей части, то идём сразу дальше.
3) 136-137. Округление числа, о котором мы говорили выше. По сути x раскладывается на 2 части. Большая часть u и остаток x. Для гипотетического примера, предположим что у нас есть число 0.345678. После данной операции данное число разложится на два u=0.34, а х станет 0.005678.
4) 140-142. Вычисление синуса (s) и косинуса © от x из предыдущего пункта рядом Тейлора. Прошу заметить, что cos(x)=1-c, потому что в разложении нет первого члена 1.0, а остальные члены идут с противоположными знаками (см. исходный файл), чем в ряду Тейлора для косинуса.
5) 143. Извлекаются табличные значения для переменной u. Если использовать наш гипотетический пример сверху, при u=0.34 берётся элемент таблицы под номером 34. sin(u)=sn+ssn, cos(u)=cs+ccs. sn и cs — «большие» части значения синуса и косинуса в точке u, а snn и ccs — малые.
6) 144-145. Используя вторую формулу из этой статьи получим sin(u+x)=(sn+ssn)*(1-c)+(cs+ccs)*s. Уже зная, как правильно складывать и умножать числа с плавающей точкой, раскройте скобки и полученное выражение сравните с 144-145. Кто проведёт вычисления — получит в конце небольшой сюрприз.)))
На самом деле я описал лишь самую простую часть вычисления синуса данным способом. За «бортом» осталось очень много математики. Например, как вычислить размер таблицы и сами элементы в ней? Откуда взялись магические числа 0.126 и 0.855469? Когда обрубать вычисление рядом Тейлора? Поправки в коэффициенты ряда Тейлора, для уточнения результата.
Всё это, конечно, интересно, но, объективно, представленный способ имеет много недостатков: необходимо вычисление синуса (s) и косинуса © одновременно, что требует в два раза больше вычислений ряда Тейлора1. Умножение на табличные значения, как мы видим, тоже не бесплатно. Так же хранение таблицы в 3520 байт в оперативной памяти, конечно, не проблема, а вот обращаться к ней (даже в кеше) может быть накладно.
Поэтому в следующей части мы попытаемся избавится от таблички и посчитать синус в интервале [0.126, 0.855469] напрямую, но более точно, чем в первой главе.
Прежде чем закончить — вопрос на сообразительность. Число big в данном примере 52776558133248=3*2^44. Откуда взялось такое число, на не, например, 2^45? Сформулирую вопрос более точно. Почему при округлении чисел оптимально число 3*2^N, а не, например, 2^(N+1)? Другой вопрос, какое N нужно выбрать, чтобы округлить число до целого?
1Стоит отметить, что существенное преимущество у такого подхода может появится при одновременном вычислении синуса и косинуса от одного и того же угла. Вторую функцию можно вычислить почти бесплатно.
