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

Question

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

Ожидаемая цель функции предназначена для графики, поэтому двойная точность не нужна; Однако на функциях целевой платформы, Trig (и SQRT), которые ведут поплавки (SINF, COSF, ATAN2F, ASINF, ACOSF и SQRTF) работают быстрее на удвоении, чем на плавании (вероятно, потому что инструкция для расчета таких значений может фактически требовать двойной; если пропущено поплавок, значение должно быть отброшено в двойную, что требует копирования его в область с большим количеством памяти - то есть накладным расходом). В результате все переменные, участвующие в функции, являются двойной точностью.

Вот проблема: я пытаюсь оптимизировать мою функцию, чтобы его можно было назвать более раз в секунду. Поэтому я заменил звонки на SIN, COS, SQRT, ET CETERA с вызовами в версиях плавающих точек этих функций, поскольку они приводят к увеличению скорости в 3-4 разах в целом. Это работает практически для всех входов; Однако, если входные векторы близки к параллельному стандартным блокам векторов (i, j или k), ошибках округлая для различных функций, которые накапливаются достаточно, чтобы вызвать более поздние вызовы функций SQRTF или обратного трига (ASINF, ASINF, ASOSF, atan2f) пройти аргументы, которые едва за пределами области этих функций.

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

Итак, любые идеи?

Редактировать: Кто-то выразил путаницу над моим использованием двойных операций с плавающей запятой. Функция намного быстрее, если я фактически хранишь все мои значения в двойных емковых контейнерах (т.е. переменные двойных типов), чем если я храним их в контейнерах по плаванию размера. Тем не менее, операции Trig Trig Floating Point находятся быстрее, чем у двойных точных тригристов для очевидных причин.

Answer 1

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

Answer 2

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

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

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 Scalar Math, где это возможно.

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

// 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 Scalar Math включен. Вы можете сделать это с флагом компилятора /arch:sse2 на MSVC или -mfpmath=sse на GCC.

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

Answer 3

Вы смотрели на Графическая программирование Черная книга Или, возможно, передавать расчеты к вашему грастру?