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

Question

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

Answer 1

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

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

Answer 2

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

---

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

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

ч = √ (б ² + а ² )
= 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))
• ч ≈ b (1 + 0,428 (a / b) ² ) = b + 0,428 a ² / b

Answer 3

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

Answer 4

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

#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 с, но при условии, что вы реализовали основные операции, преобразование его для работы с фиксированной точкой не должно быть ужасно сложным.

Answer 5

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

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

Ресурсы AVR

• Примечание к приложению Atmel AVR200: процедуры умножения и деления (pdf)
• Эта sqrt функция на форуме AVR Freaks
• Другой AVR Freaks post

Answer 6

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

Answer 7

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