?

Log in

численные методы - Поклонник деепричастий [entries|archive|friends|userinfo]
Anatoly Vorobey

[ website | Website ]
[ userinfo | livejournal userinfo ]
[ archive | journal archive ]

Links
[Links:| English-language weblog ]

численные методы [дек. 2, 2006|03:46 am]
Anatoly Vorobey

Исследование одного очень красивого хака (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:
[User Picture]From: juan_gandhi
2006-12-02 01:59 am
Kind of obvious, but cool.
(Ответить) (Thread)
[User Picture]From: avva
2006-12-02 02:05 am
Наверное, нетривиальное в этом не сам код, а изначальная идея того, что можно будет быстро (две операции с целыми числами) подобрать хорошее начальное приближение, и после этого всего одной итерации хватит, чтобы получить очень высокую точность. Мне, скажем, это совсем не очевидно, даже и в голову не пришло бы это пробовать. Но я совсем не имел дело на практике с численными методами и у меня нет интуиции насчет скорости схождения приближений.
(Ответить) (Parent) (Thread)
From: dmpogo
2006-12-02 02:43 am
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.
(Ответить) (Parent) (Thread)
[User Picture]From: ygam
2006-12-02 02:39 am
Ой.
(Ответить) (Thread)
[User Picture]From: shkrobius
2006-12-02 03:34 am

???

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.
(Ответить) (Thread)
[User Picture]From: avva
2006-12-02 10:39 pm

Re: ???

Ага.
(Ответить) (Parent) (Thread)
[User Picture]From: neo_techno
2006-12-02 05:32 am
а вот история поиска кто написал эту функцию

http://www.beyond3d.com/articles/fastinvsqrt/
(Ответить) (Thread)
[User Picture]From: nm_work
2006-12-02 06:39 am
супер :)
(Ответить) (Parent) (Thread)
[User Picture]From: faceted_jacinth
2006-12-02 02:00 pm
прикольно
(Ответить) (Parent) (Thread)
From: ext_21754
2006-12-02 09:19 am

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

при рассмотрении этой функции,впервые рассматривал a piece of source на си, без всяких дурных мыслей о том, как бы ее выразить на любимой яве?
(Ответить) (Thread)
[User Picture]From: iratus
2006-12-02 09:43 am
Таки 3DFX ;-)
насколько я помню, magic number просто подбирали долго...
(Ответить) (Thread)
[User Picture]From: dimchansky
2006-12-02 05:23 pm
Так Нюьтоном ещё при делении чисел пользуются, например. Если нужно найти u/v, то находят приближение 1/v по ньютону, а потом множат на u. Это при условии, что есть быстрые алгоритмы умножения.
(Ответить) (Thread)
[User Picture]From: ygam
2006-12-02 10:16 pm
Именно так это обычно и реализуется в железе. Только начальные биты смотрят по таблице, а не находят путем хака.
(Ответить) (Parent) (Thread)
[User Picture]From: dimchansky
2006-12-03 02:51 am
Кстати, Кнут во втором томе (4.3.3) приводит алгоритм R, основанный на алгоритме Кука, вычисляющий обратное N-битовое число.
Если v=(0.v1v2v3v4...)2, где v1=1, то начальное приближение z=1/v выбирается таким:
z0 = 1/4 * [32/(4*v1 + 2*v2 + v3)],
а далее проводятся итерации по Ньютону для k=0,...; 2k < N:
Vk=(0.v1v2v3...v2k+1+3)2
zk+1 = 2*zk - Vk*zk2 + r, где r, удовлетворяющее неравенству, 0<=r<2-2k+1-1, прибавляется при необходимости, чтобы zk+1 было кратным 2-2k+1-1.

А с таблицами - это как?
(Ответить) (Parent) (Thread)
[User Picture]From: ygam
2006-12-03 03:00 am
Таблица - это, например, программируемая логика, которая находит первые несколько бит частного по первым нескольким битам делимого.

Деление через метод Ньютона - это, например, http://www.ecs.umass.edu/ece/koren/arith/slides/Part8-div-mlt.ppt
(Ответить) (Parent) (Thread)
[User Picture]From: beyba
2006-12-02 06:29 pm
красота какая, спасибо!
(Ответить) (Thread)
[User Picture]From: avva
2006-12-02 10:39 pm
Всегда пожалуйста :)
(Ответить) (Parent) (Thread)
[User Picture]From: lordakryl
2006-12-02 06:45 pm
Интересно.
А как Вы находите такие ссылки? Я бы почитал ещё что-нибудь подобное
(Ответить) (Thread)
[User Picture]From: avva
2006-12-02 10:38 pm
Эту прислал знакомый, но он ее, кажется, увидел на programming.reddit.com (я это не читаю, нет времени).
(Ответить) (Parent) (Thread)
[User Picture]From: lordakryl
2006-12-03 02:59 am
ок, спасибо
(Ответить) (Parent) (Thread)
[User Picture]From: ltwood
2006-12-04 03:57 pm
О подобных вещах есть целая книга -- Hacker's Delight, Henry S. Warren Jr. [http://www.hackersdelight.org/]. Кстати, переведена на русский.
(Ответить) (Parent) (Thread)