Calcul rapide de racines carrées

[MrCool]

Voici une astuce tirée du moteur de Quake 3 (quoiqu'on ne sache pas si John Carmack est l'auteur réel ou pas). Merci à Gawen de l'avoir postée sur le chan :)


Voici une fonction qui calcule 1/sqrt():

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 y0
x = *(float *)&i; // convert bits back to float
x = x * (1.5f - xhalf * x * x); // Newton step, repeating increases accuracy
return x;
}


Evidemment: x * (1 / sqrt(x)) = sqrt(x)

Il est donc facile de récupérer sqrt() à partir de cette fonction.

En gros, cette fonction réalise une interpolation linéaire de la fonction 1/sqrt(). Un PDF contenant l'explication complète... (http://www.math.purdue.edu/~clomont/Math/Papers/2003/InvSqrt.pdf)

Il s'agit donc d'une approximation des racines carrées. L'erreur relative maximum est de 0.00175228.

[Thelvyn]

waaaa ... c cool ca :00000020:

[Comtois]

je viens de faire une version PureBasic , alors je ne sais pas où est mon erreur , mais l'écart est plus important que 0.00175228.
ça ne fait rien , c'est plus rapide, je vais l'adopter , merci pour l'info :)

[EDIT]
Après vérification , ce n'est pas plus rapide que Sqr() , au contraire !! Donc je garde mon Sqr()

Voici le code PureBasic

Procedure.f InvSqrt(x.f)
xhalf.f = 0.5*x;
i = PeekL(@x); // get bits for floating value
i = $5f3759df - (i >> 1); // gives initial guess y0
x = PeekF(@i); // convert bits back to float
x = x * (1.5 - xhalf * x * x); // Newton step, repeating increases accuracy
ProcedureReturn x
EndProcedure

For i=1 To 20
a.f=i*InvSqrt(I)
b.f=Sqr(i)
c.f=b-a
Debug StrF(a,4) + " / " + StrF(b,4) + " = " + StrF(c,4)
Next i

Et les 20 résultats

0.9983 / 1.0000 = 0.0017
1.4139 / 1.4142 = 0.0004
1.7305 / 1.7321 = 0.0015
1.9966 / 2.0000 = 0.0034
2.2357 / 2.2361 = 0.0004
2.4461 / 2.4495 = 0.0034
2.6421 / 2.6458 = 0.0036
2.8277 / 2.8284 = 0.0007
2.9966 / 3.0000 = 0.0034
3.1569 / 3.1623 = 0.0054
3.3114 / 3.3166 = 0.0052
3.4611 / 3.4641 = 0.0030
3.6051 / 3.6056 = 0.0005
3.7410 / 3.7417 = 0.0007
3.8663 / 3.8730 = 0.0067
3.9932 / 4.0000 = 0.0068
4.1206 / 4.1231 = 0.0025
4.2422 / 4.2426 = 0.0004
4.3589 / 4.3589 = 0.0000
4.4714 / 4.4721 = 0.0007

[MrCool]

Bizarre, tu es sur de tes fonctions Peek?

D'apres l'auteur du pdf, c'est cense etre 4x plus rapide que la libm... enfin le sqrt de pure basic est peut etre tres optimise ;)

[Comtois]

yes ça doit être ça Sqr() est hyper optimisé .

Pour les Peek j'en sais rien , mais j'obtiens des résultats cohérents ,enfin l'échelle de grandeur est correcte , et puis je ne sais pas comment faire autrement :)

Est-ce que quelqu'un a fait un test de vitesse en C ? c'est à dire avec le code original ? Et pourrait aussi copier le résultat des 20 premières valeurs ? histoire de comparer , je suis curieux :)

[Comtois]

et ça c'est le code assembleur que génère PureBasic , si quelqu'un s'y connait en assembleur pour proposer des optimisations , je suis preneur :)

