LINUX.ORG.RU

Отношение расстояний с нормализацией без float и math.h

 


0

3

Как можно отсюда убрать необходимость во float и math.h?

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

uint16_t SmlGradientRadialValue(SmlLine line, SmlPoint target)
{
  /*                          ┌──────┐
                         ┌────┘      └────┐
                       ┌─┘                └─┐
                     ┌─┘                    └─┐
                   ┌─┘                        └─┐
                   │                            │
                   │   __ Gradient line         │
                 ┌─┘  /                         └─┐
                 │   /             A(x1, y1)      │
                 │──────────────╴()               │
                 │                                │
                 └─┐                            ┌─┘
                   │                            │
                   │                T(tx, ty)   │
                   └─┐            ()          ┌─┘
                     └─┐                    ┌─┘
                       └─┐   B(x2, y2)    ┌─┘
                         └─()─┐      ┌────┘
                              └──────┘
   * Returns gradient value [SML_GRAD_MIN..SML_GRAD_MAX] from A to B relative
     to the point T.
   Solution:
   I.  If the target point is equal to A, return SML_GRAD_MIN
   II. Calculate the AB & AT distances.
   III. If AT is bigger than AB, return the SML_GRAD_MAX
   IV. Return the of the AT length to AB length and normalise it. */
 
   if ((target.x == line.p1.x) && (target.y == line.p1.y))
       return SML_GRAD_MIN;
   float AB = sqrt(SML_SQR(line.p1.x - line.p2.x) +
                   SML_SQR(line.p1.y - line.p2.y));
   float AT = sqrt(SML_SQR(line.p1.x - target.x) +
                   SML_SQR(line.p1.y - target.y));
   if (AT > AB)
       return SML_GRAD_MAX;
   return SML_GRAD_MAX * AT / AB;
}
★★
Ответ на: комментарий от sambist

Подозреваю, что внутри sqrt math.h лежит та же аппроксимация методом Ньютона, только оптимизированная под использование конвееров, fpu и прочих плюшек процессора. Поэтому надо смотреть на конкретное железо, что бы добиться хотя бы такой же производительности.

four_str_sam
()
Ответ на: комментарий от false

Вы здесь вообще другие слова кроме этих знаете? Когда я вижу во цифрам, что оно будет лагать при отрисовке, я оптимизирую. Капиш?

sambist ★★
() автор топика
Ответ на: комментарий от hey-propagandist

Это код радиального градиента. Функция возвращает число от 0 до MAX_GRAD, которое означает степень удаления от первой точки (как, например, когда вы в ГИМПе рисуете радиальный градиент). Потом это число отправляется в другую функцию, которая параллельно принимает индекс градиента в хранилище. Результатом является цвет и...

int y, x;
for (y = 0; y < 1000; y++)
for (x = 0; x < 1000; x++)
{
    uint16_t val = SmlGradientRadialValue(line, (SmlPoint){x, y});
    SmlColour col;
    SmlGradientGetColour(grd1, val, &col);
    SmlImagePixelSet(img1, col, (SmlPoint){x, y});
}

Результат.

sambist ★★
() автор топика
Последнее исправление: sambist (всего исправлений: 1)
Ответ на: комментарий от hey-propagandist

Тред почти не читал, за лишнюю подлинкованную библиотеку почему так трясешься? Или у тебя просадка в производительности?

Воображаю себе что когда-нибудь код будет работать на ардуинках.

sambist ★★
() автор топика
Ответ на: комментарий от hey-propagandist

Отказываться от float есть смысл, когда у тебя мало умножений и нет делений, как правило

Если все это сводимо без потери качества к int, который по определению быстрее, то почему нет?

sambist ★★
() автор топика
Ответ на: комментарий от sambist

Тогда и попытки оптимизировать sqrt тоже бессмысленны. Надо заниматься алгоритмами, например тебе уже подсказали, что можно избавиться от одного sqrt в функции.

four_str_sam
()
Ответ на: комментарий от four_str_sam

Где такое определение?

В тестах. Каждый раз когда я успешно меняю float на int скорость вырастает в 3-4 раза. Сейчас застопорилось из-за sqrt, ибо он скорее всего как-то оптимизирован не по героновской формуле как у меня, ибо каогда у меня числа типа 30000, то герону требуется очень много итераций.

