Rating@Mail.ru

Форум по операционной системе GNU/Linux и свободному программному обеспечению


Текущее время: 24 ноя 2017, 08:16

Часовой пояс: UTC + 3 часа




Начать новую тему Ответить на тему  [ Сообщений: 7 ] 
Автор Сообщение
 Заголовок сообщения: алгоритмические трюки
Непрочитанное сообщениеДобавлено: 29 авг 2017, 21:07 
Не в сети
Писатель
Аватара пользователя

Зарегистрирован: 24 сен 2011, 14:22
Сообщения: 10228
Откуда: Харьков
Есть такие не очевидные трюки, которые позволяют иногда на порядок ускорить вычисление результата.
Обычно, они (такие трюки) относятся к очень часто используемых вычислений.
Профессиональному программисту нужно бы такие вычисления-трюки знать наперёд.
Поэтому для их фиксации - вот такая отдельная тема.
Отдельная тема ещё потому, что такие трюки не зависят от языка программирования, могут использоваться везде ... и було бы неразумно включать их в существующие темы форума, посвящённые конкретным языкам программирования.


Вернуться к началу
 Профиль Отправить личное сообщение  
 
 Заголовок сообщения: Re: алгоритмические трюки
Непрочитанное сообщениеДобавлено: 29 авг 2017, 21:12 
Не в сети
Писатель
Аватара пользователя

Зарегистрирован: 24 сен 2011, 14:22
Сообщения: 10228
Откуда: Харьков
Olej писал(а):
Есть такие не очевидные трюки, которые позволяют иногда на порядок ускорить вычисление результата.
Обычно, они (такие трюки) относятся к очень часто используемых вычислений.

«Магическая константа» 0x5f3759df
Цитата:
В этой статье мы поговорим о «магической» константе 0x5f3759df, лежащей в основе элегантного алгоритмического трюка для быстрого вычисления обратного квадратного корня.

Вот полная реализация этого алгоритма:
Код:
float FastInvSqrt(float x) {
  float xhalf = 0.5f * x;
  int i = *(int*)&x;  // представим биты float в виде целого числа
  i = 0x5f3759df - (i >> 1);  // какого черта здесь происходит ?
  x = *(float*)&i;
  x = x*(1.5f-(xhalf*x*x));
  return x;
}


Зачем?:
Цитата:
Зачем вообще нам может понадобиться считать обратный квадратный корень, да ещё и пытаться делать это настолько быстро, что нужно реализовывать для этого специальные хаки? Потому, что это одна из основных операций в 3D программировании. При работе с 3D графикой используются нормали поверхностей. Это вектор (с тремя координатами) длиной в единицу, который нужен для описания освещения, отражения и т.д. Таких векторов нужно много. Нет, даже не просто много, а МНОГО. Как мы нормализуем (приводим длину к единице) вектор? Мы делим каждую координату на текущую длину вектора. Ну или, перефразируя, умножаем каждую координату вектора на величину:
Изображение
Расчёт image относительно прост и работает быстро на всех типах процессоров. А вот расчёт квадратного корня и деление – дорогие операции. И вот поэтому – встречайте алгоритм быстрого извлечения обратного квадратного корня — FastInvSqrt.


А то, почему такой "дикий" код, смешивающий int и float представления, работает - подробно анализируется в статье.


Вернуться к началу
 Профиль Отправить личное сообщение  
 
 Заголовок сообщения: Re: алгоритмические трюки
Непрочитанное сообщениеДобавлено: 30 авг 2017, 19:14 
Не в сети
Писатель
Аватара пользователя

Зарегистрирован: 24 сен 2011, 14:22
Сообщения: 10228
Откуда: Харьков
Ещё один очень частый случай:
- В обработке временных рядов, цифровой обработке сигналов и для других целей часто вычисляется среднее и дисперсия (и, далее, среднеквадратичное отклонение) числовой последовательности X[ i ]: mean = 1 / N * ∑( X[ i ] ), disp = 1 / N * ∑( ( X[ i ] - mean )^2 ). Но прямое вычисление характеристик по математическим формулам требует 2-х проходов: сначала вычисление среднего, а затем уже — дисперсии.
- В этом есть много минусов: в 2 раза больше перебор, промахи попаданий в кэш при повторном проходе ... но самое худшее при этом то, что нужно хранить в памяти задачи всю последовательность чисел, что при больших N, например, 10M double ... или просто непрерывный поток входных отсчётов, неприемлемо.
Т.е. такое вычисление "в лоб", из учебника, выглядит примерно как-то так:
Код:
#define n 1000000
double x[ n ], s1 = 0.0, s2 = 0.0;
for( int i = 0; i < n; i++ )
   s1 += x[ i ];
s1 /= n;
for( int i = 0; i < n; i++ )
   s2 += ( x[ i ] - s1 ) * ( x[ i ] - s1 );
s2 = s2 / n

Здесь: s1 - среднее, постоянная составляющая сигнала, s2 - дисперсия, мера разброса, отклонения от среднего, мощность сигнала.

