05.11.2020 03:10
dorex 26 1700 Источник
Продолжаю цикл статей по работе с плавающей точкой. В первой статье я дал небольшое математическое введение и показал самый простой и очевидный способ вычисления синуса с примерами программ с разными «подводными камнями». Сегодня статья будет немного другая по стилю. Здесь не будет практики, зато мы копнём глубже математику и залезем в святая-святых — код стандартной библиотеки. Так же я дам ответ на вопрос в конце первой статьи. Итак, поехали.

Снова немного математики


Очевидно, чем быстрее убывает по модулю ряд Тейлора, тем меньше слагаемых нужно, чтобы добиться необходимой точности. И тем, кажется, точнее будет результат (ниже это будет обсуждено подробнее). Для сравнения, например, возьмём член седьмой степени ряда Тейлора ($\frac{x^7}{7!}$) при $x=1.0$ и $x=0.1$. Значения выражения будут $1.198\cdot 10^{-4}$ и $1.198\cdot 10^{-11}$ соответственно. Огромная разница, не правда ли? Так давайте попробуем найти способ уменьшить верхний предел интервала расчёта функции синуса.

Разложение в ряд вокруг заданных значений


Для того, чтобы понять этот способ нам нужно вернуться на первый курс института и вспомнить определения ряда Тейлора (wiki). В двух словах: зная функцию и её производные в какой-то точке можно найти значения функции в окрестностях этой точки разложением в ряд Тейлора. Для функции синуса это означает следующее

$\sin(x_0+\Delta x) \approx \sin(x_0)+\sin'(x_0)\frac{\Delta x}{1!}+\sin''(x_0)\frac{\Delta x^2}{2!}+\sin'''(x_0)\frac{\Delta x^3}{3!}+\cdots$


Что даёт нам этот подход с практической точки зрения? Представьте себе, что у нас есть интервал от $0$ до $\pi/2$. Выберем на этом интервале 10 линейно распределённых точек (выбор не оптимальный): $x_0=0$, $x_1=\pi/20$, $x_2=2\pi/20$, $x_i=i\cdot \pi/20$. Для каждой точки посчитаем табличку с синусом и его производными в данной точке. Теперь можно модифицировать функцию так, что при получении значения $x$ функция берёт ближайшее значение $x_i$ и раскладывает в ряд вокруг точки $x_i$, а не вокруг нуля ($\Delta x = x - x_i$).

Использование тригонометрических преобразований


Если вернуться ещё дальше, в старшие классы школы, то можно вспомнить одну очень важную формулу:

$\sin(x_0+\Delta x)=\sin(x_0)\cos(\Delta x)+\cos(x_0)\sin(\Delta x) $


А дальше всё так же как и в предыдущем параграфе. Выбираем точки внутри интервала, считаем для них синус и косинус, а при вызове функции синуса ищем ближайшую и, используя формулу выше, вычисляем синус используя малое значение $\Delta x$.
Подумайте, какой из этих двух способов лучше выбрать, а пока мы перейдём от математики к практическим вычислениям.

Распределительное свойство умножения в мире чисел с плавающей точкой


Пришлось спросить совета у интернета, как же всё-таки называется $a\cdot(b+c)=a\cdot b + a\cdot c$. Оказывается распределительное свойство. Вернёмся к вопросу, который я задал в конце первой части. А именно, почему математически эквивалентные выражения $a\cdot(b+c)$ и $a\cdot b + a\cdot c$ могут давать разные результаты в вычислениях с плавающей точкой? Проще всего это проиллюстрировать на примере. Возьмём гипотетическую систему, которая работает с числами с плавающей точкой в десятичном формате с 4-мя знаками точности. Предположим, что $b=1.000\text{E}0$, $c=1.234\text{E-}2$, а $a=5.678\text{E-}2$. Для начала возьмём выражение со скобками и посчитаем его шаг за шагом, не забывая округлять на каждом шаге:
1) $b+c=1.000\text{E}0 + 1.234\text{E-}2=1.012\text{E}0$
2) $a\cdot(b+c)=5.678\text{E-}2\times 1.012\text{E}0=5.746\text{E-}2$
Полученный ответ $y=5.746\text{E-}2$
Теперь давайте так же шаг за шагом вычислим второе выражение:
$a\cdot b + a\cdot c=5.678\text{E-}2\times 1.000\text{E}0 + 5.678\text{E-}2\times 1.234\text{E-}2=5.678\text{E-}2 + 7.007\text{E-}4=5.748\text{E-}2$
Полученный ответ $y=5.748\text{E-}2$
Истинный ответ 0.0574806652.