;
; PureBasic v3.94 (Windows - x86) generated code
;
; (c) 2005 Fantaisie Software
;
; The header must remain intact for Re-Assembly
;
; StringExtension
; MemoryExtension
;
format MS COFF
;
extrn _ExitProcess@4
extrn _HeapCreate@12
extrn _HeapDestroy@4
;
extrn PB_PeekF
extrn PB_PeekL
extrn _PB_StrF2@8
extrn _memset
extrn _SYS_FreeStructureStrings@8
extrn SYS_FreeArray
extrn SYS_FreeStringStructuredArray
extrn _PB_StringBase
extrn PB_StringBase
extrn SYS_AllocateString
extrn SYS_CopyString
extrn SYS_FastAllocateString
extrn SYS_FastAllocateStringFree
extrn SYS_FreeString
extrn SYS_InitString
extrn SYS_StringEqual
extrn SYS_StringSup
extrn SYS_StringInf
public PB_NullString
;
public _PB_Instance
public _PBV_ExecutableType
public _PB_MemoryBase
public PB_Instance
public PB_MemoryBase
public _PB_EndFunctions
public _PB_DEBUGGER_LineNumber
public _PB_DEBUGGER_IncludedFiles
macro align value { rb (value-1) - ($-_PB_DataSection + value-1) mod value }
macro bssalign value { rb (value-1) - ($-_PB_BSSSection + value-1) mod value }
public _WinMain@16
;
section '.text' code readable executable
;
;
_WinMain@16:
;
PUSH dword I_BSSEnd-I_BSSStart
PUSH dword 0
PUSH dword I_BSSStart
CALL _memset
ADD esp,12
MOV eax,[esp+4]
MOV [_PB_Instance],eax
PUSH dword 0
PUSH dword 4000
PUSH dword 0
CALL _HeapCreate@12
MOV [PB_MemoryBase],eax
CALL SYS_InitString
; :
; Procedure.f InvSqrt(x.f)
JMP _EndProcedure0
_Procedure0:
PUSH ebx
PUSH ecx
PUSH ebp
PUSH esi
PUSH edi
MOV esi,esp
SUB esp,12
MOV eax,esp
MOV edx,eax
ADD edx,12
_ClearLoop0:
MOV dword [eax],0
ADD eax,4
CMP eax,edx
JNE _ClearLoop0
MOV eax, dword [esi+24]
MOV dword [esp+0],eax
; xhalf.f = 0.5*x;
FLD dword [esp]
MOV dword [esp-4],1056964608
FMUL dword [esp-4]
FSTP dword [esp+4]
; i = PeekL(@x); // get bits for floating value
LEA eax,[esp]
CALL PB_PeekL
MOV dword [esp+8],eax
; i = $5f3759df - (i >> 1); // gives initial guess y0
MOV ebx,dword [esp+8]
SAR ebx,1
NEG ebx
ADD ebx,1597463007
MOV dword [esp+8],ebx
; x = PeekF(@i); // convert bits back to float
LEA eax,[esp+8]
CALL PB_PeekF
FSTP dword [esp]
; x = x * (1.5 - xhalf * x * x); // Newton step, repeating increases accuracy
FLD dword [esp]
FLD dword [esp+4]
FMUL dword [esp]
FMUL dword [esp]
FCHS
MOV dword [esp-4],1069547520
FADD dword [esp-4]
FMULP ST1,ST0
FSTP dword [esp]
; ProcedureReturn x
FLD dword [esp]
JMP _EndProcedure1
; EndProcedure
XOR eax,eax
_EndProcedure1:
ADD esp,12
POP edi
POP esi
POP ebp
POP ecx
POP ebx
RET 4
_EndProcedure0:
;
; For i=1 To 20
MOV dword [v_i],1
_For1:
MOV eax,20
CMP eax,dword [v_i]
JL _Next2
; a.f=i*InvSqrt(I)
FILD dword [v_i]
PUSH dword [v_i]
FILD dword [esp]
FSTP dword [esp]
CALL _Procedure0
FMULP ST1,ST0
FSTP dword [v_a]
; b.f=Sqr(i)
FILD dword [v_i]
FSQRT
FSTP dword [v_b]
; c.f=b-a
FLD dword [v_b]
FSUB dword [v_a]
FSTP dword [v_c]
; Debug StrF(a,4) + " / " + StrF(b,4) + " = " + StrF(c,4)
; Next i
_NextContinue2:
INC dword [v_i]
JMP _For1
_Next2:
; IDE Options = PureBasic v3.94 (Windows - x86)
; CursorPosition = 14
; Folding = -
_PB_EOP_NoValue:
PUSH dword 0
_PB_EOP:
CALL _PB_EndFunctions
PUSH dword [PB_MemoryBase]
CALL _HeapDestroy@4
CALL _ExitProcess@4
_PB_EndFunctions:
RET
;
;
section '.data' data readable writeable
;
_PB_DataSection:
_PB_DEBUGGER_LineNumber: dd -1
_PB_DEBUGGER_IncludedFiles: dd 0
_PBV_ExecutableType: dd 0
PB_NullString: db 0
align 4
s_s:
dd 0
dd -1
;
section '.bss' readable writeable
_PB_BSSSection:
bssalign 4
;
I_BSSStart:
_PB_MemoryBase:
PB_MemoryBase: rd 1
_PB_Instance:
PB_Instance: rd 1
;
bssalign 4
PB_DataPointer rd 1
v_a rd 1
v_b rd 1
v_c rd 1
v_i rd 1
bssalign 4
bssalign 4
bssalign 4
bssalign 4
I_BSSEnd:
section '.data' data readable writeable
SYS_EndDataSection:

