Calcul rapide de 1/Sqr()

[comtois]

Je viens de lire ce post

http://forum.games-creators.org/showthr ... #post20898

Je me suis dit , chic , 4 fois plus rapide , c'est parfait , je l'utilise bcp en 3D.

J'ai donc codé ça :

Code : Tout sélectionner

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

Mais d'après mes tests , c'est plus lent que Sqr()

C'est là qu'une optimisation en assembleur serait la bienvenue, mais je ne sais pas faire ça , quelqu'un pourrait se pencher sur la question ?
Ou alors il y a une erreur dans mon code ? il y a plus rapide à faire ?
Je ne sais pas comment remplacer les peekl et peekf ,c'est sûrement ça qui ralentit ?

[filperj]

Petit essai (j'ai pas testé la vitesse)

Code : Tout sélectionner

Procedure.f InvSqrt(x.f)
   Protected i.f  ;x est à [esp], i est à [esp+4]
   !mov eax , [esp]
   !mov ebx , 5F3759DFh
   !sar eax , 1
   !sub ebx , eax
   !mov [esp+4] , ebx
   ProcedureReturn i * (1.5 - 0.5 * x * i * i)
EndProcedure

For i=1 To 20
   Debug StrF(InvSqrt(i)) + " : " + StrF(1/Sqr(i))
Next i

Ceci dit, "1/Sqr(x)" ne génère pas d'appel de proc, mais directement des instructions FPU (dont FSQRT, car le FPU sait faire une racine), alors doit falloir s'accrocher pour faire mieux!
Mais c'est surement possible d'optimiser pour une tâche donnée...

[Fred]

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).

[comtois]

Merci filperj , intéressant , il faudra que je m'exerce à ce genre de manipulations aussi.

Fred ,merci pour les explications , comme ça c'est clair que notre Sqr() est bien le plus rapide

[Dr. Dri]

bah tant qu'à faire autant quelle devienne native de PB ^^

Dri

[Fred]

Au lieu de 1/Sqr() ? Je vois pas trop l'interet la.

[filperj]

Je crois que le docteur a lu un peu en diagonale

J'avais rien à faire et y fait pas beau, alors j'ai essayé cette méthode avec les instructions SSE2, pour calculer tout un paquet de racines à la fois, mais je n'obtiens rien de mieux qu'avec l'instruction SSE "sqrtps" toute simple toute bête:

Code : Tout sélectionner

Procedure.f diffmax(*src1.float, *src2.float, longueur)
   max.f = 0.0
   While longueur
      diff.f = Abs(*src1\f - *src2\f)
      If diff > max
         max = diff
      EndIf
      *src1 + 4
      *src2 + 4
      longueur - 1
   Wend
   ProcedureReturn max
EndProcedure


Procedure SQR_FPU(source, destination, longueur)
   !mov esi , [esp]
   !mov edi , [esp+4]
   !mov eax , [esp+8]
   !mov ebx , 4
   !test eax , eax
   !jz .Fin
   !.Bouclette:
      !fld dword[esi]
      !fsqrt
      !fstp dword[edi]
      !add esi , ebx
      !add edi , ebx
      !dec eax
   !jnz .Bouclette
   !.Fin:
EndProcedure


; Attention: le addresses source et destination doivent être alignées
; sur 16 octets, et la longueur(cad le nombre de flottants 32bits à traiter)
; doit être multiple de 4
Procedure SQR_SSE_simple(source, destination, longueur)
   !mov esi , [esp]
   !mov edi , [esp+4]
   !mov eax , [esp+8]
   !mov ebx , 16
   !shr eax , 2
   !jz .Fin
   !.Bouclette:
      !movaps xmm0 , [esi]
      !sqrtps xmm1 , xmm0
      !movaps [edi] , xmm1
      !add esi , ebx
      !add edi , ebx
      !dec eax
   !jnz .Bouclette
   !.Fin:
EndProcedure

