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

Точнее, корень можно только один раз вычислять: при нормализации в конце извлекать из частного AT^2/AB^2.

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

У вас вся формула герона на целочисленной арифметике будет постоянно выдавать либо 0, либо 1, ибо там, где коэффициенты будут хоть немного близкими их отношение сразу вместо 0.3 или 0.5 сползет к 0 или 1, плавали, знаем.

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

зависит от требуемой точности. Например определить функцию расстояния как большая из dX,dY плюс 2/5 меньшей - ну и можно уточнять по потребностям.

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

Лучше даже (точнее) итерационно считать корень из (SML_GRAD_MAX^2 * AT^2 / AB^2). Нужно только условие остановки цикла продумать.

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

Нужно только условие остановки цикла продумать.

Элементарное - abs(Xn - Xn-1) > 1

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

Помножь квадрат на 100, ёлы! Результат подели на 10.

pS
()

Есть алгоритм быстрого извлечения обратного квадратного корня. Работает в целых числах. Естественно, точность пострадает.

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

Кармаковский который? У себя реализовал, но оно на float'ах. Причем строго завязано.

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

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

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

Оно. Алгоритм можно допилить, что бы он пользовался только целочисленной арифметикой. Вкратце - представить флоат в виде дроби и отдельно посчитать числитель и знаменатель.

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

Промежуточные вычисления нужно проводить в длинных целых.

Там не будет переполнения, 6553[5..6]^2 укладывается в 2^32

а полученный результат делить на нормализующий коэффициент.

Как получить нормализующий коэффициент? Я тут прикинул пару вариантов - он везде разный, зависимость пока не нашел.

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

Ну она либо не изменится, либо изменится на некоторый коэффициент и будет К*0x5f3759df.

Я о том, что она привязана в формату float из IEEE754.

hateyoufeel ★★★★★
()
Ответ на: комментарий от hateyoufeel
The calculation of y = 1/√x is based on the identity

    \log_2(y) = - \frac{1}{2}\log_2(x)

Using the approximation of the logarithm above, applied to both x and y, the above equation gives:

    I_y \approx \frac{3}{2} L (B - \sigma) - \frac{1}{2} I_x

which is written in the code as

i  = 0x5f3759df - ( i >> 1 );

The first term above is the “magic number”

    \frac{3}{2} L (B - \sigma) = \text{0x5f3759df}

from which it can be inferred σ ≈ 0.0450466. The second term, ½ Ix, is calculated by shifting the bits of Ix one position to the right.

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

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

Вычисли на бумажке sqrt(UNSIGNED_LONG_MAX/(SML_GRAD_MAX^2)).
При умножении на квадрат нормализующего коэффициента нужно расширяться до unsigned long.

pS
()

Вот ещё некто Николай Гарбуз рассуждает о проблеме квадратных корней.

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

Значит, 1 и есть. Корень из (SML_GRAD_MAX^2 / AB^2 * AT^2) обеспечит нужную точность. Порядок операций именно такой.

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

Корень из (SML_GRAD_MAX^2 / AB^2 * AT^2) обеспечит нужную точность. Порядок операций именно такой.

И чем это отличается от первоначальной задачи?

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

Ещё раз: только для малых величин подкоренного выражения. Для этих случаев можно предусмотреть коррекцию путём нормализации (умножения на квадрат некоего числа, деление результата на это число). Задачу нахождения критического значения подкоренного выражения и вычисления соответствующего коэффициента нормализации возлагаю на ваши плечи.

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

Нет проблем вообще. Проверил лично для малых чисел (от 1 до 9). Алгоритм либо стабилизируется на ближайшем снизу значении, либо колеблется вокруг точного значения корня. Следовательно, для проверки окончания вычислений необходимо проверять условие x(i+1)>=x(i) и возвращать значение x(i).

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

Мда... В общем написал реализацию, проверил, работает. Вот только выигрыш... Кхем. Немного не в ту сторону:

float: 100000x1000 cycles took 9.280323 s
sqrti: 100000x1000 cycles took 42.152074 s
Реализовал так:
uint32_t sqrti(uint32_t value)
{
    uint32_t p;
    uint32_t n = value;
    do
    {
        p = n;
        n = (p + value / p) >> 1;
    }
    while (p > n);

    return n;
}

Все что менял, это:

   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) +
   uint32_t AB = sqrti(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) +
   uint32_t AT = sqrti(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 ★★
() автор топика
Последнее исправление: sambist (всего исправлений: 1)
Ответ на: комментарий от sambist

Два раза корень извлекается. AB и AT надо оставить квадратами. В конце return sqrti(SML_GRAD_MAX * SML_GRAD_MAX / AB * AT);.

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

Проверял, при значениях меньше 1000 - на 2х итерациях, при 1500 примерно 4-5, выше экспотенциально.

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

Вот те варианты вычислений что пробовал:

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

    return n;

//    return n;

    //fast
//       unsigned int m, y, b;
//       m = 0x4000;
//       y = 0;
//       while (m != 0){
//          b = y | m;
//          y = y >> 1;
//          if (value >= b) {
//             value = value - b;
//             y = y | m;
//          }
//          m = m >> 2;
//       }
//       return y;

    // chld
//    uint32_t c = 0;
//    uint32_t d = 1;
//    while (value)
//    {
//        value -= d;
//        c++;
//        d += 2;
//    }
//    return c;

    //root a few iterations
//    uint32_t x;
//    x = (value / 0x3f + 0x3f) >> 1;
//    x = (value / x + x) >> 1;
//    x = (value / x + x) >> 1;
//    return x;
}

Все на двух корнях показывают 42, на одном 30 секунд.

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

Что вообще за железо?

Сейчас x86. Есть задумка потом этот код портировать на что-то более слабое, но там, где можно поднять минимальные иксы.

Почему так хочется отказаться от math.h?

Потому что math.h в этот модуль тащится только из-за этого одного сраного sqrt. Да и очень неплохо было бы нафиг отказаться от float'ов в функции, которая теоретически может вызываться около миллиона раз в секунду.

sambist ★★
() автор топика

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

Кстати, что это за «градиент» такой у тебя? Не очень понимаю, но что будет, если корень с квадратом заменить на модуль?

hey-propagandist
()

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

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