Как компьютер вычисляет значение трансцендентных функций?

[alex18]

Добрый день! Помогите, пожалуйста, ответить на вопрос по теме машинного вычисления значения трансцендентных функций $\sin(x)$, $e^x$, $\ln(x)$. Преподаватель озвучил следующую проблему: если аргумент очень велик, то приближение рядом Тейлора сопряжено с большими погрешностями (тем более, у $\sin(x)$ члены разложения знакопеременные, а если где-то происходит вычитание близких чисел, то за этим следует большая потеря точности), что делать? Также нужно дать комментарии по идее в случае с функцией $\sin(x)$ воспользоваться периодичностью и вычесть из аргумента необходимое количество периодов, чтобы сдвинуть его, насколько возможно, к нулю, но, если аргумент был огромен, то, учитывая погрешность округления числа $\pi$, вычитание из аргумента большого количества $\pi$ также приведёт к потере точности - так как и когда можно достаточно точно вычислить $\sin(x)$?

У меня есть несколько общих идей: перегруппировать сложение членов ряда, чтобы минимизировать погрешность (складывать близкие числа и избегать вычитания близких чисел), для $e^x$ домножить и поделить $x$ на $10^{-k}$, чтобы $x \cdot 10^{-k}$ было достаточно мало (и тогда экспонента с помощью приближения первыми членами ряда Тейлора будет вычисляться достаточно точно), и результат возвести в степень $10^k$ (что тоже сопряжено с погрешностями, выиграем ли точность таким образом?), а для задачи с логарифмом в случае большого $x$ перебором найти такое натуральное $k$, чтобы $x^\frac{1}{k}$ было достаточно близко к единице, для него посчитать приближённое значение логарифма c помощью нескольких первых членов ряда Тейлора в единице, и потом умножить результат на $k$.
Для $\sin(x)$ и идеи использовать периодичность можно, например, заранее вычислить "сетку" из некоторых значений, кратных $\pi$, (например, каждые $1000 \pi n$) с машинной точностью (для этого специально выделив при вычислении дополнительную точность), и тогда погрешность при сдвиге на периоды удастся сделать поменьше.

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

[wrest]

По поводу синусов, мне кажется что $\pi$ известно как предвычисленное с такой огромной точностью, что вычесть из аргумента синуса $2 \pi n$ можно без ущерба для любых практически требуемых $n$

[mihaild]

Это огромная тема, и на практике для разных функций приходится думать. Но ряд Тейлора работает плохо всегда. Общая схема - как-то заменить аргумент значением поменьше и использовать какую-нибудь хорошую полиномиальную аппроксимацию. Для экспоненты, например, можно написать $e^x = e^a \cdot 2^b$ где $b$ целое (и на $2^b$ соответственно можно умножать без потери точности), а $|a| \leqslant \frac{\ln 2}{2}$. Для логарифма иногда можно написать $x = \frac{1 - y}{1 + y}$ и вместо $\ln x$ считать $\ln(1 \pm y)$.
Для точного вычисления синуса в любом случае нужно точно знать далекие разряды аргументов. Если аргумент в 32-битном числе с плавающей запятой и по модулю больше чем $2^{26}$, то синус от него может быть любым числом от $-1$ до $1$, и с этим ничего не поделать.

[worm2]

О, я как раз писал такую библиотеку.
С синусом так и не решил проблему (я сам предвычислял $\pi$, готовыми результатами не пользовался). У меня была мысль воспользоваться формулами двойного угла: $(\sin\alpha, \cos\alpha) \to (\sin 2\alpha, \cos 2\alpha)$, но я её не довёл до логического завершения. Поскольку в моей библиотеке максимальное представимое число было $10^{10^8}$, то аргумент функции $e^x$ был ограничен сравнительно небольшим значением. С остальными функциями ещё проще, тут уже написали.

[ewert]

mihaild в сообщении #1531606 писал(а):

Для логарифма иногда можно написать $x = \frac{1 - y}{1 + y}$ и вместо $\ln x$ считать $\ln(1 \pm y)$.

Только непонятно, зачем. Число в машине хранится как комбинация двоичного порядка и мантиссы, причём мантисса нормализована.
Благодаря нормализации двоичный логарифм мантиссы эффективно считается, скажем, чебышёвским приближением, после чего надо просто извлечь порядок из соответствующего поля и прибавить.

[Null]

(Оффтоп)

Чую программирование 1-2 курса мехмата МГУ.
Если мне память не изменяет, придется писать длинную арифметику и делить на $2\pi$ столбиком. Типа $1000000000$ вполне себе точное число и синус должен считаться с любой требуемой точностью доступной в типе double(C/C++). Преподу надо сказать что погрешность результата не ниже погрешности исходного числа в типе double. То есть если $\varepsilon$ - машинный эпсилон, то посчитать $\sin(x)$ с точностью меньше $\max(\varepsilon x,\varepsilon)$ не получиться(не требуется). Если $x>1/\varepsilon$ выдаем ноль, типа точнее не возможно.

[Sender]

Null в сообщении #1531634 писал(а):

делить на $2\pi$ столбиком

Боюсь, если требуется считать быстро, столбик - не самый эффективный метод. Обычно деление заменяется умножением, а методы быстрого умножения известны.

[Null]

Sender в сообщении #1531639 писал(а):