; Attention: le addresses source et destination doivent être alignées
; sur 16 octets, et la longueur(cad le nombre de flottants 32bits à traiter)
; doit être multiple de 4
Procedure SQR_SSE2_plus_truc(source, destination, longueur)
   !mov eax , 0.5
   !push eax
   !push eax
   !push eax
   !push eax
   !movups xmm1 , [esp]
   !_0p5 equ xmm1
   !mov eax , 5F3759DFh
   !mov [esp] , eax
   !mov [esp+4] , eax
   !mov [esp+8] , eax
   !mov [esp+12] , eax
   !movups xmm2 , [esp]
   !magic equ xmm2
   !mov eax , 1.5
   !mov [esp] , eax
   !mov [esp+4] , eax
   !mov [esp+8] , eax
   !mov [esp+12] , eax
   !movups xmm3 , [esp]
   !_1p5 equ xmm3
   !mov eax , 1.0
   !mov [esp] , eax
   !mov [esp+4] , eax
   !mov [esp+8] , eax
   !mov [esp+12] , eax
   !movups xmm4 , [esp]
   !_1p0 equ xmm4
   !add esp , 16
   !mov esi , [esp]
   !mov edi , [esp+4]
   !mov eax , [esp+8]
   !mov ebx , 16
   !shr eax , 2
   !jz .Fin
   !.Bouclette:
      !movaps xmm0 , [esi]
      !movaps xmm5 , _0p5
      !mulps xmm5 , xmm0
      !psrld xmm0 , 1
      !movaps xmm6 , magic
      !psubd xmm6 , xmm0
      !mulps xmm5 , xmm6
      !mulps xmm5 , xmm6
      !movaps xmm7 , _1p5
      !subps xmm7 , xmm5
      !mulps xmm7 , xmm6
      !movaps xmm0 , _1p0
      !divps xmm0 , xmm7
      !movaps [edi] , xmm0
      !add esi , ebx
      !add edi , ebx
      !dec eax
   !jnz .Bouclette
   !.Fin:
EndProcedure

Structure AlignedMemBlock
   AlignedAddr.l
   AllocAddr.l
EndStructure

Procedure.l AllocateAlignedMemBlock(*BlockInfo.AlignedMemBlock, size, Alignment)
   r = AllocateMemory(size + alignment)
   If r
      *BlockInfo\allocaddr = r
      *BlockInfo\alignedaddr = r + alignment - r % alignment
     Else
      *BlockInfo\allocaddr = 0
      *BlockInfo\alignedaddr = 0
   EndIf
   ProcedureReturn r
EndProcedure


DefType.AlignedMemBlock Source, DestFPU, DestSSE_simple, DestSSE2_plus_truc

