2014 dxdy logo

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

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




Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 10:59 


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

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

 Профиль  
                  
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 12:33 
Заслуженный участник


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

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


17/09/09
5
Хм...

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

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

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


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

 Профиль  
                  
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 15:51 
Аватара пользователя


31/10/08
1244
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 


17/09/09
5
Цитата:
Но оптимизировать можно
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 
Заслуженный участник


27/04/09
28128
Вот как для двумерного вектора оптимизировали в 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 
Аватара пользователя


31/10/08
1244
arseniiv
Это трюк для повышения точности, а не оптимизация по скорости.
Цитата:
- уменьшается кол-во возведений в квадрат на 1.
нет как было 2 умножения так и осталось, а еще и деление добавилось.

 Профиль  
                  
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение17.09.2009, 19:12 
Заслуженный участник


27/04/09
28128
Ааа, а я думал-думал, в чём же оптимизация. А это и не оптимизация... Э :oops:

 Профиль  
                  
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение18.09.2009, 10:28 
Заслуженный участник


11/05/08
32166
Pavia в сообщении #244191 писал(а):
нет как было 2 умножения так и осталось, а еще и деление добавилось.

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

 Профиль  
                  
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение18.09.2009, 20:50 
Заслуженный участник


26/07/09
1559
Алматы
2flerant
Цитата:
Есть ли какие-то алгоритмы, которые могли бы сделать это быстрее, чем банальное вычисление квадратного корня из суммы квадратов?

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

 Профиль  
                  
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение18.09.2009, 20:55 
Заслуженный участник


04/05/09
4587
Circiter в сообщении #244522 писал(а):
2flerant
Цитата:
Есть ли какие-то алгоритмы, которые могли бы сделать это быстрее, чем банальное вычисление квадратного корня из суммы квадратов?

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

 Профиль  
                  
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение18.09.2009, 22:15 


17/09/09
5
Circiter в сообщении #244522 писал(а):
Иногда можно корень просто игнорировать. Например, если вам нужно сравнить длины двух векторов, то необязательно сравнивать корни из сумм квадратов координат, достаточно сравнить сами суммы квадратов... Можно и другие действия с длинами оптимизировать повыбрасывав квадратные корни (корни нужно будет извлекать только в самом конце, при выводе данных)...


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

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


Аналогично.

 Профиль  
                  
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение19.09.2009, 01:57 
Заслуженный участник


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

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

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

 Профиль  
                  
 
 Re: Алгоритм быстрого вычисления длины вектора
Сообщение19.09.2009, 05:26 
Заслуженный участник


26/07/09
1559
Алматы
Цитата:
Иногда возможно эти другие выражения преобразовать, чтобы уменьшить количество медленных операций.

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

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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 16 ]  На страницу 1, 2  След.

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group