[Loulou]

Avec les processeurs actuels, et surtout le temps de plus en plus grand que les applis passent sur le GPU, je me demande si ce genre d'algorithmes "bidouille" sont encore utiles.

Bref, pour ceux que ça intéresse, nVidia a fait une compil plus ou moins intéressante : Fast Math Routines (http://developer.nvidia.com/object/fast_math_routines.html)

Il me semble que Ogre intègre également un module de maths optimisé.

Sur les forums de Gamedev / Devmaster / Flipcode, il y a aussi des posts avec ce genre de challenge, et les posteurs qui s'y donnaient de bon coeur pour trouver la feinte la plus rapide. Je crois même que quelqu'un en a fait une compil, mais impossible de la retrouver.

[LXS]

Et via le développement limité on arriverait pas à calculer ça encore plus vite?

[Atréides]

Hmm.. il me semble que j'ai lu cet algo dans un bouquin vieux de ... heu ... au moins 30 ans.
Je vais y jeter un oeil, si je le retrouve :)

[grob1212]

Serait-il possible de comparer sur une plus grande échelle de nombres et d'établir la précision de la méthode par rapport au sqrf standard ?

[Loulou]

Et via le développement limité on arriverait pas à calculer ça encore plus vite?
A mon avis c'est justement la méthode "standard" utilisée, donc pas la plus rapide.

Il y a aussi autre chose : il ne faut pas oublier que les CPU intègrent de plus en plus de calculs sur les flottants directement, peut-être même une racine carrée. Mais ça c'est à vérifier, parce que je n'en sais pas plus que ça :)

[Comtois]

J'ai posé la question sur le forum de PureBasic Fr.

Forum PureBasic (http://purebasic.hmt-forum.com/viewtopic.php?t=4113)

Et voici ce que dit Fred (l'auteur de PureBasic)

Un des problemes vient de ce qui est ajouté autour de la procedure pour la gerer (ClearLoop, push/pop des registres). Mais bon, c'est sur que FSQRT est plutot rapide, surtout quand on considere que 5 multiplications d'entiers sont quand meme necessaires. D'apres ma doc, FSQRT prend environ 80 cycles, et un seul MULI prends en moyenne 30 cycles. Un DIVI prend 40 cycles. Donc grosso modo, vu le nombre reduite d'instruction, 1/Sqr() devrait etre toujours plus rapide (a tester quand meme sur differents proc, les intels/amd ne reagissant pas du tout pareils).

Conclusion , je garde mon Sqr() :)

[grob1212]

En Intel SSE, il y a des instructions pour calculer les racines carrés et chez AMD 3D Now il me semble.

[LXS]

Sqrt existe sur certains GPU oui, enfin c'est que j'en ai compris en lisant ça (http://www.vis.uni-stuttgart.de/vis03_tutorial/weiskopf_intro_hw.pdf).