c - teorema - hipotenusa formula

¿Algoritmo rápido de hipotenusa para procesador integrado?

Question

¿Existe un algoritmo inteligente / eficiente para determinar la hipotenusa de un ángulo (es decir, sqrt(a² + b²) ), al usar la matemática de punto fijo en un procesador integrado sin multiplicar el hardware?

Answer 1

A menos que esté haciendo esto a> 1 kHz, multiplicar incluso en una MCU sin hardware MUL no es terrible. Lo que es mucho peor es el sqrt . Intentaría modificar mi aplicación para que no tenga que calcularla en absoluto.

Las bibliotecas estándar probablemente serían las mejores si realmente las necesita, pero podría considerar el uso del método de Newton como una posible alternativa. Sin embargo, requeriría varios ciclos de multiplicación / división.

Recursos AVR

• Nota de la aplicación Atmel AVR200: Multiplicar y dividir rutinas (pdf)
• Esta función sqrt en el foro AVR Freaks
• Another publicación de AVR Freaks

Answer 2

Considere el uso de métodos CORDIC. El Dr. Dobb''s tiene un artículo y una fuente de biblioteca asociada here . La raíz cuadrada, la multiplicación y la división se tratan al final del artículo.

Answer 3

Para el registro, aquí hay algunas aproximaciones más, que se enumeran en orden de complejidad y precisión. Todos estos suponen 0 ≤ a ≤ b.

• h = b + 0.337 * a // max error ≈ 5.5 %
• h = max(b, 0.918 * (b + (a>>1))) // max error ≈ 2.6 %
• h = b + 0.428 * a * a / b // max error ≈ 1.04 %

Edición : para responder a la pregunta de Ecir Hana, aquí es cómo derivé estas aproximaciones.

Primer paso Aproximar una función de dos variables puede ser un problema complejo. Así, primero transformé esto en el problema de aproximar una función de una variable. Esto se puede hacer eligiendo el lado más largo como un factor de "escala", de la siguiente manera:

h = √ (b ² + a ² )
= b √ (1 + (a / b) ² )
= bf (a / b) donde f (x) = √ (1 + x ² )

Agregar la restricción 0 ≤ a ≤ b significa que solo nos interesa la aproximación de f (x) en el intervalo [0, 1].

A continuación se muestra la gráfica de f (x) en el intervalo relevante, junto con la aproximación dada por Matthew Slattery (a saber (√2−1) x + 1).

Segundo paso El siguiente paso es mirar esta trama, mientras se hace la pregunta "¿Cómo puedo aproximarme a esta función a bajo costo?". Como la curva parece aproximadamente parabólica, mi primera idea fue usar una función cuadrática (tercera aproximación). Pero como esto todavía es relativamente caro, también observé aproximaciones lineales y por partes lineales. Aquí están mis tres soluciones:

Las constantes numéricas (0.337, 0.918 y 0.428) fueron inicialmente parámetros libres. Los valores particulares se eligieron para minimizar el error absoluto máximo de las aproximaciones. La minimización sin duda podría hacerse con algún algoritmo, pero simplemente lo hice "a mano", trazando el error absoluto y ajustando la constante hasta que se minimice. En la práctica esto funciona bastante rápido. Escribir el código para automatizar esto hubiera tomado más tiempo.

El tercer paso es volver al problema inicial de aproximar una función de dos variables:

• h ≈ b (1 + 0.337 (a / b)) = b + 0.337 a
• h ≈ b max (1, 0.918 (1 + (a / b) / 2)) = max (b, 0.918 (b + a / 2))
• h ≈ b (1 + 0.428 (a / b) ² ) = b + 0.428 a ² / b

Answer 4

Puede comenzar por reevaluar si necesita el sqrt . Muchas veces está calculando la hipotenusa para compararla con otro valor; si ajusta el valor que está comparando, puede eliminar la raíz cuadrada por completo.

Answer 5

Si el resultado no tiene que ser particularmente preciso, puede obtener una aproximación simple de manera simple:

Tome los valores absolutos de a y b , e intercambie si es necesario para que tenga a <= b . Entonces:

h = ((sqrt(2) - 1) * a) + b

Para ver de forma intuitiva cómo funciona esto, considere la forma en que se traza una línea en ángulo poco profunda en una pantalla de píxeles (por ejemplo, utilizando el algoritmo de Bresenham). Se ve algo como esto:

+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
| | | | | | | | | | | | | | | | |*|*|*|    ^
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+    |
| | | | | | | | | | | | |*|*|*|*| | | |    |
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+    |
| | | | | | | | |*|*|*|*| | | | | | | | a pixels
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+    |
| | | | |*|*|*|*| | | | | | | | | | | |    |
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+    |
|*|*|*|*| | | | | | | | | | | | | | | |    v
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
 <-------------- b pixels ----------->

Para cada paso en la dirección b , el siguiente píxel a trazar es inmediatamente a la derecha, o un píxel arriba y a la derecha.

La línea ideal de un extremo al otro se puede aproximar por la ruta que une el centro de cada píxel al centro del adyacente. Esta es una serie de segmentos de longitud sqrt(2) y ba segmentos de longitud 1 (tomando un píxel como la unidad de medida). De ahí la fórmula anterior.

Esto da claramente una respuesta precisa para a == 0 y a == b ; pero da una sobreestimación de los valores intermedios.

El error depende de la relación b/a ; el error máximo se produce cuando b = (1 + sqrt(2)) * a resulta ser 2/sqrt(2+sqrt(2)) , o aproximadamente 8.24% sobre el valor verdadero. Eso no es genial, pero si es lo suficientemente bueno para su aplicación, este método tiene la ventaja de ser simple y rápido. (La multiplicación por una constante se puede escribir como una secuencia de turnos y adiciones).

Answer 6

Tal vez podría usar algunas de las bibliotecas de ensambladores Elm Chans y adaptar la función ihypot a su ATtiny. Necesitará reemplazar el MUL y tal vez (no he revisado) algunas otras instrucciones.

Answer 7

Una posibilidad se ve así:

#include <math.h>

/* Iterations   Accuracy
 *  2          6.5 digits
 *  3           20 digits
 *  4           62 digits
 * assuming a numeric type able to maintain that degree of accuracy in
 * the individual operations.
 */
#define ITER 3

double dist(double P, double Q) {
/* A reasonably robust method of calculating `sqrt(P*P + Q*Q)''
 *
 * Transliterated from _More Programming Pearls, Confessions of a Coder_
 * by Jon Bentley, pg. 156.
 */

    double R;
    int i;

    P = fabs(P);
    Q = fabs(Q);

    if (P<Q) {
        R = P;
        P = Q;
        Q = R;
    }

/* The book has this as:
 *  if P = 0.0 return Q; # in AWK
 * However, this makes no sense to me - we''ve just insured that P>=Q, so
 * P==0 only if Q==0;  OTOH, if Q==0, then distance == P...
 */
    if ( Q == 0.0 )
        return P;

    for (i=0;i<ITER;i++) {
        R = Q / P;
        R = R * R;
        R = R / (4.0 + R);
        P = P + 2.0 * R * P;
        Q = Q * R;
    }
    return P;
}

Esto todavía hace un par de divisiones y cuatro múltiplos por iteración, pero rara vez se necesitan más de tres iteraciones (y dos a menudo son adecuadas) por entrada. Al menos con la mayoría de los procesadores que he visto, eso generalmente será más rápido que el sqrt por sí solo.

Por el momento está escrito para el double s, pero suponiendo que haya implementado las operaciones básicas, convertirlo para que funcione con un punto fijo no debería ser muy difícil.