Функция греха Apple libm
Этот вопрос касается реализации функции sin в математической библиотеке Apple libm C, доступной здесь . Основная функция sin начинается в строке 736. Короче говоря: когда x настолько мало, что ближайшее к sin(x) число с плавающей запятой равно x, почему функция выполняет умножение с плавающей запятой, результат которого никогда не используется, а затем возвращает x вместо того, чтобы просто немедленно возвращать x?
Функция использует разные методы в зависимости от величины xk входных данных x. Меня интересует случай, когда xk настолько мал, что ближайшее к sin(x) число с плавающей запятой составляет всего x, строки 809-829. Я бы предположил, что в этом случае функция просто выполнила бы return x. Вот что он делает вместо этого:
/* Make t0 volatile to force compiler to fetch it at run-time
rather than optimizing away the multiplication.
*/
static volatile const double t0 = 1/3.;
/* Get t0 once. If we wrote "t0*t0", the compiler would load it
twice, since it is volatile.
*/
const double t1 = t0;
/* Perform a multiplication and pass it into an assembly construct to
prevent the compiler from knowing we do not use the result and
optimizing away the multiplication.
*/
__asm__("" : : "X" (t1*t1));
// Return the floating-point number nearest sine(x), which is x.
return x;
Таким образом, он создает новый статический изменчивый const double t0, равный 1/3., другой const double t1 с тем же значением, умножает t1 на себя так, как он утверждает, что компилятор не может оптимизировать... затем возвращает x. Почему он делает это, а не просто возвращает x в начало?
@Ctx Я подумал об этом, но не думаю, что это правильно. Помимо всего остального, это просто математическая библиотека общего назначения.
Ответы 1
… почему функция выполняет умножение с плавающей запятой, результат которого никогда не используется, а затем возвращает
x…
Операция умножения с плавающей запятой имеет два выхода. Одним из них является числовое значение в регистре с плавающей запятой, а другим — обновленные флаги в битах состояния с плавающей запятой. В том месте кода, о котором вы спрашиваете, x настолько мало, что ближайшее к sin x представимое значение равно x, поэтому оно будет возвращаемым значением функции. Однако синус x не совсем x, поэтому мы хотели бы поднять флаг неточного. Умножение 1/3. само по себе достигает этого.
Вы заметите аналогичную заботу в return x * (1 - 0x1p-53); в строке 803. В этот момент известно, что x находится в субнормальном диапазоне (здесь написано «ненормально», но более уместно «субнормально»), поэтому x * (1 - 0x1p-53) имеет, из-за округления, то же самое. значение как x. Но это также поднимает флаг неточной точности, если x не равно нулю. (А синус нуля точен, поэтому лучше не поднимать флаг неточной точности, когда x равен нулю.)
Я думаю, ты понял, спасибо. Они включают в себя заявление об отказе от ответственности «может или не может быть точным, даже если результат точен», но, похоже, они прилагают некоторые усилия, чтобы сделать это правильно в двух самых маленьких случаях.
@MatthewTowers: «Они» — это я. Я написал эту программу.
черт возьми. Мы очень ценим ваши обширные комментарии в источнике.
@EricPostpischil Вызвал бы внутренний эквивалент fesetexcept(FE_INEXACT) то же самое или это более тонко?
@IanAbbott: Я думаю, что все по-другому; fesetexcept только поднимает флаг. Выполнение фактической операции с плавающей запятой с неточным результатом не только поднимает флаг, но и вызывает исключение, если оно включено для неточных результатов. В любом случае, простая загрузка подготовленного значения и инструкции умножения с плавающей запятой — это практически самый быстрый способ поднять флаг.
Я подозреваю, что sin(x) неточно для всех значений x радиан, кроме x == +/- 0. Теперь меня интересует попытка правильно установить биты состояния с плавающей запятой для sine_grades(x), поскольку это дает точные рациональные результаты при 0 +180°*n, +/-30 +180°*n. Раньше я выступал за уменьшение угла в восьмых долях (360°/8), затем вызов поддиапазона sine_grades_45(x) и последующее построение ответа. Я вижу, что для достижения правильных битов состояния с плавающей запятой коду может потребоваться сначала уменьшить угол на двенадцатые доли (360 °/12) (или, может быть, на шестые). Спасибо за понимание построения функций.
@chux-ReinstateMonica: Да, sin не имеет рациональных значений аргументов (представимых значений), кроме нуля. Таким образом, это будет неточно, за исключением нуля, NaN и ±∞. (Это вызывает недопустимость сигнализации NaN или ±∞). Утверждение, что это может или не может привести к неточным результатам, даже если результат точен, по сути, является частью спецификации шаблона для большинства подпрограмм математической библиотеки. Это более применимо, например, в log10, где log10(1e6) может выполнить некоторые вычисления и дать точный результат 6, но промежуточные вычисления будут неточными.
@chux-ReinstateMonica: Поскольку sin всегда вызывает неточные результаты для конечных аргументов, отличных от нуля, от этого нет особого смысла. Любой, кто звонит sin, обычно должен ожидать неточного ответа. Однако гипотетически это полезно для тех, кто пытается выполнить точные вычисления и вызывает библиотечную процедуру, которая, по его мнению, может быть точной, но вызывает внутренний вызов sin. Как правило, стандартные подпрограммы математической библиотеки рассматривались как вещи, которые мы хотели разумно уважать семантику с плавающей запятой. Для людей, которым нужна скорость, есть библиотеки Vector/SIMD, которые могут игнорировать флаги и тому подобное.
Дикая догадка: возможно, для предотвращения атак по побочным каналам времени, когда функция используется в контекстах, связанных с безопасностью?