sambist ★★
() автор топика
Ответ на: комментарий от four_str_sam

например тебе уже подсказали, что можно избавиться от одного sqrt в функции.

Избавился. Ну стало вместо 42 30 секунд. Все равно до флоатовских 9 как до китая раком.

sambist ★★
() автор топика
Ответ на: комментарий от four_str_sam

http://hijos.ru/2012/04/25/krasivaya-modifikaciya-metoda-izvlecheniya-kvadrat...

Во какую штуку выкопал. Ща вкурю и попробую реализовать. Это не просто «детский» способ, тут вместо всей 1+3+5.. последовательности попарно разбивается...

sambist ★★
() автор топика
Ответ на: комментарий от sambist

У тебя начальное приближение плохое. По ссылке , которую уже приводили, в конце статьи приводится алгоритм с адаптивным начальным приближением. В среднем получается 3 итерации, максимум 5.

four_str_sam
()
Ответ на: комментарий от sambist
:~/Documents/glibc-2.19$ cat sysdeps/i386/fpu/e_sqrtl.c 
/*
 * Written by J.T. Conklin <jtc@netbsd.org>.
 * Public domain.
 *
 * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
 */

#include <math_private.h>

#undef __ieee754_sqrtl
long double
__ieee754_sqrtl (long double x)
{
  long double res;

  asm ("fsqrt" : "=t" (res) : "0" (x));

  return res;
}
strong_alias (__ieee754_sqrtl, __sqrtl_finite)

Вот как бы и ответ на вопрос про скорость.

four_str_sam
()
Ответ на: комментарий от four_str_sam

Замечательно. Почти флоатовский результат. Я очень извиняюсь, что невнимательно прочел статью.

    //float: 100000x1000 cycles took 9.280323 s
    //heron: 100000x1000 cycles took 42.152074 s
    //fast : 100000x1000 cycles took 42.063780 s
    //chld : 100000x1000 cycles took 42.026744 s
    //prer : 100000x1000 cycles took 10.264415 s <<

Еще хочу попробовать все таки реализовать тот арифмометрический алгоритм и на его тайминг посмотреть. Уж очень понравилась фраза из статьи «этот способ легко может быть формализован и записан в виде программного кода для ЭВМ, причём время выполнения этой программы сопоставимо с временем выполнения операции деления».

Большое вам спасибо.

sambist ★★
() автор топика
Ответ на: комментарий от four_str_sam

Да ладно простыня. Там всего 200 строк. Даже с коментариями)

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

nanoolinux ★★★★
()
Ответ на: комментарий от four_str_sam

Подумал немного и модифицировал код из статьи. В общем рекорд в 9 секунд от float побит на целочисленном вычислении:
Код:

uint32_t sqrti(uint32_t value)
{
    uint32_t div;
    if (value & 0xFFFF0000U)
        if (value & 0xFF000000U)
            div = 0x3FFF;
        else
            div = 0x3FF;
    else
        if (value & 0x0000FF00U)
            div = 0x3F;
        else div = (value > 4) ? 0x7 : value;

    uint32_t p;
    uint32_t n = (value / div + div) >> 1;
    do
    {
        p = n;
        n = (p + value / p) >> 1;
    }
    while (p > n);

    return n;
}

Результаты:

    //float : 100000x1000 cycles took 9.280323 s
    //heron : 100000x1000 cycles took 42.152074 s
    //fast  : 100000x1000 cycles took 42.063780 s
    //chld  : 100000x1000 cycles took 42.026744 s
    //prer  : 100000x1000 cycles took 10.264415 s
    //premod: 100000x1000 cycles took 7.893720 s <<

sambist ★★
() автор топика
Последнее исправление: sambist (всего исправлений: 1)
Ответ на: комментарий от four_str_sam

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

sambist ★★
() автор топика
Ответ на: комментарий от four_str_sam

А если вернутся к первоначальной задаче и не отделять вычисление корня из произвольного числа, а использовать как начальное приближение max(abs(x1-x2),abs(y1-y2)), возможно получится быстрее

Elyas ★★★★★
()
Ответ на: комментарий от sambist

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