Если раскрыть выражение (квадрат разности) для дисперсии и произвести дальнейшие упрощения, то можно прийти к выражению: disp = 1 / N * ∑( X[ i ] 2 ) - mean2. Теперь мы можем накапливать сумму квадратов последовательности чисел в потоке, а вычисление характеристик отложить на завершение:
Код:
#define n 1000000
double x[ n ], s1 = 0.0, s2 = 0.0;
for( int i = 0; i < n; i++ ) {
   s1 += d;
   s2 += d * d;
}
s1 /= n;
s2 = s2 / n - s1 * s1;


Если же вам нужно СКО, линейное среднеквадратичное отклонение ряда, то оно:
Код:
s2 = sqr( s2 );


Вернуться к началу
 Профиль Отправить личное сообщение  
 
 Заголовок сообщения: Re: алгоритмические трюки
Непрочитанное сообщениеДобавлено: 01 сен 2017, 19:33 
Не в сети
Писатель
Аватара пользователя

Зарегистрирован: 24 сен 2011, 14:22
Сообщения: 10228
Откуда: Харьков
Olej писал(а):
А то, почему такой "дикий" код, смешивающий int и float представления, работает - подробно анализируется в статье.

Код:
#include <stdio.h>
#include <math.h>

float FastInvSqrt( float x ) {
  float xhalf = 0.5f * x;
  int i = *(int*)&x;                    // представим биты float в виде целого числа
  i = 0x5f3759df - ( i >> 1 );          // какого черта здесь происходит ?
  x = *(float*)&i;
  x = x * ( 1.5f - ( xhalf * x * x ) ); // одна итерация уточнения (метод Ньютона)
  return x;
}

int main( void ) {
   float x[] = { .1, 1., 20, 50, 100, 200 };
   for( int i = 0; i < sizeof( x ) / sizeof( *x ); i++ ) {
      float s1 = 1. / sqrt( x[ i ] ), s2 = FastInvSqrt( x[ i ] );
      printf( "%8.2f:\t%f\t%f\t[%f%]\n", x[ i ], s1, s2, ( s1 / s2 - 1. ) * 100 );
   }
   return 0;
}

Код:
[olej@dell sqrt]$ gcc -lm fsqrti.c -o fsqrti

[olej@dell sqrt]$ ./fsqrti
    0.10:   3.162278   3.157232   [0.159812%]
    1.00:   1.000000   0.998307   [0.169575%]
   20.00:   0.223607   0.223571   [0.016224%]
   50.00:   0.141421   0.141356   [0.046277%]
  100.00:   0.100000   0.099845   [0.155365%]
  200.00:   0.070711   0.070678   [0.046277%]

Как легко видеть - точность быстрого вычисление не хуже 0.15% относительно точного значения.


Вернуться к началу
 Профиль Отправить личное сообщение  
 
 Заголовок сообщения: Re: алгоритмические трюки
Непрочитанное сообщениеДобавлено: 01 сен 2017, 19:38 
Не в сети
Писатель
Аватара пользователя

Зарегистрирован: 24 сен 2011, 14:22
Сообщения: 10228
Откуда: Харьков
Olej писал(а):
А то, почему такой "дикий" код, смешивающий int и float представления, работает - подробно анализируется в статье.

По мотивам этой статьи можно построить и вычисление корня квадратного (без инверсии, 1/sqrt()):
Код:
#include <stdio.h>
#include <math.h>

float FastInvSqrt( float x ) {
  int i = *(int*)&x;                 // представим биты float в виде целого числа
  i = 0x1fbd1df5 + ( i >> 1 );
  float y = *(float*)&i;
  y = y - ( y * y - x ) / 2. / y;    // одна итерация уточнения (метод Ньютона)
  return y;
}

int main( void ) {
   float x[] = { .1, 1., 20, 50, 100, 200 };
   for( int i = 0; i < sizeof( x ) / sizeof( *x ); i++ ) {
      float s1 = sqrt( x[ i ] ), s2 = FastInvSqrt( x[ i ] );
      printf( "%8.2f:\t%f\t%f\t[%f%]\n", x[ i ], s1, s2, ( s1 / s2 - 1. ) * 100 );
   }
   return 0;
}

Код:
[olej@dell sqrt]$ gcc -lm fsqrt.c -o fsqrt

[olej@dell sqrt]$ ./fsqrt
    0.10:       0.316228        0.316243        [-0.004882%]
    1.00:       1.000000        1.000064        [-0.006413%]
   20.00:       4.472136        4.472575        [-0.009817%]
   50.00:       7.071068        7.071161        [-0.001317%]
  100.00:       10.000000       10.000242       [-0.002420%]
  200.00:       14.142136       14.142322       [-0.001317%]

Если даже исключить (закомментировать) 1 шаг уточнения итерационным алгоритмом Ньютона:
Код:
  y = y - ( y * y - x ) / 2. / y;    // одна итерация уточнения (метод Ньютона)

... результат остаётся достаточно удовлетворительным для многих целей:
Код:
[olej@dell sqrt]$ ./fsqrt
    0.10:       0.316228        0.319369        [-0.983626%]
    1.00:       1.000000        0.988738        [1.138985%]
   20.00:       4.472136        4.409907        [1.411116%]
   50.00:       7.071068        7.034907        [0.514019%]
  100.00:       10.000000       10.069814       [-0.693297%]
  200.00:       14.142136       14.069814       [0.514019%]


