¿Cómo se puede explicar que para los errores de redondeo en aritmética de punto flotante para trig inversa (y raíz cuadrada) en funciones (C)?

Question

Tengo una función bastante complicado que tiene varios valores dobles que representan dos vectores en 3-espacio de la forma (magnitud, latitud, longitud) en la latitud y longitud son en radianes, y un ángulo. El propósito de la función es para hacer girar el primer vector alrededor de la segunda por el ángulo especificado y devolver el vector resultante. Ya he verificado que el código es lógicamente correcto y funciona.

El objetivo esperado de la función es para los gráficos, por lo doble precisión no es necesaria; Sin embargo, en la plataforma de destino, trig (y sqrt) funciones que toman flotadores (senf, COSF, atan2f, asinf, acosf y sqrtf específicamente) el trabajo más rápido en dobles que en los flotadores (probablemente porque la instrucción para calcular estos valores en realidad puede requerir una dobles; si se pasa un flotador, el valor debe ser fundido a una doble, lo que requiere de copiarlo en un área con más memoria - es decir, encima de la cabeza). Como resultado, todas las variables que intervienen en la función son de doble precisión.

Aquí está el problema: Estoy tratando de optimizar mi función para que se le puede llamar más veces por segundo. Por tanto, he reemplazado las llamadas a sin, cos, sqrt, etcétera con llamadas a las versiones de punto flotante de esas funciones, tal como resultan en un aumento de velocidad de 3-4 veces en total. Esto funciona para casi todas las entradas; Sin embargo, si los vectores de entrada están cerca de paralelas con los vectores unitarios estándar (i, j, o K), de redondeo errores para las diversas funciones de acumular suficiente para causar las llamadas posteriores a sqrtf o inversa funciones trigonométricas (asinf, acosf, atan2f) para pasar argumentos que son apenas fuera del dominio de esas funciones.

Por lo tanto, me quedo con este dilema: o bien que sólo puedo llamar funciones de doble precisión y evitar el problema (y terminar con un límite de alrededor de 1.300.000 operaciones vectoriales por segundo), o puede tratar de llegar a algo más . En última instancia, me gustaría una manera de desinfectar la entrada a la inversa trig funciones para cuidar de casos extremos (es trivial para hacerlo por sqrt: uso justo abs). La ramificación no es una opción, ya que incluso una sola sentencia condicional añade mucha sobrecarga que cualquier beneficio de rendimiento se pierden.

Por lo tanto, cualquier idea?

Editar: alguien expresado confusión sobre mis utilizando dobles frente a operaciones de punto flotante. La función es mucho más rápido si realmente almacenar todos mis valores en recipientes de tamaño doble (variables de tipo POR EJEMPLO el doble) que si los almacena en recipientes de flotación de tamaño. Sin embargo, las operaciones flotantes punto de precisión trigonométricas son más rápidas que las operaciones de doble precisión trig por razones obvias.

Answer 1

Básicamente, es necesario encontrar un algoritmo href="http://en.wikipedia.org/wiki/Numerical_stability" rel="nofollow"> estable numéricamente número de condición si los pasos individuales. Y de hecho, puede ser imposible si el problema subyacente es en sí mal condicionada.

Answer 2

coma flotante de precisión simple introduce inherentemente error. Por lo tanto, es necesario construir su matemáticas de manera que todas las comparaciones tienen un cierto grado de "decantación" mediante el uso de un factor épsilon, y que necesita para higienizar las entradas a las funciones con dominios limitados.

El anterior es lo suficientemente fácil cuando de ramificación, por ejemplo

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

pero que de desordenado. Apriete entradas de dominio es un poco más difícil, pero mejor. La clave es usar condicional operadores movimiento , que en hacer algo general, como

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

en una sola op, sin incurrir en una rama.

Estos varían dependiendo de la arquitectura. En la unidad de coma flotante x87 puede hacerlo con el FCMOV condicional movimiento op , pero que es torpe, porque depende de indicadores de condición de ser establecido previamente, por lo que es lento. Además, no es un compilador intrínseca consistente para cmov. Esta es una de las razones por las que evitamos x87 punto flotante a favor de matemáticas SSE2 escalar cuando sea posible.

movimiento condicional es mucho mejor soportado en SSE por el emparejamiento de un operador de comparación con un AND bit a bit. Esto es incluso preferible para las matemáticas escalar:

// 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
    }
}

Los compiladores son mejores, pero no excelente, sobre emisión de este tipo de patrón para ternarios cuando se habilita SSE2 matemáticas escalar. Usted puede hacer eso con la bandera /arch:sse2 compilador MSVC o en -mfpmath=sse en GCC.

En el PowerPC y muchas otras arquitecturas RISC, fsel() es un código de operación del hardware y por lo tanto por lo general un compilador intrínseca también.

Answer 3

¿Has mirado en la gráfico de programación Libro Negro o tal vez de entregar los cálculos fuera de la GPU?