Cómo hacer el cálculo de hypot2 (x, y) cuando los números pueden desbordar

Question

Me gustaría hacer un cálculo de hypot2 en un procesador de 16 bits.

La fórmula estándar es c = sqrt((a * a) + (b * b)). El problema con esto es con grandes entradas que se desborda. P.ej. 200 y 250, multiplique 200 * 200 para obtener 90,000, que es mayor que el valor máximo de 32,767 firmado, por lo que se desborda, como hace b, los números se agregan y el resultado puede ser inútil; incluso podría indicar una condición de error debido a un sqrt negativo.

En mi caso, estoy tratando con números de 32 bits, pero la multiplicación de 32 bits en mi procesador es muy rápida, alrededor de 4 ciclos. Estoy usando un microcontrolador dsPIC. Prefiero no tener que multiplicar con números de 64 bits porque eso está desperdiciando memoria preciosa e indudablemente será más lento. Además, solo tengo sqrt para números de 32 bits, por lo que los números de 64 bits requerirían otra función. Entonces, ¿cómo puedo calcular una hipótesis cuando los valores pueden ser grandes?

Tenga en cuenta que realmente solo puedo usar números enteros para esto. Usar algo como matemáticas de coma flotante implica un golpe de velocidad que prefiero evitar. Mi procesador tiene una rutina rápida integer/punto fijo atan2, alrededor de 130 ciclos; ¿podría usar esto para calcular la longitud de la hipotenusa?

Answer 1

Según la precisión que necesite, puede evitar los cuadrados y la operación de raíz cuadrada. Hay una sección sobre este tema en Understanding Digital Signal Processing by Rick Lyons (sección 10.2, "Aproximación de Magnitud de Vector de Alta Velocidad", comenzando en la página 400 en mi edición).

la aproximación es esencialmente:

magnitude = alpha * min + beta * max 

donde max y min son los valores absolutos mínimos de los componentes reales e imaginarios máximo y, y alfa y beta son dos constantes que se eligen para dar una distribución de error razonable sobre el rango de interés.Estas constantes se pueden representar como fracciones con potencia de 2 divisores para mantener el aritmético simple/eficiente. En el libro que sugiere alpha = 15/16, beta = 15/32, y a continuación, puede simplificar la fórmula a:

magnitude = (15/16) * (max + min/2) 

que podría ser implementado como sigue usando operaciones de enteros:

magnitude = 15 * (max + min/2)/16 

y, por supuesto, podemos usar los cambios para las divisiones:

magnitude = (15 * (max + (min >> 1))) >> 4 

El error es de +/- 5% en un cuadrante.

Más información sobre esta técnica aquí: http://www.dspguru.com/dsp/tricks/magnitude-estimator

Answer 2

¿Necesita la precisión completa? Si no lo hace, puede aumentar un poco su rango descartando algunos bits menos significativos y multiplicándolos después.

Can a y b ¿nada? ¿Qué tal una tabla de búsqueda si solo tiene unos pocos a y b que debe calcular?

Answer 3

Una solución sencilla para evitar el desbordamiento es dividir tanto a y b por a+b antes de la cuadratura, y después multiplicar la raíz cuadrada por a+b. O haga lo mismo con max(a,b).

Answer 4

Puede hacer un poco de álgebra simple para devolver los resultados al rango.

sqrt((a * a) + (b * b)) 
= 2 * sqrt(((a * a) + (b * b))/4) 
= 2 * sqrt((a * a)/4 + (b * b)/4) 
= 2 * sqrt((a/2 * a/2) + (b/2 * b/2)) 

Answer 5

Aniko y John, me parece que no ha tratado el problema del OP. Si ayb son números enteros, es probable que un * a + b * b se desborde, porque se están realizando operaciones enteras. La solución obvia es convertir ayb en valores de coma flotante antes de calcular a * a + b * b. Pero el OP no nos ha permitido saber qué idioma deberíamos usar, así que estamos un poco bloqueados.

Answer 6

Como esencialmente no puedes hacer ninguna multiplicación sin desbordamiento, es probable que pierdas algo de precisión.

para obtener los números en un rango aceptable, sacar algún factor x y utilizar

c = x*sqrt((a/x)*(a/x) + (b/x)*(b/x)) 

Si x es un factor común, que no perderá precisión, pero si no es así, perderá precisión .

Actualización: Aún mejor, considerando que se puede hacer un trabajo suave con números de 64 bits, con una sola adición de 64 bits, se puede hacer el resto de este problema en 32-bits con sólo una pequeña pérdida de precisión. Para hacer esto: haga las dos multiplicaciones de 32 bits para darle dos números de 64 bits, agréguelos y luego cambie los bits según sea necesario para volver a reducir la suma a 32 bits antes de tomar la raíz cuadrada. Si siempre mueve un bit en 2 bits, simplemente multiplique el resultado final por 2^(la mitad del número de cambios de bit), según la regla anterior.El truncamiento solo debería causar una pérdida de precisión muy pequeña, no más de 2^31 o 0.00000005% de error.

Answer 7

esto es tomado textualmente de this @John D. Cook blog post, por lo tanto, CW:

Aquí es cómo calcular sqrt(x*x + y*y) sin correr el riesgo de desbordamiento.

1. max = maximum(|x|, |y|)
2. min = minimum(|x|, |y|)
3. r = min/max
4. return max*sqrt(1 + r*r)

Si @ John D. Cook llega y post con los cuales se le debe dar al aceptar :)

Answer 8

La fórmula estándar es c = sqrt ((a * a) + (b * b)). El problema con esto es con entradas grandes que se desborda.

La solución para los desbordamientos (aparte de arrojar un error) es saturar sus cálculos intermedios.

Calcular C = a * a + b * b. Si a y b tienen números de 16 bits firmados, nunca tendrá un desbordamiento. Si son números sin signo, primero deberá desplazar hacia la derecha las entradas para que la suma encaje en un número de 32 bits.

Si C> (MAX_RADIUS)^2, regrese MAX_RADIUS, donde MAX_RADIUS es el valor máximo que puede tolerar antes de contar un desbordamiento.

De lo contrario, utilice sqrt() o CORDIC algorithm, que evita el costo de las raíces cuadradas en favor de iteración de bucle + suma + desplazamiento, para recuperar la amplitud del vector (a, b).

Answer 9

Si puede restringir que a y b tengan como máximo 7 bits, no obtendrá ningún desbordamiento. Puedes usar una instrucción de contar cero para determinar cuántos bits debes tirar.

Supongamos que a> = b.

int bits = 16 - count_leading_zeros(a); 
if (bits > 7) { 
    a >>= bits - 7; 
    b >>= bits - 7; 
} 
c = sqrt(a*a + b*b); 
if (bits > 7) { 
    c <<= bits - 7; 
} 

Las porciones de los procesadores tienen esta instrucción hoy en día, y si no, se pueden utilizar otros fast techniques.

Aunque esto no le dará la respuesta exacta, será muy cerca (a lo sumo ~ 1% bajo).