Вернуться к началу
 Профиль Отправить личное сообщение  
 
 Заголовок сообщения: Re: алгоритмические трюки
Непрочитанное сообщениеДобавлено: 01 сен 2017, 19:42 
Не в сети
Писатель
Аватара пользователя

Зарегистрирован: 24 сен 2011, 14:22
Сообщения: 10228
Откуда: Харьков
Olej писал(а):
По мотивам этой статьи можно построить и вычисление корня квадратного (без инверсии, 1/sqrt()):

Аналогично можно извлекать и корень кубический (степень 0.333):
Код:
#include <stdio.h>
#include <math.h>

float FastInvSqrt( float x ) {
  int i = *(int*)&x;                 // представим биты float в виде целого числа
  i = (int)( 0x2a517d3c + ( 0.333f * i ) );
  x = *(float*)&i;
  return x;
}

int main( void ) {
   float x[] = { .1, 1., 20, 50, 100, 200 };
   for( int i = 0; i < sizeof( x ) / sizeof( *x ); i++ ) {
      float s1 = powf( x[ i ], 0.333f ),
            s2 = FastInvSqrt( x[ i ] );
      printf( "%8.2f:\t%f\t%f\t[%f%]\n", x[ i ], s1, s2, ( s1 / s2 - 1. ) * 100 );
   }
   return 0;
}

Код:
[olej@dell sqrt]$ gcc -lm fsqrt3.c -o fsqrt3

[olej@dell sqrt]$ ./fsqrt3
    0.10:       0.464515        0.448858        [3.488183%]
    1.00:       1.000000        0.963818        [3.754067%]
   20.00:       2.711709        2.685760        [0.966132%]
   50.00:       3.679231        3.559906        [3.351903%]
  100.00:       4.634469        4.451782        [4.103673%]
  200.00:       5.837716        5.783813        [0.931954%]

Это показано без последнего шага итерационного уточнения.
Можете это доделать сами ... метод Ньютона описан здесь.


Вернуться к началу
 Профиль Отправить личное сообщение  
 
 Заголовок сообщения: Re: алгоритмические трюки
Непрочитанное сообщениеДобавлено: 01 сен 2017, 19:52 
Не в сети
Писатель
Аватара пользователя

Зарегистрирован: 24 сен 2011, 14:22
Сообщения: 10228
Откуда: Харьков
Olej писал(а):
А то, почему такой "дикий" код, смешивающий int и float представления, работает - подробно анализируется в статье.

В комментариях публикации обстоятельно обсуждается, что это всё работает, если компилятору С/C++ установлена опция -no-strict-aliasing (или по умолчанию установлено). Или на другой аппаратной платформе, отличной от x86.
Это достаточно интересно...
Но у GCC нет этой опции:
Код:
[olej@dell sqrt]$ gcc -no-strict-aliasing -lm fsqrt.c -o fsqrt
gcc: ошибка: unrecognized command line option «-no-strict-aliasing»

Скорее всего, обсуждение касается LLVM.
Подробно разъяснения по опции см. По следам C++ Siberia: дракон в мешке.
В конце концов ... если есть опасения, проблема решается (один из способов) так:
Код:
#include <stdio.h>
#include <math.h>
#include <string.h>

float FastInvSqrt( float x ) {
  int i;
  memcpy( &i, &x, sizeof( float ) ); // представим биты float в виде целого числа
  i = 0x1fbd1df5 + ( i >> 1 );
  float y;
  memcpy( &y, &i, sizeof( float ) );
  y = y - ( y * y - x ) / 2. / y;    // одна итерация уточнения (метод Ньютона)
  return y;
}

int main( void ) {
   float x[] = { .1, 1., 20, 50, 100, 200 };
   for( int i = 0; i < sizeof( x ) / sizeof( *x ); i++ ) {
      float s1 = sqrt( x[ i ] ), s2 = FastInvSqrt( x[ i ] );
      printf( "%8.2f:\t%f\t%f\t[%f%]\n", x[ i ], s1, s2, ( s1 / s2 - 1. ) * 100 );
   }
   return 0;
}

Код:
[olej@dell sqrt]$ gcc -lm fsqrtc.c -o fsqrtc

[olej@dell sqrt]$ ./fsqrtc
    0.10:       0.316228        0.316243        [-0.004882%]
    1.00:       1.000000        1.000064        [-0.006413%]
   20.00:       4.472136        4.472575        [-0.009817%]
   50.00:       7.071068        7.071161        [-0.001317%]
  100.00:       10.000000       10.000242       [-0.002420%]
  200.00:       14.142136       14.142322       [-0.001317%]


Вернуться к началу
 Профиль Отправить личное сообщение  
 
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 7 ] 

Часовой пояс: UTC + 3 часа


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

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


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

Найти:
Перейти:  
cron
Создано на основе phpBB® Forum Software © phpBB Group
Русская поддержка phpBB
[ Time : 0.306s | 19 Queries | GZIP : On ]