Как я могу учесть ошибки округления в арифметике с плавающей точкой для обратных функций trig (и sqrt) (в C)?

Question

У меня есть довольно сложная функция, которая принимает несколько двойных значений, которые представляют два вектора в 3-х пространствах формы (величина, широта, долгота), где широта и долгота указаны в радианах, и угол.Цель функции - повернуть первый вектор вокруг второго на заданный угол и вернуть результирующий вектор.Я уже проверил, что код логически корректен и работает.

Ожидаемое назначение функции для графики, поэтому двойная точность не требуется;однако на целевой платформе функции trig (и sqrt), которые принимают числа с плавающей запятой (в частности, sinf, cosf, atan2f, asinf, acosf и sqrtf), работают быстрее с двойными числами, чем с числами с плавающей запятой (вероятно, потому что инструкция для вычисления таких значений может фактически потребоватьdouble; если передано значение с плавающей запятой, значение должно быть приведено к двойному значению, что требует копирования его в область с большим объемом памяти - т. е. накладные расходы).В результате все переменные, включенные в функцию, имеют двойную точность.

Вот проблема: я пытаюсь оптимизировать свою функцию так, чтобы ее можно было вызывать больше раз в секунду.Поэтому я заменил вызовы sin, cos, sqrt и так далее вызовами версий этих функций с плавающей запятой, поскольку они приводят к общему увеличению скорости в 3-4 раза.Это работает почти для всех входов;однако, если входные векторы близки к параллельным со стандартными единичными векторами (i, j или k), ошибок округления для различных функций накапливается достаточно, чтобы вызвать последующие вызовы функций sqrtf или обратного триггера (asinf, acosf,atan2f) передать аргументы, которые едва находятся вне области действия этих функций.

Итак, я остаюсь с этой дилеммой: либо я могу вызывать только функции двойной точности и избежать проблемы(и в конечном итоге с ограничением около 1 300 000 векторных операций в секунду), или я могу попытаться придумать что-то еще.В конечном счете, я хотел бы найти способ санировать входные данные для функций обратного триггера, чтобы позаботиться о крайних случаях (это тривиально для sqrt: просто используйте abs).Ветвление не вариант, так как даже один условный оператор добавляет столько накладных расходов, что любой прирост производительности теряется.

Итак, есть какие-нибудь идеи?

Редактировать: кто-то выразил недоумение по поводу того, что я использую удвоение по сравнению соперации с плавающей запятой.Функция намного быстрее, если я на самом деле храню все свои значения в контейнерах двойного размера (переменные двойного типа IE), чем если бы я сохранял их в контейнерах плавающего размера.Однако триггерные операции с плавающей запятой выполняются быстрее, чем триггерные операции с двойной точностью по очевидным причинам.

Answer 1

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

Answer 2

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

Первое достаточно легко при ветвлении, например

bool IsAlmostEqual( float a, float b ) { return fabs(a-b) < 0.001f; } // or
bool IsAlmostEqual( float a, float b ) { return fabs(a-b) < (a * 0.0001f); } // for relative error

но это грязно. Зажимать входные данные домена немного сложнее, но лучше. Ключ должен использовать операторы условного перемещения , которые обычно делают что-то вроде

float ExampleOfConditionalMoveIntrinsic( float comparand, float a, float b ) 
{ return comparand >= 0.0f ? a : b ; }

в одном операторе без ветвления.

Они различаются в зависимости от архитектуры. На модуле с плавающей запятой x87 вы можете сделать это с помощью операции условного перемещения FCMOV , но это неуклюже, потому что это зависит от ранее установленных флагов условий, поэтому это медленно. Кроме того, нет единого встроенного компилятора для cmov. Это одна из причин, почему мы избегаем x87 с плавающей точкой в ​​пользу скалярной математики SSE2, где это возможно.

Условное перемещение намного лучше поддерживается в SSE путем объединения оператора сравнения с побитовым AND. Это предпочтительно даже для скалярной математики:

// assuming you've already used _mm_load_ss to load your floats onto registers 
__m128 fsel( __m128 comparand, __m128 a, __m128 b ) 
{
    __m128 zero = {0,0,0,0};
    // set low word of mask to all 1s if comparand > 0
    __m128 mask = _mm_cmpgt_ss( comparand, zero );  
    a = _mm_and_ss( a, mask );    // a = a & mask 
    b = _mm_andnot_ss( mask, b ); // b = ~mask & b
    return _mm_or_ss( a, b );     // return a | b
    }
}

Компиляторы лучше, но не очень хороши в генерации такого рода шаблонов для троичных файлов, когда включена скалярная математика SSE2. Вы можете сделать это с флагом компилятора /arch:sse2 на MSVC или -mfpmath=sse на GCC.

На PowerPC и многих других архитектурах RISC fsel() - это аппаратный код операции и, следовательно, обычно также встроенный компилятор.

Answer 3

Вы смотрели на "Черную книгу" по программированию графики или, возможно, передали вычисления на свой графический процессор?