Elyas ★★★★★
()
Последнее исправление: Elyas (всего исправлений: 1)
Ответ на: комментарий от Elyas
expec: 100000x100 cycles took 11.950104 s

Таки не вышло.

sambist ★★
() автор топика
Ответ на: комментарий от nanoolinux

И такие как я написали PaX, Stack-Smashing Protector и StackGuard. А вы и дальше пытайтесь играть в лагающий тетрис на 16 ГБ памяти и последнем проце.

sambist ★★
() автор топика
Ответ на: комментарий от sambist

Перечёл тред снова и возникло такое соображение: оптимизация вычисления корня для произвольного числа — это хорошо, но нельзя ли при обходе всего массива применять для вычисления расстояний дифференциальный метод? Подобный подход (алгоритм Брезенхама) применялся для слабых видюх при поточечном построении прямых и даже эллипсов.

pS
()
Ответ на: комментарий от pS

Черканул в блокноте, подумаю на досуге над отдельной функцией. Правда там придется float'ами крутить...

sambist ★★
() автор топика
Ответ на: комментарий от pS

Оперируют то целыми, сам писал брезенхема для линии, но учитывая, что разница между соседними точками может быть формально меньше единицы - то тут надо думать. Может быть получится, может нет. По крайней мере на процедуру рисования градиентного прямоугольника пригодится.

sambist ★★
() автор топика
Ответ на: комментарий от anonymous

Большое спасибо, анон. Правда, придется кого-то мне теребить чтобы помог переписать под gcc'шный ассемблер, ибо я в него ни гугу почти.

sambist ★★
() автор топика
Ответ на: комментарий от anonymous

Забавно. На сайте http://hijos.ru/2012/04/25/krasivaya-modifikaciya-metoda-izvlecheniya-kvadrat... который ведет преподавательница из Питера она говорит, что этот алгоритм придумал Савич С.В., который по ее словам, связался с ней чтобы опубликовать этот алгоритм. Цитата:

Вот что пишет сам Сергей Валентинович: “Схема вычисления квадратных корней была придумана мной в конце 70-х г., и мне хотелось даже опубликовать её описание. Но время арифмометров закончилось, и этот алгоритм остался только в моей памяти.

Дата публикации статьи - 2012й год.

А по ссылке, который вы дали тот же самый алгоритм опубликован в Inferno #04, вышедшей в 2003м.

sambist ★★
() автор топика
Ответ на: комментарий от sambist

В общем сделал реализацию кода Савича (хотя я сомневаюсь теперь что это его алгоритм, ибо я нашел его так же в компиляторе 1983 года), 18 секунд. Придется остановится на модифицированном Героне.

Спасибо всем.

P.S. Кроме «ненужнистов»

sambist ★★
() автор топика
Ответ на: комментарий от sambist

А можете выложить код бенчмарка? Хотябы сами функции(насколько я понял premod - это sqrti, которую вы привели выше? Тогда 5).

anonymous
()
Ответ на: комментарий от anonymous

Ну в общем ситуация такая:

В основном приложении:


SmlPoint p1 = (SmlPoint){500, 500};
SmlPoint p2 = (SmlPoint){  0,   0};

void test2()
{
    int32_t i;
    SmlLine line = (SmlLine){p1, p2};
    uint16_t val;
    for (i = 0; i < 1000; i++)
    {
        val += SmlGradientRadialValue(line, (SmlPoint){i, 250});
    }
}

///
SML_PROFILER_INIT;
SML_PROFILER_MEAS(100000, test2(), "expec");