Как видите, полученный ответ во втором случае гораздо ближе к истинному, чем в первом. Если объяснять это на пальцах, то представьте себе, что когда в первом случае мы добавляем к 1.0 число $c=1.234\text{E-}2$ мы просто выкидываем последние две цифры. Их больше нет. Во-втором же случае, выкидывание происходит в самом конце, уже после умножения. Т.е. во втором случае операция(-ии) умножения точнее.

Вроде на этом можно и закончить, но присмотритесь ещё раз к первому способу и скажите какой будет результат вычисления $b+c-b$. И… У нас есть способ округления чисел с плавающей точкой! Не проходите мимо этого примера. Дайте себе время в нём разобраться. Округление чисел будет нами очень интенсивно использоваться далее в этой и следующих статьях.

Подметим ещё одну особенность данного выражения. Представьте себе, что нам точности в 4 цифры у переменной недостаточно. Что делать? И тут ответ у нас уже есть — представить число в виде $b+c$ и хранить её в памяти в виде суммы двух цифр. И, соответственно, проводить операции (например умножения) отдельно для обоих слагаемых. Этот приём более подробно описан в статье Сложение двух чисел с плавающей запятой без потери точности.

В предыдущей статье я так же написал, что у способа $a\cdot(b+c)$ есть одна неприятная особенность. И она заключается вот в чём. Число $c$ всегда отрезается по последней значащей цифре числа $b$. Это значит, что в независимости от числа $c$, если $b+c \neq b$, то ошибка в последнем знаке возможна всегда даже при малых $c$. А это недопустимо в подходе в следующей главе.

Как это работает на примере библиотеки GNU


Ну как? Выбрали какой из двух способов описанных в начале статьи Вы выбрали для точного расчёта синуса? Какой бы способ вы не выбрали — оба способа правильные. Более того — абсолютно идентичные. Хоть поверьте, хоть проверьте. Ниже я буду использовать школьные формулы. Они проще для объяснения.

Вооружившись знаниями полученными в предыдущей и в этой статье можно легко понять код стандартной библиотеки. Откроем файл s_sin.c и найдём там функцию __sin:
image
Код её достаточно простой. Легко понять, что она вызывает разный набор функций в зависимости от пределов входной переменной. В этой статье обсудим участок кода 218-224 для углов 2^-26<|x|< 0.855469. Видно что в данном участке кода вызывается функция do_sin (x, 0). На данной функции мы остановимся подробнее:

image

  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 ) и косинуса ( c ) от 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, а ssn и ccs — малые.
  6. 144-145. Используя вторую формулу из этой статьи получим sin(u+x)=(sn+ssn)*(1-c)+(cs+ccs)*s. Уже зная, как правильно складывать и умножать числа с плавающей точкой, раскройте скобки и полученное выражение сравните с 144-145. Кто проведёт вычисления — получит в конце небольшой сюрприз.

На самом деле я описал лишь самую простую часть вычисления синуса данным способом. За «бортом» осталось очень много математики. Например, как вычислить размер таблицы и сами элементы в ней? Откуда взялись магические числа 0.126 и 0.855469? Когда обрубать вычисление рядом Тейлора? Поправки в коэффициенты ряда Тейлора, для уточнения результата.

Всё это, конечно, интересно, но, объективно, представленный способ имеет много недостатков: необходимо вычисление синуса ( s ) и косинуса ( c ) одновременно, что требует в два раза больше вычислений ряда Тейлора1. Умножение на табличные значения, как мы видим, тоже не бесплатно. Так же хранение таблицы в 3520 байт в оперативной памяти, конечно, не проблема, а вот обращаться к ней (даже в кеше) может быть накладно.

Поэтому в следующей части мы попытаемся избавится от таблички и посчитать синус в интервале [0.126, 0.855469] напрямую, но более точно, чем в первой главе.

Прежде чем закончить — вопрос на сообразительность. Число big в данном примере 52776558133248=3*244. Откуда взялось такое число, на не, например, 245? Сформулирую вопрос более точно. Почему при округлении чисел оптимально число 3*2N, а не, например, 2N+1? Другой вопрос, какое N нужно выбрать, чтобы округлить число до целого?

1Стоит отметить, что существенное преимущество у такого подхода может появится при одновременном вычислении синуса и косинуса от одного и того же угла. Вторую функцию можно вычислить почти бесплатно.