Вы когда-нибудь задумывались, как компьютеры вычисляют такие тригонометрические функции, как sin и cos? Возможно, вы помните из исчисления, что любая дифференцируемая функция может быть выражена в виде бесконечной полиномиальной суммы, известной как ряд Тейлора. Я всегда предполагал, что реализации тригонометрических функций используют ряды Тейлора под капотом, потому что процессоры способны только к базовой арифметике.

Недавно я решил погрузиться в исходный код glibc (библиотека GNU C), чтобы проверить, действительно ли функция sin была реализована с использованием ряда Тейлора. Читая исходный код для расчета ряда Тейлора, я заметил некоторые небольшие улучшения, которые я мог бы сделать. В этой статье обсуждаются изменения, которые я внес в glibc, и базовые знания, необходимые для их понимания.

Обзор серии Тейлора

Ряд Тейлора позволяет аппроксимировать любую дифференцируемую функцию полиномом. Эта полиномиальная аппроксимация становится все более точной по мере того, как вы включаете в нее больше членов. Каждый ряд Тейлора также имеет центральную точку, и аппроксимация становится менее точной для значений, удаленных от центра. Все приближения Тейлора на рисунке ниже сосредоточены вокруг 0.

Ниже представлено разложение в ряд Тейлора для общей функции f(x) вокруг центральной точки b, которое можно доказать с помощью интегрирования по частям.

Мы можем использовать рисунок 1, чтобы получить следующее расширение для функции sin с центром вокруг 0.

glibc Синусоидальная реализация

Если вы перейдете к файлу s_sin.c до моего изменения, вы увидите комментарий о том, что код реализует следующий вариант расширения sin Тейлора.

Это именно то, что мы имеем на рис. 2, за исключением загадочного d_aterm в конце. Я попытаюсь объяснить, что такое d_a и откуда взялся этот последний термин.

Синусоидальная реализация glibc интерпретирует свой ввод x как x = a + d_a, где d_a — небольшое число порядка 10^-19. Этот метод выражения одного 64-битного значения в виде суммы двух 64-битных значений известен как двойная двойная арифметика. Это распространенный метод эмуляции более высокой точности, когда у вас есть только 64-битные примитивы для работы. Реализация требует более высокой точности для борьбы с ошибками округления, которые возникают при каждой арифметической операции ряда Тейлора.

Давайте проверим последний член для себя, подставив a + d_a в формулу на рисунке 2. Нам нужно только включить d_a в первые 2 члена разложения, потому что при более высоких мощностях эффект d_a становится незначительным. Поскольку мы вычисляем эффект d_a алгебраически, а затем включаем его в конец, это гарантирует, что d_a не будет стерто последующими арифметическими операциями с плавающей запятой.

После расширения через Wolfram Alpha и упрощения получаем следующее соотношение.

Мы можем опустить два последних члена в круглых скобках, потому что они являются более высокими степенями крошечного значения d_a, поэтому член в круглых скобках на рисунке 5 упрощается до d_a — a^2 * d_a / 2. Это отличается от термина d_a * (1 — a^2)/2, включенного в комментарий. Если вы посмотрите на исходный код, вы увидите, что макрос TAYLOR_SIN фактически вычисляет нашу версию d_a — a^2 * d_a / 2 при оценке ряда Тейлора. Так что видимо комментарий некорректен.

Также странно, что приведенный ниже код макроса TAYLOR_SIN предполагает, что a всегда равно x, потому что он вычисляет ³ путем умножения xx (что равно x²) и a. Вы можете идентифицировать термин a³, найдя, какой термин имеет старший коэффициент s1, который является предварительно вычисленным значением для -1/3!. Вы также увидите, что a действительно равно x, если вы проверите единственный сайт вызова макроса TAYLOR_SIN.

Теперь, когда вы видели и поняли весь соответствующий код, вот краткое изложение 2 изменений, которые я объединил в glibc.

  1. Обновите комментарий, чтобы включить правильный последний термин d_a — a^2 * d_a / 2
  2. Уточните, что макрос TAYLOR_SIN на самом деле ожидает x в качестве второго параметра, переименовав a → x (и аналогично da → dx). Также обновите формулу в документации макроса, чтобы отразить это.

Вы можете возразить, что мой вклад был бессмысленным, потому что ни одно из этих улучшений фактически не приводит к изменению поведения в библиотеке. Но, надеюсь, вы узнали что-то, прочитав этот пост, и нашли исходный код glibc немного более доступным.

Заключительные мысли

После попытки разобраться в математической библиотеке glibc я вновь проникся уважением к ее создателям и сопровождающим. Я не знал, что для точного вычисления основных математических функций требуется так много продвинутых численных вычислительных методов. Я до сих пор удивлен, что любопытство в конце концов привело меня к небольшому улучшению широко используемой базовой системной библиотеки. Работать над этим было очень весело, потому что это требовало объединения концепций исчисления и компьютерной архитектуры. Я надеюсь сделать больше вкладов с открытым исходным кодом в будущем.

Я также хотел бы поблагодарить Siddhesh Poyarekar и Paul Zimmermann за ответы на все мои вопросы по glibc и проверку моего кода.

Найдите меня в Твиттере и загляните на мой личный сайт.

Первоначально опубликовано на https://www.awelm.com 18 декабря 2021 г.