#TestLen = 1000000
AllocateAlignedMemBlock(@Source, #TestLen * 16, 16)
AllocateAlignedMemBlock(@DestFPU, #TestLen * 16, 16)
AllocateAlignedMemBlock(@DestSSE_simple, #TestLen * 16, 16)
AllocateAlignedMemBlock(@DestSSE2_plus_truc, #TestLen * 16, 16)


*initieur.float = Source\alignedaddr
For i = 1 To #testlen * 4
   *initieur\f = Random(10000)
   *initieur + 4
Next

For i = 1 To 30

   depchrono = ElapsedMilliseconds()
   SQR_FPU(source\alignedaddr, DestFPU\alignedaddr, #testlen * 4)
   finchrono = ElapsedMilliseconds()
   chronoFPU + finchrono - depchrono
   
   depchrono = ElapsedMilliseconds()
   SQR_SSE_simple(source\alignedaddr, DestSSE_simple\alignedaddr, #testlen * 4)
   finchrono = ElapsedMilliseconds()
   chronoSSE_simple + finchrono - depchrono
   
   depchrono = ElapsedMilliseconds()
   SQR_SSE_simple(source\alignedaddr, DestSSE2_plus_truc\alignedaddr, #testlen * 4)
   finchrono = ElapsedMilliseconds()
   chronoSSE2_plus_truc + finchrono - depchrono

Next

MessageRequester("", "FPU: " + Str(chronoFPU) + #CR$ + "SSE simple: " + Str(chronoSSE_simple) + #CR$ + "SSE2 + truc: " + Str(chronoSSE2_plus_truc))

MessageRequester("différence max FPU-SSE simple", StrF(diffmax(DestFPU\alignedaddr, DestSSE_simple\alignedaddr, #testlen * 4)))

MessageRequester("différence max FPU-SSE2 + truc", StrF(diffmax(DestFPU\alignedaddr, DestSSE2_plus_truc\alignedaddr, #testlen * 4)))

Pour calculer un paquet d'inverses de la racine carrée, l'écart FPU-SSE se creuse de manière intéressante(tout simplement parceque SSE a une instruction exprès pour ça), mais je n'arrive toujours pas à tirer parti de la "méthode magique":

Code : Tout sélectionner

Procedure.f diffmax(*src1.float, *src2.float, longueur)
   max.f = 0.0
   While longueur
      diff.f = Abs(*src1\f - *src2\f)
      If diff > max
         max = diff
      EndIf
      *src1 + 4
      *src2 + 4
      longueur - 1
   Wend
   ProcedureReturn max
EndProcedure


Procedure INVSQR_FPU(source, destination, longueur)
   !mov esi , [esp]
   !mov edi , [esp+4]
   !mov eax , [esp+8]
   !mov ebx , 4
   !test eax , eax
   !jz .Fin
   !.Bouclette:
      !fld1
      !fld dword[esi]
      !fsqrt
      !fdivp st1,st0
      !fstp dword[edi]
      !add esi , ebx
      !add edi , ebx
      !dec eax
   !jnz .Bouclette
   !.Fin:
EndProcedure


; Attention: le addresses source et destination doivent être alignées
; sur 16 octets, et la longueur(cad le nombre de flottants 32bits à traiter)
; doit être multiple de 4
Procedure INVSQR_SSE_simple(source, destination, longueur)
   !mov esi , [esp]
   !mov edi , [esp+4]
   !mov eax , [esp+8]
   !mov ebx , 16
   !shr eax , 2
   !jz .Fin
   !.Bouclette:
      !movaps xmm0 , [esi]
      !rsqrtps xmm1 , xmm0
      !movaps [edi] , xmm1
      !add esi , ebx
      !add edi , ebx
      !dec eax
   !jnz .Bouclette
   !.Fin:
EndProcedure

; Attention: le addresses source et destination doivent être alignées
; sur 16 octets, et la longueur(cad le nombre de flottants 32bits à traiter)
; doit être multiple de 4
Procedure INVSQR_SSE2_plus_truc(source, destination, longueur)
   !mov eax , 0.5
   !push eax
   !push eax
   !push eax
   !push eax
   !movups xmm1 , [esp]
   !_0p5 equ xmm1
   !mov eax , 5F3759DFh
   !mov [esp] , eax
   !mov [esp+4] , eax
   !mov [esp+8] , eax
   !mov [esp+12] , eax
   !movups xmm2 , [esp]
   !magic equ xmm2
   !mov eax , 1.5
   !mov [esp] , eax
   !mov [esp+4] , eax
   !mov [esp+8] , eax
   !mov [esp+12] , eax
   !movups xmm3 , [esp]
   !_1p5 equ xmm3
   !add esp , 16
   !mov esi , [esp]
   !mov edi , [esp+4]
   !mov eax , [esp+8]
   !mov ebx , 16
   !shr eax , 2
   !jz .Fin
   !.Bouclette:
      !movaps xmm0 , [esi]
      !movaps xmm5 , _0p5
      !mulps xmm5 , xmm0
      !psrld xmm0 , 1
      !movaps xmm6 , magic
      !psubd xmm6 , xmm0
      !mulps xmm5 , xmm6
      !mulps xmm5 , xmm6
      !movaps xmm7 , _1p5
      !subps xmm7 , xmm5
      !mulps xmm7 , xmm6
      !movaps [edi] , xmm7
      !add esi , ebx
      !add edi , ebx
      !dec eax
   !jnz .Bouclette
   !.Fin:
EndProcedure

Structure AlignedMemBlock
   AlignedAddr.l
   AllocAddr.l
EndStructure

Procedure.l AllocateAlignedMemBlock(*BlockInfo.AlignedMemBlock, size, Alignment)
   r = AllocateMemory(size + alignment)
   If r
      *BlockInfo\allocaddr = r
      *BlockInfo\alignedaddr = r + alignment - r % alignment
     Else
      *BlockInfo\allocaddr = 0
      *BlockInfo\alignedaddr = 0
   EndIf
   ProcedureReturn r
EndProcedure


DefType.AlignedMemBlock Source, DestFPU, DestSSE_simple, DestSSE2_plus_truc

#TestLen = 1000000
AllocateAlignedMemBlock(@Source, #TestLen * 16, 16)
AllocateAlignedMemBlock(@DestFPU, #TestLen * 16, 16)
AllocateAlignedMemBlock(@DestSSE_simple, #TestLen * 16, 16)
AllocateAlignedMemBlock(@DestSSE2_plus_truc, #TestLen * 16, 16)


*initieur.float = Source\alignedaddr
For i = 1 To #testlen * 4
   *initieur\f = Random(10000)
   *initieur + 4
Next

For i = 1 To 30

   depchrono = ElapsedMilliseconds()
   INVSQR_FPU(source\alignedaddr, DestFPU\alignedaddr, #testlen * 4)
   finchrono = ElapsedMilliseconds()
   chronoFPU + finchrono - depchrono
   
   depchrono = ElapsedMilliseconds()
   INVSQR_SSE_simple(source\alignedaddr, DestSSE_simple\alignedaddr, #testlen * 4)
   finchrono = ElapsedMilliseconds()
   chronoSSE_simple + finchrono - depchrono
   
   depchrono = ElapsedMilliseconds()
   INVSQR_SSE_simple(source\alignedaddr, DestSSE2_plus_truc\alignedaddr, #testlen * 4)
   finchrono = ElapsedMilliseconds()
   chronoSSE2_plus_truc + finchrono - depchrono

Next

MessageRequester("", "FPU: " + Str(chronoFPU) + #CR$ + "SSE simple: " + Str(chronoSSE_simple) + #CR$ + "SSE2 + truc: " + Str(chronoSSE2_plus_truc))

MessageRequester("différence max FPU-SSE simple", StrF(diffmax(DestFPU\alignedaddr, DestSSE_simple\alignedaddr, #testlen * 4)))

MessageRequester("différence max FPU-SSE2 + truc", StrF(diffmax(DestFPU\alignedaddr, DestSSE2_plus_truc\alignedaddr, #testlen * 4)))

Bon, c'était toujours plus rigolo que de regarder tomber les gouttes

Ce qui pourrait être intéressant - mais un poil plus couillu - ce serai d'intercaler des instructions CPU qui font autrechose en parallèle...
mais là j'ai peur que ça dépasse mes capacités de concentration!

[djes]

J'ai retrouvé une vieille routine dont je me servais sur Amiga en Blitz. Si j'ai le temps, je la convertirai; en attendant je la laisse là. Je ne sais pas si elle serait plus rapide.

Code : Tout sélectionner

;D'abord la création de la table
  For ii.f=0 To 2048                        ; create square root lookuptable
    vv.f=Sqr(ii/4096)*65535
    Poke.w ?sqrlup+ii*2,Int(vv)
  Next 

;Ensuite la procédure proprement dite

;*****************************************************************
;* Andrews' quick square root routine
Function.q qsqr{a.q}
  LEA sqrlup,a0
  MOVEQ#18,d2:MOVE.l #2048,d3:MOVEQ#0,d4
  loop:CMP.l d3,d0:BLT done
  LSR.l#1,d0:ROXR#1,d4:LSR.l#1,d0:ROXR#1,d4
  DBRA d2,loop:done:ADD d0,d0
  MOVEM 0(a0,d0),d0/d5:MULU d4,d5:NOT d4:MULU d4,d0
  ADD.l d5,d0:LSR.l d2,d0
  AsmExit
End Function 

[filperj]

Houlà, GROSSE BOURDE sur mon code plus haut!
En plus d'être lent, c'est surtout très boggué.

Leq routines simplement SSE ont l'air correctes, si ça intéresse. Mais le bénéfice est pas énorme, faut vraiement avoir beaucoup de racines à calculer.