function Q_rsqrt de John Carmak !

[Thyphoon]

Bon tout le monde connait John Carmak ! Et en consultant sa fiche sur Wiki pedia
http://fr.wikipedia.org/wiki/John_Carmack
j'ai découvert une de ses fonctions (C) pour calculer l'inverse d'une racine.Cette fonction est pour le moins particulière puisqu'elle ne contient aucune boucle, elle ne fait intervenir qu'une suite de calculs élémentaires! Pour autant, celle-ci fournit des approximations tout à fait acceptables (de l'ordre de 10-3). Elle est même jusqu'à 4 fois plus rapide que la fonction sqrt(float x) et se révèle ainsi parfaite pour un jeu vidéo.
et j'ai donc voulu la transposé (pour m'amuser) en PB

Mais je bloc sur un certain nombre de truc ... voici la function d'origine

Code : Tout sélectionner

float Q_rsqrt( float number ){
    long i;
    float x2, y;
    const float threehalfs = 1.5F;
 
    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;  // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 ); // what the fuck?
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
    // y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
 
    #ifndef Q3_VM
    #ifdef __linux__
      assert( !isnan(y) ); // bk010122 - FPE?
    #endif
    #endif
    Return y;
}

et voilà ce que j'ai réussi a faire .... Comme vous pouvez le remarquer il y a pas mal de chose que je n'ai pas su transposer !
il y a le 1.5F et 0.5F et les deuxième calcul de i et de y ! Si quelqu'un a une idée ...

Code : Tout sélectionner

Procedure Q_rsqrt(number.f)
  Protected i.l,x2.f,y.f
  #threehalfs=1.5
  x2 = number * 0.5
  y=number
  i= ????????????
  y= $5f375a86 - ( i >> 1 ); // what the fuck? better magic key than john Carmak
  y  = ????????
  y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
    // y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
ProceduRereturn y
EndProcedure

enfin voilà je voulais testé ! je me dis qu'en regardant dans le code source de Quake3 on pourrait trouver quelques bonnes fonction a transposer ..

[zaphod]

Salut,
Je sais le faire en Emergence basic et en Blitz max mais pas en PB (honte à moi).
Je te donne quand même ce que j'ai (j'ai commenté) :

Emergence :

Code : Tout sélectionner

OPENCONSOLE
float  x =0.9 'on cherche à calculer sqr(x) 
float xhalf = 0.5 * x 
pointer f2i = &x  'le pointeur f2i pointe sur l'adresse de x
int i =*<int>f2i[0]  'i vaut loword du pointeur f2i
*<int>f2i[0]=0x5f3759df - ( i >> 1) 'le calcul se fait sur loword du pointeur f2i
x=x*(1.5 - xhalf * x * x)
x=x*(1.5 - xhalf * x * x) '2 eme tour pour meilleure précision 
'x=x*(1.5 - xhalf * x * x) '3 tour pour améliorer encore 
print 1/x ' valeur de sqr(x)
WAITCON
CLOSECONSOLE
END

et Blitz max :

Code : Tout sélectionner

Local x:Float = 0.9 'on cherche à calculer sqr(x)

Local xhalf:Float = 0.5 * x
Local f2i:Int Ptr = Int Ptr(Varptr x)  'le pointeur f2i pointe sur l'adresse de x
Local i : Int = f2i[0]  'i vaut loword f2i
f2i[0]=$5f3759df - ( i Shr 1)
	
x=x*(1.5 - xhalf * x * x)
x=x*(1.5 - xhalf * x * x)   '2 eme tour pour meilleure précision 
'x=x*(1.5 - xhalf * x * x)  '3 tour pour améliorer encore 

Print 1/x

Je sais pas si ca peut t'aider... mais bon.

[zaphod]

En PB, (je crois que c'est bon) :

Code : Tout sélectionner

x.f =2 ;on cherche à calculer sqr(x) 
xhalf.f = 0.5 * x 
*f2i = @x  ;le pointer f2i pointe sur l'adresse de x
i.i =PeekI(*f2i)  
PokeI(*f2i,$5f3759df - ( i >> 1)) 
x=x*(1.5 - xhalf * x * x)
x=x*(1.5 - xhalf * x * x) ;2 eme tour pour meilleure précision 
;x=x*(1.5 - xhalf * x * x) ;3 tour pour améliorer encore 
Debug 1/x ;' valeur de sqr(x)

[comtois]

ça me rappelle des trucs

http://www.purebasic.fr/french/viewtopi ... &hilit=sqr

[Thyphoon]

Merci a tous !! je comprends mieux le fonctionnement maintenant que je vois le code en purebasic!

Effectivement comtois j'étais passé a côté de ce sujet

mais au moins j'aurais appris quelques choses !

[Fig]

Je peux me tromper mais je crois que cette fonction n'est pas vraiment du génialissime Mr Carmak, mais plutôt du citoyen Héron d'Alexandrie 1er siècle après Jésus Christ
Évidemment il ne programmait pas aussi bien

[Thyphoon]

tu as sans doute raison !
Mais cette fonction avait fait parler d'elle lorsqu'elle avait été découverte dans le code de Quake ... je crois que c'est là d'ou vient le fait qu'on l'assimile a John Carmak !

[Backup]

@Zaphod :

aucun interet

ta fonction est plus lente que le sqr() de purebasic :

comparatif a lancer en mode debug
ps : si vous obtenez 2 fois 0 millisecondes , ajoutez un zero a a boucle t
For t =1 To 20000(0)

Code : Tout sélectionner

 ;*********************************************************************
TempsDepart = ElapsedMilliseconds()  ; Récupère la valeur actuelle
For t =1 To 200000
x.f =2 ;on cherche à calculer sqr(x)
xhalf.f = 0.5 * x
*f2i = @x  ;le pointer f2i pointe sur l'adresse de x
i.i =PeekI(*f2i) 
PokeI(*f2i,$5F3759DF - ( i >> 1))
x=x*(1.5 - xhalf * x * x)
x=x*(1.5 - xhalf * x * x) ;2 eme tour pour meilleure précision
x=x*(1.5 - xhalf * x * x) ;3 tour pour améliorer encore
a.f= 1/x ;' valeur de sqr(x)
Next t

TempsEcoule = ElapsedMilliseconds()-TempsDepart  
Debug "Temps écoulé : "+Str(TempsEcoule)+" millisecondes" 
Debug a.f
;*********************************************************************

TempsDepart = ElapsedMilliseconds()  ; Récupère la valeur actuelle
For t =1 To 200000
    x.f=2 
    a.f=Sqr(x.f)
Next t
TempsEcoule = ElapsedMilliseconds()-TempsDepart  
Debug "Temps écoulé : "+Str(TempsEcoule)+" millisecondes" 
Debug a.f
;*********************************************************************

[djes]

Dobro> tu peux aussi simplement remplacer tes lignes "debug" par des messagerequester pour faire un test qui ne soit pas en mode de débogage (ça fausse tout, surtout quand les codes à comparer sont aussi différents)

[Backup]

oui aussi

[Le Soldat Inconnu]

y = * ( float * ) &i;

Je penche pour

y = peekf(@i)

et pour i = * ( long * ) &y
je dirais

i = peekl(@y)