2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




На страницу 1, 2  След.
 
 Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 10:59 
Есть такая задача. Имеются координаты вектора в 3х или 4х мерном пространстве. Необходимо как можно эффективнее (с точки зрения скорости) вычислить его длину. Есть ли какие-то алгоритмы, которые могли бы сделать это быстрее, чем банальное вычисление квадратного корня из суммы квадратов?

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

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

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 12:51 
Хм...

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

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

код: [ скачать ] [ спрятать ]
Используется синтаксис 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

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 14:25 
Аватара пользователя
Маленький тестик показал, что на моём P4 fsqrt (pi) выполняется примерно за 40 тактов.

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 15:51 
Аватара пользователя
ewert в сообщении #244042 писал(а):
Вряд ли. Даже операция извлечения корня по нынешним временам занимает не более чем несколько тактов (а может, и всего один, не знаю). На чём тут можно сэкономить?...

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

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

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

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 17:14 
Цитата:
Но оптимизировать можно
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 функция проигрывает вдвое.

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 18:16 
Вот как для двумерного вектора оптимизировали в 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)} $, но это наверно хуже будет, наоборот, чем для двумерного случая - теперь слишком много делений. Но двумерная формула точно лучше работает

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 19:07 
Аватара пользователя
arseniiv
Это трюк для повышения точности, а не оптимизация по скорости.
Цитата:
- уменьшается кол-во возведений в квадрат на 1.
нет как было 2 умножения так и осталось, а еще и деление добавилось.

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 19:12 
Ааа, а я думал-думал, в чём же оптимизация. А это и не оптимизация... Э :oops:

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение18.09.2009, 10:28 
Pavia в сообщении #244191 писал(а):
нет как было 2 умножения так и осталось, а еще и деление добавилось.

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

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение18.09.2009, 20:50 
2flerant
Цитата:
Есть ли какие-то алгоритмы, которые могли бы сделать это быстрее, чем банальное вычисление квадратного корня из суммы квадратов?

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

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение18.09.2009, 20:55 
Circiter в сообщении #244522 писал(а):
2flerant
Цитата:
Есть ли какие-то алгоритмы, которые могли бы сделать это быстрее, чем банальное вычисление квадратного корня из суммы квадратов?

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

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение18.09.2009, 22:15 
Circiter в сообщении #244522 писал(а):
Иногда можно корень просто игнорировать. Например, если вам нужно сравнить длины двух векторов, то необязательно сравнивать корни из сумм квадратов координат, достаточно сравнить сами суммы квадратов... Можно и другие действия с длинами оптимизировать повыбрасывав квадратные корни (корни нужно будет извлекать только в самом конце, при выводе данных)...


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

venco в сообщении #244528 писал(а):
Я ещё и деления стараюсь умножениями заменить.


Аналогично.

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

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

flerant в сообщении #244569 писал(а):
Корни нужны прямо здесь и сейчас и именно их значения - они далее входят в другие выражения.
Иногда возможно эти другие выражения преобразовать, чтобы уменьшить количество медленных операций.

 
 
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение19.09.2009, 05:26 
Цитата:
Иногда возможно эти другие выражения преобразовать, чтобы уменьшить количество медленных операций.

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

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

 
 
 [ Сообщений: 16 ]  На страницу 1, 2  След.


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group