Боюсь, если требуется считать быстро, столбик - не самый эффективный метод. Обычно деление заменяется умножением, а методы быстрого умножения известны.

Разрядов там мало для быстрых способов. Но да умножить можно.

[Евгений Машеров]

Тейлор для вычисления функций плох. Если только аргумент $\Delta x$ не слишком мал. Да и тогда можно найти получше. Обычно находят способ приведения аргумента в определённый интервал, и ищут наилучшее в этом интервале приближение. Естественным способом приведения для тригонометрических функций является вычитание $ 2\pi n$. Погрешность задания числа Пи здесь, конечно, может сказаться, но только при очень больших n: даже если использовать тип данных float (одинарная точность), существенная ошибка появится, учитывая, что для этого типа в среднем 7.2 десятичных знака, при n порядка сотен тысяч, для наиболее распространённого double при n порядка сотни триллионов (15.9 знака) и для extended (long double), также аппаратно поддерживаемого, это уже должно быть что-то около сотни миллионов миллиардов (19.2 знака). Но если n настолько велико - надо пересматривать алгоритм. Что-то в нём не так. Очевидно приведение для гамма-функции, для экспоненты и логарифма.
Помимо представления в виде полиномов (Чебышева, в частности), могут использоваться дробно-рациональные функции, есть методы вычисления "цифра за цифрой", особенно хорошие для тригонометрических функций.
Могу указать на справочник Вычисление функций на ЭВМ. Попов Б.А., Теслер Г.С. Наукова думка 1984
и более ранний Благовещенского и Теслера.

[Someone]

Ю. Люк. Специальные математические функции и их аппроксимации. "Мир", Москва, 1980.

[Евгений Машеров]

Вот ещё можно
Muller J.-M. etc. Handbook of Floating-Point Arithmetic

[Евгений Машеров]

Ещё есть
Иванов В.В. Методы вычислений на ЭВМ
К.: Наукова думка, 1986.
Но там вообще численные методы, вычисление функций лишь в одной главе книги.
И по поводу Тейлора. Он вещь не бесполезная, но напрямую для вычисления функций плох. Он "узкий специалист согласно Бернарду Шоу", который "знает всё ни о чём". Даёт наилучшее приближение не только для значения функции, но и всех её производных в точке. А нам надо в отдалении от точки. Можно, например, взяв ряд Тейлора с очень большим числом членов, непрактичным для расчёта, но гарантирующим нужную точность в выбранном интервале, затем его "телескопировать", представив через Чебышева, отбросив члены с совсем уж малыми коэффициентами, и вернувшись к обычным полиномам, но теперь уже равномерное приближение в отрезке.
Или вообще не полиномы, а, скажем, CORDIC. Очень хорош для как раз тригонометрии (гиперболические тоже).
А возвращаясь к "синусу от очень большого аргумента" - вычитание $2\pi$ и вообще приведение замечательно помогает, а значения аргумента, при которых доступная точность числа Пи недостаточна, это либо чистая казуистика, в реальности не бываемая, либо грандиозный прокол в разработке алгоритма. Хотя тут надо посоветовать использовать заложенное в библиотеках значение этой константы с полутора десятками знаков, довольствоваться привычным 3.14 или даже "Это я знаю и помню прекрасно" недостаточно, даже продолжение "пи многие знаки мне лишни-напрасны" может не хватить.

[ewert]

(Оффтоп)

Евгений Машеров в сообщении #1531646 писал(а):

extended (long double), также аппаратно поддерживаемого, это уже должно быть что-то около сотни миллионов миллиардов (19.2 знака)

Ньюанец. Не знаю, что такое long double, но extended -- это стандартный десятибайтовый формат внутреннего представления в регистрах процессора (в отличие от выноса double наружу), и хранит он всего лишь 17-18 десятичных цифр.

Да, раз уж пошла такая пьянка. Естественно, универсально оптимальных алгоритмов действительно не бывает. (далее оффтопик)

Вот, скажем, классическая задачка -- как извлечь квадратный корень.

Есть очень древний алгоритм; не помню, какого года и месяца издания -- но по существу это просто (более поздний формально, видимо) метод Ньютона: $x_{k+1}=x_k-\frac1{2x_k}$. И конкретно для этой задачки он оптимален.

Вопрос: а с какого начального приближения начинать?...
Сходится-то он с любого начального приближения; но какое оптимально?...

Почему-то мне только вчера пришла в голову банальная мысль. Надо просто разделить пополам двоичный порядок, мантиссу же усреднить с единицей.

А ещё лучше не усреднить, а просто заменить на единицу, т.к. деление порядка пополам, реализуемое двоичным сдвигом, этот порядок занижает.

Можно ли это дальше оптимизировать? -- сильно вряд ли. Ну допустим, улучшим мы начальное приближение. Максимум, на что можно в данном случае рассчитывать -- это на уменьшение количества итераций примерно на единицу. Т.е. примерно на три арифметические операции. Любая попытка оптимизации займёт заведомо больше.

Да и замена квадратичной скорости сходимости на, допустим, кубическую ситуацию тоже качественно не меняет.

Ну это так, лирика была; прошу пардону.

[Евгений Машеров]

long double - употребляемое в некоторых компиляторах название 10-байтного плавающего. Из 80 бит на мантиссу приходится у него 65, что и позволяет заявлять о "19.2 десятичных разрядах". 17-18 это уже "разумная перестраховка".