В библиотеке
uint16_t SmlGradientRadialValue(SmlLine line, SmlPoint target)
{
  /*                          ┌──────┐
                         ┌────┘      └────┐
                       ┌─┘                └─┐
                     ┌─┘                    └─┐
                   ┌─┘                        └─┐
                   │                            │
                   │   __ Gradient line         │
                 ┌─┘  /                         └─┐
                 │   /             A(x1, y1)      │
                 │───────────────()               │
                 │                                │
                 └─┐                            ┌─┘
                   │                            │
                   │                T(tx, ty)   │
                   └─┐            ()          ┌─┘
                     └─┐                    ┌─┘
                       └─┐   B(x2, y2)    ┌─┘
                         └─()─┐      ┌────┘
                              └──────┘
   * Returns gradient value [SML_GRAD_MIN..SML_GRAD_MAX] from A to B relative
     to the point T.
   Solution:
   I.  If the target point is equal to A, return SML_GRAD_MIN
   II. Calculate the AB & AT distances.
   III. If AT is bigger than AB, return the SML_GRAD_MAX
   IV. Return the of the AT length to AB length and normalise it. */

   if ((target.x == line.p1.x) && (target.y == line.p1.y))
       return SML_GRAD_MIN;
   uint32_t AB = (SML_SQR(line.p1.x - line.p2.x) +
                  SML_SQR(line.p1.y - line.p2.y));
   uint32_t AT = (SML_SQR(line.p1.x - target.x) +
                  SML_SQR(line.p1.y - target.y));
   if (AT > AB)
       return SML_GRAD_MAX;
   return sqrti(SML_SQR((uint32_t)SML_GRAD_MAX) / AB * AT);
}

профайлер:
/* E+ */
#define SML_PROFILER_INIT            \
    clock_t sml_prof_time_start = 0; \
    clock_t sml_prof_time_stop  = 0; \
    int32_t sml_prof_counter

/// @optimize use system getclock (phone note).
#define SML_PROFILER_MEAS(COUNT, FUNC, MESSAGE)   \
    sml_prof_time_start = clock();                \
    for (sml_prof_counter = 0;                    \
         sml_prof_counter < (COUNT);              \
         sml_prof_counter++)                      \
        (FUNC);                                   \
    sml_prof_time_stop  = clock();                \
    printf(MESSAGE ": %d cycles took %f s\n",     \
        (int32_t)(COUNT),                         \
        (float)(sml_prof_time_stop -              \
                sml_prof_time_start) /            \
               CLOCKS_PER_SEC)
/* E- */
sambist ★★
() автор топика
Последнее исправление: sambist (всего исправлений: 1)
Ответ на: комментарий от sambist

Что-то у меня не получается воспроизвести ваши результаты:

sqrti: 100000 cycles took 1.411691 s
sqrtf(lm): 100000 cycles took 0.536597 s

Т.е. я заменил sqrti(Отношение расстояний с нормализацией без float и math.h (комментарий) отсюда) на корень из lm.

anonymous
()
Ответ на: комментарий от anonymous

Да у меня теперь тоже - 6.4 от fsqrt vs 7.6 sqrti. Видимо в тот раз сказалось какое-то подтормаживание системы. Видимо теперь придется тесты делать дольше.

sambist ★★
() автор топика
Ответ на: комментарий от sambist

В тестах.

Скорее проблема в том, что вы мерили так же, как и в предыдущем измерении.

Каждый раз когда я успешно меняю float на int скорость вырастает в 3-4 раза.

Если я применяю вместо молотка бензопилу - это не значит, что бензопила пилит медленнее ножовки/ручной пилы. Просто я её применяю не так и не там.

Поэтому вся ваша разница - это в 99% случаев следствие ваших реализаций/измерений, а не реального проявления разницы между производительностью флоата и интов.

Все ваши заключение строятся на догадках за которыми никаких объективных матчасть-обоснований попросту нет. Толку от этих рассуждений?

Сейчас застопорилось из-за sqrt, ибо он скорее всего как-то оптимизирован не по героновской формуле как у меня, ибо каогда у меня числа типа 30000, то герону требуется очень много итераций.

А почему вы решили, что цена вешей итерации равна цене итерации в sqrt()?

Так же, у вас изначально тупиковый путь оптимизации.

Допустим тот же градиент, судя по скриншоту:

II. Calculate the AB & AT distances.

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

int y, x;
for (y = 0; y < 1000; y++)
for (x = 0; x < 1000; x++)
{
    uint16_t val = SmlGradientRadialValue(line, (SmlPoint){x, y});
    SmlColour col;
    SmlGradientGetColour(grd1, val, &col);
    SmlImagePixelSet(img1, col, (SmlPoint){x, y});
}

