Алгоритм быстрого алгоритма Hypotenuse для встроенного процессора?

Есть ли умный/эффективный алгоритм определения гипотенузы угла (т.е. sqrt(a² + b²)), используя математику с фиксированной точкой во встроенном процессоре без аппаратного умножения?

Ответ 1

Если вы делаете это при > 1 кГц, умножить даже на MCU без аппаратного обеспечения MUL не страшно. Что еще хуже, это sqrt. Я попытался бы изменить мое приложение, поэтому его вообще не нужно вычислять.

Стандартные библиотеки, вероятно, были бы лучше всего, если бы вы на самом деле нуждались в этом, но вы могли бы использовать метод Newton как возможную альтернативу. Однако для выполнения потребуется несколько циклов умножения/деления.

Ресурсы AVR

Заметка приложения Atmel AVR200: операции умножения и разделения (pdf)
Эта функция sqrt на форуме AVR Freaks
Другой пост AVR Freaks

Ответ 2

Если результат не должен быть особенно точным, вы можете получить сырую нефть приближение довольно просто:

Возьмите абсолютные значения a и b и поменяйте, если необходимо, чтобы у вас a <= b. Тогда:

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

Чтобы интуитивно понять, как это работает, рассмотрите способ построения мелкоугольной линии на пиксельном дисплее (например, с использованием алгоритма Брешенема). Это выглядит примерно так:

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

Для каждого шага в направлении b следующий пиксел, который должен быть нанесен, либо сразу справа, либо один пиксель вверх и вправо.

Идеальная линия от одного конца до другого может быть аппроксимирована путем, который соединяет центр каждого пикселя с центром соседнего. Это серия a сегментов длины sqrt(2) и b-a сегментов длины 1 (взятие пикселя в качестве единицы измерения). Следовательно, приведенная выше формула.

Это явно дает точный ответ для a == 0 и a == b; но дает переоценку значений между ними.

Ошибка зависит от отношения b/a; максимальная ошибка возникает, когда b = (1 + sqrt(2)) * a и оказывается 2/sqrt(2+sqrt(2)), или около 8,24% над истинным значением. Это не очень хорошо, но если это достаточно хорошо для вашего приложения, этот метод имеет то преимущество, что он прост и быстр. (Умножение на константу можно записать в виде последовательности сдвигов и добавляет.)

Ответ 3

Для записи здесь приведены еще несколько приближений, возрастающий порядок сложности и точности. Все это предполагает, что 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 %

Изменить: ответить на вопрос Эрира Ханы, вот как я получил эти приближения.

Первый шаг. Приближение функции двух переменных может быть сложная проблема. Таким образом, я впервые превратил это в проблему аппроксимируя функцию одной переменной. Это можно сделать, выбирая самая длинная сторона как "масштабный" коэффициент, следующим образом:

h = √ (b ² + a ²)
= b √ (1 + (a/b) ²)
= b f (a/b) , где f (x) = √ (1 + x ²)

Добавление ограничения 0 ≤ a ≤ b означает, что мы имеем дело только с аппроксимируя f (x) в интервале [0, 1].

Ниже приведен график f (x) в соответствующем интервале, вместе с аппроксимация, приведенная Мэтью Слэттери (а именно (√2-1) x + 1).

Второй шаг. Следующий шаг - посмотреть на этот сюжет, задавая сам вопрос: "Как я могу приблизить эту функцию дешево?". Поскольку кривая выглядит грубо параболической, моя первая идея заключалась в том, чтобы использовать квадратичная функция (третье приближение). Но так как это все еще относительно дорого, я также посмотрел на линейные и кусочно-линейные приближения. Вот мои три решения:

Численные константы (0.337, 0.918 и 0.428) первоначально были свободны параметры. Конкретные значения были выбраны для минимизации максимальная абсолютная ошибка приближений. Минимизация могла конечно, нужно сделать с помощью некоторого алгоритма, но я просто сделал это "вручную", построение абсолютной ошибки и настройку константы до тех пор, пока она не будет сведено к минимуму. На практике это работает довольно быстро. Написание кода для автоматизировать это заняло бы больше времени.

Третий шаг - вернуться к исходной проблеме аппроксимации функция двух переменных:

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

Ответ 4

Рассмотрим методы CORDIC. У доктора Добба есть статья и связанный с ней источник библиотеки здесь. Квадратный корень, умножение и деление рассматриваются в конце статьи.

Ответ 5

Одна возможность выглядит так:

#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;
}

Это все еще делает пару делений и четыре мультипликатора на итерацию, но вам редко требуется больше трех итераций (а два - достаточно адекватных) для каждого ввода. По крайней мере, с большинством процессоров, которые я видел, это будет, как правило, быстрее, чем sqrt будет само по себе.

На данный момент он написан для double s, но если вы выполнили основные операции, преобразование его в работу с фиксированной точкой не должно быть ужасно трудным.

Ответ 6

Вы можете начать с переоценки, если вам нужен sqrt вообще. Много раз вы вычисляете гипотенузу, чтобы сравнить ее с другим значением - если вы сравните значение, которое вы сравниваете, вы можете полностью исключить квадратный корень.

Ответ 7

Возможно, вы могли бы использовать некоторые из Elm Chans Ассемблерные библиотеки и адаптировать ihypot-функцию к вашему ATtiny. Вам нужно будет заменить MUL и, возможно, (я еще не проверил) некоторые другие инструкции.