Алгоритм быстрого вычисления длины вектора

[flerant]

Есть такая задача. Имеются координаты вектора в 3х или 4х мерном пространстве. Необходимо как можно эффективнее (с точки зрения скорости) вычислить его длину. Есть ли какие-то алгоритмы, которые могли бы сделать это быстрее, чем банальное вычисление квадратного корня из суммы квадратов?

Особый интерес для меня представляет вычисление не самой длины, а величины, обратной ей.

[ewert]

Вряд ли. Даже операция извлечения корня по нынешним временам занимает не более чем несколько тактов (а может, и всего один, не знаю). На чём тут можно сэкономить?...

[flerant]

Хм...

Мне удавалось экономить даже на замене деления на умножение. А взятие квадратного корня, насколько мне известно, вообще весьма затратная вещь.

Вот простой пример:

Используется синтаксис C++

#include <cmath>
#include <iostream>
#include <iomanip>
#include <ctime>

using namespace std;

int main()
{
        double x[500000];
        double y[500000];
        clock_t t0 = clock();
        for (int j = 1; j <= 1000; j++)
        for (int i = 1; i <= 500000; i++)
                x[i-1] = sqrt(1. + (double)i*1.e-10);
        clock_t t1 = clock();
        for (int j = 1; j <= 1000; j++)
        for (int i = 1; i <= 500000; i++)
                y[i-1] = 1. + (double)i*.5e-10;
        clock_t t2 = clock();
        
        double a = 0.;
        for (int i = 0; i < 500000; i++)
                a += y[i] - x[i];
        
        cout << setw(17) << a << setw(8)<< (t1-t0)*invcps << setw(8) << (t2-t1)*invcps << endl;
        
        return 0;
}
 

Вывод:

      5.20825e-05   14.62    0.96

[worm2]

Маленький тестик показал, что на моём P4 fsqrt (pi) выполняется примерно за 40 тактов.

[Pavia]

ewert в сообщении #244042 писал(а):

Вряд ли. Даже операция извлечения корня по нынешним временам занимает не более чем несколько тактов (а может, и всего один, не знаю). На чём тут можно сэкономить?...

Везет, на 45нм core процессоре FSQRT занимает 6 тактов. На всех остальных 40-60 тактов SQRT вычисляется по битно, Это значит чтобы вычислить 1 бит нужно 1 умножение или 1 такт.

Но оптимизировать можно
http://en.wikipedia.org/wiki/Fast_inverse_square_root

Можно взять MKL там уже оптимизировано.

[flerant]

Цитата:

Но оптимизировать можно
http://en.wikipedia.org/wiki/Fast_inverse_square_root

Спасибо - я посмотрю

Цитата:

Можно взять MKL там уже оптимизировано.

Возможно, я что-то делаю неправильно, но почему-то у меня эта библиотека считает медленнее, чем обычный корень из суммы квадратов Вот мой код

Используется синтаксис C++

#include <cmath>
#include <iostream>
#include <iomanip>
#include <ctime>
#include <mkl_blas.h>

using namespace std;

float invcps = 1./(float)CLOCKS_PER_SEC;

int main()
{
        double x[50000];
        double y[50000];
        double a[4];
        a[0] = 1.;
        MKL_INT n = 4;
        MKL_INT d = 1;

        clock_t t0 = clock();
        for (int j = 1; j <= 1000; j++)
        {
                for (int i = 1; i <= 50000; i++)
                {
                        a[1] = (double)(i-j);
                        a[2] = (double)(j+i);
                        a[3] = (double)(j/i);
                        x[i-1] = sqrt(1. + a[1]*a[1] + a[2]*a[2] + a[3]*a[3]);
                }
        }
        clock_t t1 = clock();

        for (int j = 1; j <= 1000; j++)
        {
                for (int i = 1; i <= 50000; i++)
                {
                        a[1] = (double)(i-j);
                        a[2] = (double)(j+i);
                        a[3] = (double)(j/i);
                        y[i-1] = DNRM2(&n, a, &d);
                }
        }
        clock_t t2 = clock();
        
        cout << setw(8)<< (t1-t0)*invcps << setw(8) << (t2-t1)*invcps << endl;
        
        return 0;
}
 

Компилирую командой

icpc -O2 -no-multibyte-chars -lmkl -openmp main.cpp -o out

Вывод на экран

    0.24     2.9

Т.е. MKL функция считает чуть ли не на порядок хуже! Конечно, это объясняется тем, что компилятор оптимизирует как-то первый цикл, поскольку он довольно простой. Но даже при его искусственном усложнении MKL функция проигрывает вдвое.

[arseniiv]

Вот как для двумерного вектора оптимизировали в Delphi, причём описан код и ассемблерные команды тоже, для быстроты; может, пригодится:

function Hypot(const X, Y: Extended): Extended;
{ formula: Sqrt(X*X + Y*Y)
  implemented as:  |Y|*Sqrt(1+Sqr(X/Y)), |X| < |Y| for greater precision
var
  Temp: Extended;
begin
  X := Abs(X);
  Y := Abs(Y);
  if X > Y then
  begin
    Temp := X;
    X := Y;
    Y := Temp;
  end;
  if X = 0 then
    Result := Y
  else         // Y > X, X <> 0, so Y > 0
    Result := Y * Sqrt(1 + Sqr(X/Y));
end;
}
asm
        FLD     Y
        FABS
        FLD     X
        FABS
        FCOM
        FNSTSW  AX
        TEST    AH,$45
        JNZ      @@1        // if ST > ST(1) then swap
        FXCH    ST(1)      // put larger number in ST(1)
@@1:    FLDZ
        FCOMP
        FNSTSW  AX
        TEST    AH,$40     // if ST = 0, return ST(1)
        JZ      @@2
        FSTP    ST         // eat ST(0)
        JMP     @@3
@@2:    FDIV    ST,ST(1)   // ST := ST / ST(1)
        FMUL    ST,ST      // ST := ST * ST
        FLD1
        FADD               // ST := ST + 1
        FSQRT              // ST := Sqrt(ST)
        FMUL               // ST(1) := ST * ST(1); Pop ST
@@3:    FWAIT
end;

Формула приведена так: $\sqrt {x^2 + y^2 } = \sqrt {x^2 \left( {1 + \left( {y/x} \right)^2 } \right)} = \left| x \right|\sqrt {1 + \left( {y/x} \right)^2 }$ - уменьшается кол-во возведений в квадрат на 1.
Очевидно, можно что-то и в этом духе сделать: $\sqrt {x_1^2 + x_2^2 + \ldots + x_n^2 } = \left| {x_1 } \right|\sqrt {\left( {1 + \ldots + \left( {x_n /x_1 } \right)^2 } \right)} $, но это наверно хуже будет, наоборот, чем для двумерного случая - теперь слишком много делений. Но двумерная формула точно лучше работает

[Pavia]

arseniiv
Это трюк для повышения точности, а не оптимизация по скорости.

Цитата:

- уменьшается кол-во возведений в квадрат на 1.

нет как было 2 умножения так и осталось, а еще и деление добавилось.

[arseniiv]

Ааа, а я думал-думал, в чём же оптимизация. А это и не оптимизация... Э

[ewert]

Pavia в сообщении #244191 писал(а):

нет как было 2 умножения так и осталось, а еще и деление добавилось.

а ещё и взятие модуля -- тоже не безобидно, между прочим

[Circiter]

2flerant

Цитата:

Есть ли какие-то алгоритмы, которые могли бы сделать это быстрее, чем банальное вычисление квадратного корня из суммы квадратов?

Хоть и глупость, но все же... Иногда можно корень просто игнорировать. Например, если вам нужно сравнить длины двух векторов, то необязательно сравнивать корни из сумм квадратов координат, достаточно сравнить сами суммы квадратов... Можно и другие действия с длинами оптимизировать повыбрасывав квадратные корни (корни нужно будет извлекать только в самом конце, при выводе данных)...

[venco]

Circiter в сообщении #244522 писал(а):

2flerant

Цитата:

Есть ли какие-то алгоритмы, которые могли бы сделать это быстрее, чем банальное вычисление квадратного корня из суммы квадратов?

Хоть и глупость, но все же...

Почему же глупость. Очень полезная оптимизация. Я ещё и деления стараюсь умножениями заменить.

[flerant]

Circiter в сообщении #244522 писал(а):

Иногда можно корень просто игнорировать. Например, если вам нужно сравнить длины двух векторов, то необязательно сравнивать корни из сумм квадратов координат, достаточно сравнить сами суммы квадратов... Можно и другие действия с длинами оптимизировать повыбрасывав квадратные корни (корни нужно будет извлекать только в самом конце, при выводе данных)...

Это я понимаю. Но это не тот случай. Корни нужны прямо здесь и сейчас и именно их значения - они далее входят в другие выражения.

venco в сообщении #244528 писал(а):

Я ещё и деления стараюсь умножениями заменить.

Аналогично.

[venco]

Для сравнения, на современных Intel процессорах сложение/вычитание выполняются за 1 такт, умножение - за 1-2 такта, а деление и корень - за ~40 тактов.
И деление, и корень по прежнему вычисляются итерациями по Ньютону, поэтому и медленно.
Специально для желающих сделать итерации самостоятельно, возможно, не до конца, если не нужна полная точность, в SSE добавили быстрые инструкции для получения начального приближения и для обратного корня, и для обратной величины.

-- Пт сен 18, 2009 19:00:36 --

flerant в сообщении #244569 писал(а):

Корни нужны прямо здесь и сейчас и именно их значения - они далее входят в другие выражения.

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

[Circiter]

Цитата:

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

А ведь действительно, как-правило большинство алгоритмов можно слегка модифицировать; например, для моей идеи с выбрасыванием корней, можно все-таки попробовать заставить алгоритмы производить операции не с $1/|\cdot|$, а с $|\cdot|^2$, это того стоит.

Ещё одна идея: можно хранить длины векторов отдельно от самих векторов и пересчитывать по необходимости.