численные методы

[avva]

Исследование одного очень красивого хака (PDF, англ.) - быстрого вычисления функции 1/sqrt(x) (обратный квадратный корень):

float InvSqrt(float x)
{
    float xhalf = 0.5f*x;
    int i = *(int*)&x;       // get bits for floating value
    i = 0x5f3759df - (i>>1); // gives initial guess 
    x = *(float*)&i;         // convert bits back to float
    x = x*(1.5f-xhalf*x*x);  // Newton step, repeating increases accuracy
    return x;
}

Статья объясняет, почему и как это работает, и откуда взялась таинственная магическая константа 0x5f3759df.

Очень красиво.

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

Comments

[ivan-gandhi.livejournal.com]

Kind of obvious, but cool.

[avva.livejournal.com]

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

[ygam.livejournal.com]

Ой.

[dmpogo.livejournal.com]

Well, Newton method is quadratic in convergence near the solution, so if intial guess is good, no suprise first step gives not a bad number, especially since demanded accuracy is not that high.

???

[shkrobius.livejournal.com]

In the end the paper does not explain where the mysterious constant is from, it only tells why it is close to the optimum choice. Perhaps the constant was found by trial and error.

[neo-techno.livejournal.com]

а вот история поиска кто написал эту функцию

http://www.beyond3d.com/articles/fastinvsqrt/

[nm-work.livejournal.com]

супер :)

Недурно. А я, кстати,

[sravniti.com (from livejournal.com)]

при рассмотрении этой функции,впервые рассматривал a piece of source на си, без всяких дурных мыслей о том, как бы ее выразить на любимой яве?

[iratus.livejournal.com]

Таки 3DFX ;-)
насколько я помню, magic number просто подбирали долго...

[faceted-jacinth.livejournal.com]

прикольно

[dimchansky.livejournal.com]

Так Нюьтоном ещё при делении чисел пользуются, например. Если нужно найти u/v, то находят приближение 1/v по ньютону (http://numbers.computation.free.fr/Constants/Algorithms/inverse.html), а потом множат на u. Это при условии, что есть быстрые алгоритмы умножения.

[beyba.livejournal.com]

красота какая, спасибо!

[lordakryl.livejournal.com]

Интересно.
А как Вы находите такие ссылки? Я бы почитал ещё что-нибудь подобное

[ygam.livejournal.com]

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

[avva.livejournal.com]

Эту прислал знакомый, но он ее, кажется, увидел на programming.reddit.com (я это не читаю, нет времени).

[avva.livejournal.com]

Всегда пожалуйста :)

Re: ???

[avva.livejournal.com]

Ага.

[dimchansky.livejournal.com]

Кстати, Кнут во втором томе (4.3.3) приводит алгоритм R, основанный на алгоритме Кука, вычисляющий обратное N-битовое число.
Если v=(0.v₁v₂v₃v₄...)₂, где v₁=1, то начальное приближение z=1/v выбирается таким:
z₀ = 1/4 * [32/(4*v₁ + 2*v₂ + v₃)],
а далее проводятся итерации по Ньютону для k=0,...; 2^k < N:
V_k=(0.v₁v₂v₃...v_2^k+1+3)₂
z_k+1 = 2*z_k - V_k*z_k² + r, где r, удовлетворяющее неравенству, 0<=r<2^-2k+1-1, прибавляется при необходимости, чтобы z_k+1 было кратным 2^-2k+1-1.

А с таблицами - это как?

[lordakryl.livejournal.com]

ок, спасибо

[ygam.livejournal.com]

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

Деление через метод Ньютона - это, например, http://www.ecs.umass.edu/ece/koren/arith/slides/Part8-div-mlt.ppt

[ltwood.livejournal.com]

О подобных вещах есть целая книга -- Hacker's Delight, Henry S. Warren Jr. [http://www.hackersdelight.org/]. Кстати, переведена на русский.