Смысла это оптимизировать вообще нет, ибо в абсолютных значениях это полностью не имеет смысла. Да, это возможно и можно сделать быстрее, но зачем? Оно априори не выдаст даже 10% того, что выдаст изначально вменяемый подход.

Две основные ошибки - попытка оптимизировать обобщенно - не имеет смысла. На х86 есть sseшный/фпушный корень - в него и будет транслироваться. Быстрее него, такими методами, ничего написать не реально. Даже такой же.

Т.е. самописный 2корень будет всегда на х86 тормазить. Но вот на чем-то, где такого «железячного» корня нету - будет всё наоборот. Уже будет тормазить дефолтный корень, а самописный можно написать и быстрее.

Поэтому уже давно такой подход для реально оптимизации не используют. Нужна оптимизация? Пишется реализация под каждый нужный таргет.

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

В данном случае обход сетки крайне плох. Везде дыры - упирается всё в латенси.

По теме - корни, 2минуты в гугле:

float qsqrt(float x) {
  float r = x;
  uint32_t * i = ((uint32_t *)&r);
  *i = (0xbe6f0000 - *i) >> 1;
  return (x * r * (1.5f - 0.5f * x * r * r));
}

Далее, -ffinite-math-only - позволяет выпилить lm. Заменяет функцию из lm на свою - раскрывается в чистый симдовый/фпушный корень. В этих функциях нет проверки на некорректное значение аргументов.

У меня функция гцц от -ffinite-math-only работает с такой же скоростью как qsqrt(), даже чутка медленней.

Есть ссешный быстрый обратный корень.

float rsqrt_fast(float value) {
  return _mm_rsqrt_ss(_mm_set_ss(value))[0];
}

Т.к. это обратный корень - надо делать как-то так:

uint16_t SmlGradientRadialValue(SmlLine line, SmlPoint target) {
  if((target.x == line.p1.x) && (target.y == line.p1.y))
    return SML_GRAD_MIN;

  float AB = (SML_SQR(line.p1.x - line.p2.x) +
                  SML_SQR(line.p1.y - line.p2.y));
  float AT = (SML_SQR(line.p1.x - target.x) +
                  SML_SQR(line.p1.y - target.y));
  if(AB < AT) return SML_GRAD_MAX;
  return rsqrt_fast(AB / AT) * SML_GRAD_MAX;
}

Точность чутка пляшет, но это у меня раза в полтора-два быстрее предыдущих вариантов.

anonymous
()
Ответ на: комментарий от sambist

7.6 sqrti

Эти корни с циклами просто априори убоги и не могут никоим образом соревноваться в железячными/глибцешными. Тем более так написанные. Первое - следствие первой основной ошибки, а второе второй.

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

6.4 от fsqrt vs 7.6 sqrti

Ваши результаты снова не верны. Соотношение не реалистично.

anonymous
()
Ответ на: комментарий от sambist

Т.е. на выходе целочисленные значения от 0 до 65535.

Всего 2^16 значений? Да вычисли ты их заранее и просто смотри по таблице. Память нынче копеечная, или ты для калькулятора пишешь?

no-such-file ★★★★★
()
Ответ на: комментарий от no-such-file

Да вычисли ты их заранее и просто смотри по таблице.

Уже пробовал. Поиск по таблице занимает 18 секунд на 100 000.

sambist ★★
() автор топика
Ответ на: комментарий от anonymous

Соотношение не реалистично.

Вот сейчас пойду и нафотошоплю что-нибудь более реалистичное.

sambist ★★
() автор топика
Ответ на: комментарий от sambist

Поиск по таблице занимает 18 секунд на 100 000.

ШТОА? Ты действительно на калькуляторе пишешь чтоли? По таблице 2^16, даже если брать поиском нужно максимум 16 итераций. У тебя один лукап в секунду чтоли?

no-such-file ★★★★★
()
Ответ на: комментарий от no-such-file

На самом деле даже 15 итераций, ибо два крайних «указателя» сойдутся на шаг раньше. Функция радиального градиента (представленная где-то тут в треде), использующая sqrti (который по профайлеру занимает 50% ее времени) 100 000 000 итераций (пардон, я нулики в прошлыйраз не дописал) занимает 18 секунд на табличном поиске.

sambist ★★
() автор топика
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.