conversion de routine PB en ASM - Mandelbrot

Pour discuter de l'assembleur
Avatar de l’utilisateur
Guillot
Messages : 532
Inscription : jeu. 25/juin/2015 16:18

conversion de routine PB en ASM - Mandelbrot

Message par Guillot »

est ce que quelqu'un pourrait me convertir la routine suivante ?
(je pense en particulier à Manababel, très actif en ce moment)

Code : Tout sélectionner

Procedure madelbrot(a.d,b.d)
    Protected.d aa,bb,c,d
    
    c = a: d = b
    For i = 0 To imax - 2
        aa = a * a
        bb = b * b
        If aa + bb > 4:Break:EndIf
        b = 2 * a * b + d
        a = aa - bb + c
    Next
    ProcedureReturn  i
EndProcedure
merci d'avance
manababel
Messages : 136
Inscription : jeu. 14/mai/2020 7:40

Re: conversion de routine PB en ASM - Mandelbrot

Message par manababel »

voici un petit programme avec la procédure écrit en asm

Code : Tout sélectionner

Global imax=255
Global lg=600
Global ht=400

  
Procedure madelbrot(a.d,b.d)

  cp.q=imax
  fo.d=4.0
  !movsd xmm7,[p.v_fo]
  !movsd xmm0,[p.v_a] ; a
  !movsd xmm1,[p.v_b] ; b
  !movupd xmm2,xmm0 ; c=a
  !movupd xmm3,xmm1 ; d=b 
  !xor rax,rax
  !boucle:
    
  !movupd xmm4,xmm0 
  !movupd xmm5,xmm1
  !mulpd xmm4,xmm4 ; aa = a * a
  !mulpd xmm5,xmm5 ; bb = b * b

  ;If (aa + bb) > 4.0:Goto fin:EndIf
  !movupd xmm6,xmm4
  !addpd xmm6,xmm5
  !VCMPGTPD xmm6,xmm6,xmm7
  !movq rdx,xmm6
  !cmp rdx,0
  !jne fin
  
  ;x1 = 2 * x0 * x1 + x3
  !mulpd xmm1,xmm0 ; b*a
  !addpd xmm1,xmm1 ; (a*b)*2
  !addpd xmm1,xmm3 ; (a*b*2)+d
  
  ;x0 = x4 - x5 + c
  !movupd xmm0,xmm4 ; a = aa
  !subpd xmm0,xmm5 ; a = aa - bb
  !addpd xmm0,xmm2 ;  a = aa - bb + c
  
  !inc rax
  !cmp rax,[p.v_cp]
  !jb boucle

  !fin:

  ProcedureReturn ; rax

EndProcedure


If OpenWindow(0, 0, 0, lg+20, ht+40, "Exemple...", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
  
  img = CreateImage(1, lg,ht,32)
  StartDrawing(ImageOutput(1))
  buf = DrawingBuffer()
  StopDrawing()
  
  ht2=ht/2
  lg2=lg/2
  lg4=lg2+lg/4
  Repeat
    Event = WindowEvent()
    t=ElapsedMilliseconds()
    For y=1 To ht-1
      py.d=(y-ht2)/ht2
      For x=1 To lg-1
        px.d=(x-lg4)/lg2
        var=madelbrot(px,py)*$10101
        PokeL((buf+((y*lg)+x)<<2),var)
      Next
    Next
    
    t=ElapsedMilliseconds()-t

    StartDrawing(WindowOutput(0))
    DrawImage(img,10,30)
    DrawText(5,5,Str(t)+"  ")
    StopDrawing()   
    Select Event
      Case #PB_Event_Gadget
        Select EventGadget()     
          Case 1
            CloseWindow(0)
            End     
        EndSelect   
    EndSelect
  Until Event = #PB_Event_CloseWindow
EndIf
manababel
Messages : 136
Inscription : jeu. 14/mai/2020 7:40

Re: conversion de routine PB en ASM - Mandelbrot

Message par manababel »

voici une version plus difficile à mettre en oeuvre
à essayer avec une structure

version corrigée

var = madelbrot(x1,x2,cy) revoie une valeur de 64 bits
(var & $0000ffff) , correspond a la couleur du pixel x1
(var & $ffff0000) , correspond a la couleur du pixel x2

Code : Tout sélectionner

Global lg=800
Global ht=600
Global infinity.f=4
Global limit=255
 


Procedure madelbrot_ok(a.d,b.d)

  cp.q=limit
  fo.d=infinity
  !movsd xmm7,[p.v_fo]
  !movsd xmm0,[p.v_a] ; a
  !movsd xmm1,[p.v_b] ; b
  !movupd xmm2,xmm0 ; c=a
  !movupd xmm3,xmm1 ; d=b
  !xor rax,rax
  !boucle:
   
  !movupd xmm4,xmm0
  !movupd xmm5,xmm1
  !mulpd xmm4,xmm4 ; aa = a * a
  !mulpd xmm5,xmm5 ; bb = b * b

  ;If (aa + bb) > 4.0:Goto fin:EndIf
  !movupd xmm6,xmm4
  !addpd xmm6,xmm5
  ;!VCMPGTPD xmm6,xmm6,xmm7
  !cmpnltsd xmm6,xmm7
  !movq rdx,xmm6
  !cmp rdx,0
  !jne fin
 
  ;x1 = 2 * x0 * x1 + x3
  !mulpd xmm1,xmm0 ; b*a
  !addpd xmm1,xmm1 ; (a*b)*2
  !addpd xmm1,xmm3 ; (a*b*2)+d
 
  ;x0 = x4 - x5 + c
  !movupd xmm0,xmm4 ; a = aa
  !subpd xmm0,xmm5 ; a = aa - bb
  !addpd xmm0,xmm2 ;  a = aa - bb + c
 
  !inc rax
  !cmp rax,[p.v_cp]
  !jb boucle

  !fin:

  ProcedureReturn ; rax

EndProcedure





Procedure.q madelbrot(a1.d,a2.d,b.d)
  
  rep.q=0
  cp.q=limit
  fo.d=infinity
  !VPBROADCASTQ xmm7,[p.v_fo]
  !MOVLPD xmm0,[p.v_a2] ; a
  !MOVHPD xmm0,[p.v_a1] ; a
  !VPBROADCASTQ xmm1,[p.v_b] ; b
  !movupd xmm2,xmm0 ; c=a
  !movupd xmm3,xmm1 ; d=b
  !xor r10,r10
  !xor r11,r11
  !xor r15,r15
  
  !xor rax,rax
  !boucle:
   
  !movupd xmm4,xmm0
  !movupd xmm5,xmm1
  !mulpd xmm4,xmm4 ; aa = a * a
  !mulpd xmm5,xmm5 ; bb = b * b

  ;If (aa + bb) > 4.0:Goto fin:EndIf
  !movupd xmm6,xmm4
  !addpd xmm6,xmm5
  !VCMPGtPD xmm8,xmm6,xmm7
  !PEXTRQ r8,xmm8,0
  !PEXTRQ r9,xmm8,1

  !cmp r8,0
  !jz saut1
    !cmp r10,0
    !jnz saut1
      !inc r10
      !mov rdx,rax
      !shl rdx,16
      !or r15,rdx
      !cmp r11,1
      !je fin
       
  !saut1:
  !cmp r9,0
  !jz saut2
    !cmp r11,0
    !jnz saut2
      !inc r11
      !or r15,rax
      !cmp r10,1
      !je fin
    
  !saut2:
  ;x1 = 2 * x0 * x1 + x3
  !mulpd xmm1,xmm0 ; b*a
  !addpd xmm1,xmm1 ; (a*b)*2
  !addpd xmm1,xmm3 ; (a*b*2)+d
 
  ;x0 = x4 - x5 + c
  !movupd xmm0,xmm4 ; a = aa
  !subpd xmm0,xmm5 ; a = aa - bb
  !addpd xmm0,xmm2 ;  a = aa - bb + c
 
  !inc rax
  !cmp rax,[p.v_cp]
  !jb boucle
  
  !cmp r10,0
  !jne saut3
    !mov rcx,rax
    !shl rcx,16
    !or r15,rcx
  !saut3:
    
  !cmp r11,0
  !jne fin
    !mov rcx,255
    !or r15,rcx
  !fin:
  !mov rax,r15
  

  ProcedureReturn 

EndProcedure

If OpenWindow(0, 0, 0, lg+20, ht+40, "Exemple...", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
 
  img = CreateImage(1, lg,ht,32)
  StartDrawing(ImageOutput(1))
  buf = DrawingBuffer()
  StopDrawing()
  
  yside.f = -2.5
  xside.f = 2.5
  tside.f = 1.25
  lside.f = -2
  xscale.f=xside/lg
  yscale.f=yside/ht
  cx1.f=0
  cx2.f=0
  cy.f=0

    t=ElapsedMilliseconds()
    For y=0 To ht-1
      cy = y * yscale + tside
      For x=0 To lg-1 Step 2
        cx1 = x * xscale + lside 
        cx2 = (x+1) * xscale + lside 
        var=madelbrot(cx1,cx2,cy)
        var1=var & $ffff
        var2=(var & $ffff0000)>>16
        PokeL((buf+((y*lg)+x)<<2),var1*$10101)
        PokeL((buf+((y*lg)+x+1)<<2),var2*$10101)
      Next
    Next
    t=ElapsedMilliseconds()-t

    StartDrawing(WindowOutput(0))
    DrawImage(img,10,30)
    DrawText(5,5,Str(t)+"  ")
    StopDrawing()   
    
    
Repeat
  Event = WindowEvent()
  Until Event = #PB_Event_CloseWindow
EndIf
Dernière modification par manababel le dim. 18/oct./2020 15:58, modifié 1 fois.
Avatar de l’utilisateur
Guillot
Messages : 532
Inscription : jeu. 25/juin/2015 16:18

Re: conversion de routine PB en ASM - Mandelbrot

Message par Guillot »

merci

au moins 2 X plus rapide

j'ai pas compris ta 2eme version

je te met vieux code pour naviguer (souris + roulette)
un peu naze, je corrigerai ça lundi

Code : Tout sélectionner

Global CBprof=200,imax,ex=1200,ey=800

Global Dim pal.l(4096)
Global Dim Bmp.l(ey - 1, ex - 1)
Global Dim t.w(ey, ex - 1)
Global Dim cx.d(ex - 1)
Global Dim cy.d(ey - 1)

Global ndt=CountCPUs(#PB_System_ProcessCPUs )
Global Dim Thread(ndt+1)
Global Event.i,z.d, fx.d, fy.d,zoom.d

Procedure ColorBlend(color1.l, color2.l, blend.f)
    Protected r.w,g.w,b.w,a.w
    r=  Red(color1) + (Red(color2)     - Red(color1)) * blend
    g=Green(color1) + (Green(color2) - Green(color1)) * blend
    b= Blue(color1) + (Blue(color2) -   Blue(color1)) * blend
    a=Alpha(color1) + (Alpha(color2) - Alpha(color1)) * blend
    ProcedureReturn  RGBA(r,g,b,a)
EndProcedure

Procedure InitPalette()
    r.b:g.b:b.b
    Dim co.l(64)
    ;co = Array(&HFF0000, &HFF00FF, &HFF&, &HFFFF&, &H4400&, &HFF0000, &HFF00FF, &HFF&, &HFFFF&, &H4400&, &HFF0000, &HFF00FF, &HFF&, &HFFFF&, &H4400&, &HFF0000, &HFF00FF, &HFF&, &HFFFF&, &H4400&)
    For i = 0 To 64: co(i) = Random( $FFFFFF): Next
    cl = 8
    For j = 0 To 63
        For i = 0 To cl - 1
            ci = j * cl + i
            pal(ci) = ColorBlend(co(j), co(j + 1), i/ cl)
        Next
    Next
    imax = CBprof: pal(imax) = 0
    ;For i = 0 To 255: CoRGB pal(i), b, g, r: Line (i * 2, 0)-(i * 2 + 1, 10), RGB(b, g, r), BF: Next
EndProcedure

Procedure madelbrot(a.d,b.d)

  cp.q=imax
  fo.d=4.0
  !movsd xmm7,[p.v_fo]
  !movsd xmm0,[p.v_a] ; a
  !movsd xmm1,[p.v_b] ; b
  !movupd xmm2,xmm0 ; c=a
  !movupd xmm3,xmm1 ; d=b
  !xor rax,rax
  !boucle:
   
  !movupd xmm4,xmm0
  !movupd xmm5,xmm1
  !mulpd xmm4,xmm4 ; aa = a * a
  !mulpd xmm5,xmm5 ; bb = b * b

  ;If (aa + bb) > 4.0:Goto fin:EndIf
  !movupd xmm6,xmm4
  !addpd xmm6,xmm5
  !VCMPGTPD xmm6,xmm6,xmm7
  !movq rdx,xmm6
  !cmp rdx,0
  !jne fin
 
  ;x1 = 2 * x0 * x1 + x3
  !mulpd xmm1,xmm0 ; b*a
  !addpd xmm1,xmm1 ; (a*b)*2
  !addpd xmm1,xmm3 ; (a*b*2)+d
 
  ;x0 = x4 - x5 + c
  !movupd xmm0,xmm4 ; a = aa
  !subpd xmm0,xmm5 ; a = aa - bb
  !addpd xmm0,xmm2 ;  a = aa - bb + c
 
  !inc rax
  !cmp rax,[p.v_cp]
  !jb boucle

  !fin:

  ProcedureReturn ; rax

EndProcedure

Procedure ___madelbrot(a.d,b.d)
    Protected.d aa,bb,c,d
    
    c = a: d = b
    For i = 0 To imax
        aa = a * a
        bb = b * b
        If aa + bb > 4:Break:EndIf
        b = 2 * a * b + d
        a = aa - bb + c
    Next
    ProcedureReturn  i
EndProcedure

Procedure fractale(num)
    Protected i,j,jdeb,jfin  ,px,py
    jdeb=ey/ndt*num
    jfin=ey/ndt*(num+1)
    
    For i = 0 To ex - 1   : cx(i) = zoom * ((2 * i) / ex -1) * ex / ey + fx: Next
    For j = jdeb To jfin-1: cy(j) = zoom * ((2 * j) / ey -1) + fy: Next
    
    dd=2
    For jj = jdeb/dd To jfin/dd:j=jj*dd 
        For ii = 0 To ex/dd - 1 :i=ii*dd
            t(j,i) = madelbrot((cx(i)), (cy(j)))
        Next
    Next
    While dd > 1
        d = dd / 2: c1 = 0: c2 = 0: c3 = 0
        
        For jj = jdeb/dd+2  To jfin/dd :j=jj*dd-d
            For ii = 2 To ex/dd - 1 :i=ii*dd-d
                pa = t(j - d,i - d)
                pb = t(j + d,i - d)
                pc = t( j + d,i + d)
                pd = t(j - d,i + d)
                If pa = pb And pa = pc And pa = pd 
                    pe = pa: c1 = c1 + 1 
                Else 
                    pe = madelbrot((cx(i)), (cy(j)))
                EndIf
                pf = t( j,i - dd)
                pg = t(j - dd,i)
                t(j,i) = pe
                If pe = pc And pe = pf And pe = pd 
                    t( j,i + d) = pe: c2 = c2 + 1 
                Else 
                    t(j,i + d) = madelbrot(cx(i + d), cy(j))
                EndIf
                If pe = pc And pe = pg And pe = pb 
                    t(j + d,i) = pe: c3 = c3 + 1 
                Else 
                    t(j + d,i) = madelbrot((cx(i)), (cy(j + d)))
                EndIf
            Next
        Next
        dd = dd / 2
        ;Debug ex * ey, c1; c2; c3
    Wend
    
    For j=jdeb To jfin-1
        For i=0 To ex-1
            bmp(j,i)=pal(t(j,i))
        Next
    Next
;     For j=jdeb To jfin-1
;         For i=0 To ex-1
;             bmp(j,i)=pal(madelbrot((cx(i)), (cy(j + d))))
;         Next
;     Next
    
EndProcedure


OpenWindow(0, 0, 0, ex, ey, "Mandelbrot", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
CanvasGadget(0, 0, 0, ex, ey, #PB_Canvas_ClipMouse | #PB_Canvas_Keyboard)
CreateImage(0, ex, ey, 32)
InitPalette()
zoom = 1:dx=1

Repeat  
    
    Event = WaitWindowEvent() 
    
    If Event = #PB_Event_Gadget
        X = GetGadgetAttribute(0, #PB_Canvas_MouseX)
        Y = GetGadgetAttribute(0, #PB_Canvas_MouseY)
        z = GetGadgetAttribute(0, #PB_Canvas_WheelDelta)
        b = GetGadgetAttribute(0, #PB_Canvas_Buttons)
        If b:dx=x-ax:dy=y-ay:Else:dx=0:dy=0:EndIf         
        ax=x:ay=y
    EndIf
    
    If z Or dx Or dy
        zoom=zoom*(1-0.1*z)
        fx-dx*zoom/ey
        fy+dy*zoom/ey
        For i=0 To ndt-1:Thread(i)=CreateThread(@fractale(),i):Next
        For i=0 To ndt-1:If Thread(i) : WaitThread(thread(i)):EndIf:Next
        ;fractale(3)
        StartDrawing(ImageOutput(0))
        CopyMemory(@bmp(0,0),DrawingBuffer(),ex*ey*4)
        StopDrawing()
        
        StartDrawing(CanvasOutput(0))
        DrawImage(ImageID(0),0,0,ex,ey)
        DrawingMode(#PB_2DDrawing_Transparent )
        DrawText(10,10,"Z:" + zoom+ "  X:"+x +"  Y:"+y)
        StopDrawing()
    EndIf
        
Until Event = #PB_Event_CloseWindow
manababel
Messages : 136
Inscription : jeu. 14/mai/2020 7:40

Re: conversion de routine PB en ASM - Mandelbrot

Message par manababel »

la première version calcule pixel par pixel
la deuxième version calcule 2 pixels en même temps
avec la première version 1 je tourne à 85 millisecondes
avec la deuxième version 2 je tourne à 50 millisecondes

pour la deuxième version, il y a eu un problème de copier coller, il manque en début de programme :

Global imax=255
Global lg=600
Global ht=400
Avatar de l’utilisateur
Micoute
Messages : 2522
Inscription : dim. 02/oct./2011 16:17
Localisation : 35520 La Mézière

Re: conversion de routine PB en ASM - Mandelbrot

Message par Micoute »

Merci pour le partage, c'est de toute beauté.
Microsoft Windows 10 Famille 64 bits : Carte mère : ASRock 970 Extreme3 R2.0 : Carte Graphique NVIDIA GeForce RTX 3080 : Processeur AMD FX 6300 6 cœurs 12 threads 3,50 GHz PB 5.73 PB 6.00 LTS (x64)
Un homme doit être poli, mais il doit aussi être libre !
manababel
Messages : 136
Inscription : jeu. 14/mai/2020 7:40

Re: conversion de routine PB en ASM - Mandelbrot

Message par manababel »

Voici une version qui permet de calculer 4 pixels en même temps.

Code : Tout sélectionner

Global lg=800
Global ht=600
Global infinity.f=4
Global limit=255
Global Dim tab.d(4)
Global Dim col.q(4)

Procedure.q madelbrot(a.d,b.d,scale.d)
  
  cp.q=limit
  fo.d=infinity
  add.q=1
  For i=0 To 3
    tab(i)=a
    a=a+scale
  Next
  
  t=@tab()
  c=@col()
  !mov rdx,[p.v_t]
  !VPBROADCASTQ ymm8,[p.v_add] ;1  1  1  1
  !VBROADCASTSD ymm7,[p.v_fo] ; 4  4  4  4
  !vmovupd ymm0,[rdx] ;         a1 a2 a3 a4
  !VBROADCASTSD ymm1,[p.v_b] ;  b  b  b  b
  !vmovupd ymm2,ymm0 ; c=a
  !vmovupd ymm3,ymm1 ; d=b
  !vpxor ymm9,ymm9,ymm9
  
  !xor rax,rax
  !boucle:
   
    !vmovupd ymm4,ymm0
    !vmovupd ymm5,ymm1
    !vmulpd ymm4,ymm4,ymm4 ; aa = a * a
    !vmulpd ymm5,ymm5,ymm5 ; bb = b * b
  
    !vmovupd ymm6,ymm4
    !vaddpd ymm6,ymm6,ymm5
    
    ;If (aa + bb) > 4.0:Goto fin:EndIf
    !vcmpltpd ymm6,ymm6,ymm7
    !vmovmskpd rcx,ymm6
    !vandpd  ymm6,ymm6,ymm8 ; and 1
    !vpaddq ymm9,ymm9,ymm6  ; add 1
    !Or rcx,rcx
    !jz fin
  
    ;x1 = 2 * x0 * x1 + x3
    !vmulpd ymm1,ymm1,ymm0 ; b*a
    !vaddpd ymm1,ymm1,ymm1 ; (a*b)*2
    !vaddpd ymm1,ymm1,ymm3 ; (a*b*2)+d
   
    ;x0 = x4 - x5 + c
    !vmovupd ymm0,ymm4 ; a = aa
    !vsubpd ymm0,ymm0,ymm5 ; a = aa - bb
    !vaddpd ymm0,ymm0,ymm2 ;  a = aa - bb + c
   
    !inc rax
    !cmp rax,[p.v_cp]
  !jb boucle
  
  !fin:
  
  !mov rdx,[p.v_c]
  !vMOVDQU [rdx],ymm9
  
  !VZEROALL

EndProcedure


If OpenWindow(0, 0, 0, lg+20, ht+40, "Exemple...", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
 
  img = CreateImage(1, lg,ht,32)
  StartDrawing(ImageOutput(1))
  buf = DrawingBuffer()
  StopDrawing()
  
  yside.d = -2.5
  xside.d = 2.5
  tside.d = 1.25
  lside.d = -2
  xscale.d=xside/lg
  yscale.d=yside/ht
  cx.d=0
  cy.d=0

    t=ElapsedMilliseconds()
    For y=0 To ht-1
      cy = y * yscale + tside
      For x=0 To lg-1 Step 4
        cx = x * xscale + lside 
        var=madelbrot(cx,cy,xscale)
        PokeL((buf+((y*lg)+x+0)<<2),col(0)*$10101)
        PokeL((buf+((y*lg)+x+1)<<2),col(1)*$10101)
        PokeL((buf+((y*lg)+x+2)<<2),col(2)*$10101)
        PokeL((buf+((y*lg)+x+3)<<2),col(3)*$10101)
      Next
    Next
    t=ElapsedMilliseconds()-t

    StartDrawing(WindowOutput(0))
    DrawImage(img,10,30)
    DrawText(5,5,Str(t)+"  ")
    StopDrawing()   
    
    
Repeat
Avatar de l’utilisateur
Ar-S
Messages : 9478
Inscription : dim. 09/oct./2005 16:51
Contact :

Re: conversion de routine PB en ASM - Mandelbrot

Message par Ar-S »

Jolie Fractale.
23ms pour ton dernier code.
~~~~Règles du forum ~~~~
⋅.˳˳.⋅ॱ˙˙ॱ⋅.˳Ar-S ˳.⋅ॱ˙˙ॱ⋅.˳˳.⋅
W11x64 PB 6.x
Section HORS SUJET : ICI
LDV MULTIMEDIA : Dépannage informatique & mes Logiciels PB
UPLOAD D'IMAGES : Uploader des images de vos logiciels
Avatar de l’utilisateur
Guillot
Messages : 532
Inscription : jeu. 25/juin/2015 16:18

Re: conversion de routine PB en ASM - Mandelbrot

Message par Guillot »

4 pixels à la fois, je sais pas comment c'est possible...

une version améliorée de mon code avec ta routine (maintenant on zoom sur le pointeur de la souris (molette))

enlever le débogueur

PS : je m'aperçois d'un probleme: avec le débogueur ça plante dans la fonction mandelbrot,
j'ai mis certaines valeurs en local pour l'utiliser en multithread
ça vient sans doute de là...

Code : Tout sélectionner

EnableExplicit

Global limit=255,ex=1280,ey=720,ascale.d,scale.d,     i,x,y,z,b,ax,ay,dx,dy
Global Event,px.d, py.d

Global Dim pal.l(limit)
Global Dim Bmp.l(ey - 1, ex - 1)

Global ndt=CountCPUs(#PB_System_ProcessCPUs )
Global Dim Thread(ndt+1)

Procedure ColorBlend(color1.l, color2.l, blend.f)
    Protected r.w,g.w,b.w,a.w
    r=  Red(color1) + (Red(color2)     - Red(color1)) * blend
    g=Green(color1) + (Green(color2) - Green(color1)) * blend
    b= Blue(color1) + (Blue(color2) -   Blue(color1)) * blend
    a=Alpha(color1) + (Alpha(color2) - Alpha(color1)) * blend
    ProcedureReturn  RGBA(r,g,b,a)
EndProcedure

Procedure InitPalette(n=8)
    Protected i,j,c1,c2
    c2=Random($ffffff)
    For j = 0 To limit/n
        c1=c2
        c2=Random($ffffff)
        For i=0 To n-1
            pal(j*n+i) = ColorBlend(c1,c2, i/ n)
        Next
    Next
    pal(limit) = 0
EndProcedure

Procedure mandelbrot(a.d,b.d,scale.d,Array col.q(1))
    Protected cp.q,fo.d,add.q,     i,t,c
    Protected infinity.f=4
    Dim tab.d(4)
    cp=limit
    fo=infinity
    add=1
    For i=0 To 3
        tab(i)=a
        a+scale
    Next

    t=@tab()
    c=@col()
    !mov rdx,[p.v_t]
    !VPBROADCASTQ ymm8,[p.v_add] ;1  1  1  1
    !VBROADCASTSD ymm7,[p.v_fo] ; 4  4  4  4
    !vmovupd ymm0,[rdx] ;         a1 a2 a3 a4
    !VBROADCASTSD ymm1,[p.v_b] ;  b  b  b  b
    !vmovupd ymm2,ymm0 ; c=a
    !vmovupd ymm3,ymm1 ; d=b
    !vpxor ymm9,ymm9,ymm9
    
    !xor rax,rax
    !boucle:
    
    !vmovupd ymm4,ymm0
    !vmovupd ymm5,ymm1
    !vmulpd ymm4,ymm4,ymm4 ; aa = a * a
    !vmulpd ymm5,ymm5,ymm5 ; bb = b * b
    
    !vmovupd ymm6,ymm4
    !vaddpd ymm6,ymm6,ymm5
    
    ;If (aa + bb) > 4.0:Goto fin:EndIf
    !vcmpltpd ymm6,ymm6,ymm7
    !vmovmskpd rcx,ymm6
    !vandpd  ymm6,ymm6,ymm8 ; and 1
    !vpaddq ymm9,ymm9,ymm6  ; add 1
    !Or rcx,rcx
    !jz fin
    
    ;x1 = 2 * x0 * x1 + x3
    !vmulpd ymm1,ymm1,ymm0 ; b*a
    !vaddpd ymm1,ymm1,ymm1 ; (a*b)*2
    !vaddpd ymm1,ymm1,ymm3 ; (a*b*2)+d
    
    ;x0 = x4 - x5 + c
    !vmovupd ymm0,ymm4 ; a = aa
    !vsubpd ymm0,ymm0,ymm5 ; a = aa - bb
    !vaddpd ymm0,ymm0,ymm2 ;  a = aa - bb + c
    
    !inc rax
    !cmp rax,[p.v_cp]
    !jb boucle
    
    !fin:
    
    !mov rdx,[p.v_c]
    !vMOVDQU [rdx],ymm9
    
    !VZEROALL
    
EndProcedure


Procedure fractale(num)
    Protected i,j,jdeb,jfin  ,x.d,y.d
    Dim col.q(4)
    jdeb=ey/ndt*num
    jfin=ey/ndt*(num+1)
    
    For j=jdeb To jfin-1
        y=py+j*scale
        x=px
        For i=0 To ex-1 Step 4
            mandelbrot(x, y,scale,col())
            bmp(j,i+0)=pal(col(0))
            bmp(j,i+1)=pal(col(1))
            bmp(j,i+2)=pal(col(2))
            bmp(j,i+3)=pal(col(3))
            x+scale*4
        Next
    Next 
EndProcedure

OpenWindow(0, 0, 0, ex, ey, "Mandelbrot - Use mouse + wheel", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
CanvasGadget(0, 0, 0, ex, ey, #PB_Canvas_ClipMouse | #PB_Canvas_Keyboard)
CreateImage(0, ex, ey, 32)
InitPalette(8)
scale=2/ey
px=-ex/2*scale
py=-ey/2*scale
dx=1

Repeat  
    
    Event = WaitWindowEvent() 
    
    If Event = #PB_Event_Gadget
        X = GetGadgetAttribute(0, #PB_Canvas_MouseX)
        Y = GetGadgetAttribute(0, #PB_Canvas_MouseY)
        z = GetGadgetAttribute(0, #PB_Canvas_WheelDelta)
        b = GetGadgetAttribute(0, #PB_Canvas_Buttons)
        If b:dx=x-ax:dy=y-ay:Else:dx=0:dy=0:EndIf         
        ax=x:ay=y
    EndIf
    
    If z Or dx Or dy
        ascale=scale
        scale=scale*(1-0.1*z)
        px-dx/2*scale
        py+dy/2*scale
        px+(-x)*(scale-ascale)
        py+(y-ey)*(scale-ascale)
        For i=0 To ndt-1:Thread(i)=CreateThread(@fractale(),i):Next
        For i=0 To ndt-1:If Thread(i) : WaitThread(thread(i)):EndIf:Next
        ;ndt=1:fractale(0)
        StartDrawing(ImageOutput(0))
        CopyMemory(@bmp(0,0),DrawingBuffer(),ex*ey*4)
        StopDrawing()
        
        StartDrawing(CanvasOutput(0))
        DrawImage(ImageID(0),0,0,ex,ey)
        DrawingMode(#PB_2DDrawing_Transparent )
        DrawText(10,10,"Scale:" + scale+ "  X:"+px +"  Y:"+py)
        StopDrawing()
    EndIf
    
Until Event = #PB_Event_CloseWindow
manababel
Messages : 136
Inscription : jeu. 14/mai/2020 7:40

Re: conversion de routine PB en ASM - Mandelbrot

Message par manababel »

je regard le problème du plantage avec le debugger

avec la version qui calcule 2 pixels en même temps , le multithreading devrait fonctionner

par contre , la version qui calcule 4 pixels en meme temps, contient un tableau .
le tableau rend cette routine incompatible avec les threads.
Avatar de l’utilisateur
Guillot
Messages : 532
Inscription : jeu. 25/juin/2015 16:18

Re: conversion de routine PB en ASM - Mandelbrot

Message par Guillot »

si on active la gestion des threads (options de compilation) ça fonctionne aussi en mode débogueur
(curieux que ça fonctionne sans...)
manababel
Messages : 136
Inscription : jeu. 14/mai/2020 7:40

Re: conversion de routine PB en ASM - Mandelbrot

Message par manababel »

voici une version sans tableaux

Code : Tout sélectionner

Global limit=255,ex=1280,ey=720,ascale.d,scale.d,     i,x.q,y.q,z.q,b,ax,ay,dx,dy
Global Event,px.d, py.d

Global Dim pal.l(limit)
Global Dim Bmp.l(ey - 1, ex - 1)

Global ndt=CountCPUs(#PB_System_ProcessCPUs )
Global Dim Thread(ndt+1)

Procedure ColorBlend(color1.l, color2.l, blend.f)
    Protected r.w,g.w,b.w,a.w
    r=  Red(color1) + (Red(color2)     - Red(color1)) * blend
    g=Green(color1) + (Green(color2) - Green(color1)) * blend
    b= Blue(color1) + (Blue(color2) -   Blue(color1)) * blend
    a=Alpha(color1) + (Alpha(color2) - Alpha(color1)) * blend
    ProcedureReturn  RGBA(r,g,b,a)
EndProcedure

Procedure InitPalette(n=8)
    Protected i,j,c1,c2
    c2=Random($ffffff)
    For j = 0 To limit/n
        c1=c2
        c2=Random($ffffff)
        For i=0 To n-1
            pal(j*n+i) = ColorBlend(c1,c2, i/ n)
        Next
    Next
    pal(limit) = 0
EndProcedure

Procedure.q mandelbrot(a.d,b.d,scale.d)
    Protected cp.q,fo.d,add.q
    Protected infinity.f=8
    cp=limit
    fo=infinity
    add=1
    
    !movsd xmm1,[p.v_scale] ; 00 sc
    !movsd xmm2,[p.v_a] ;     00 a0
    !addpd xmm2,xmm1 ; a1
    !addpd xmm2,xmm1 ; a2
    !addpd xmm2,xmm1 ; a3
    !movddup xmm0,xmm2 ; a3 a3
    !subpd xmm0,xmm1   ; a3 a2
    !VINSERTF128 ymm0,ymm0,xmm0,1 ; a3 a2 a3 a2
    !subpd xmm0,xmm1   ; a1
    !movddup xmm0,xmm0 ; a3 a2 a1 a1
    !subpd xmm0,xmm1   ; a3 a2 a1 a0
    
    !VPBROADCASTQ ymm8,[p.v_add] ;1  1  1  1
    !VBROADCASTSD ymm7,[p.v_fo] ; 4  4  4  4
    !VBROADCASTSD ymm1,[p.v_b] ;  b  b  b  b
    !vmovupd ymm2,ymm0 ; c=a
    !vmovupd ymm3,ymm1 ; d=b
    !vpxor ymm9,ymm9,ymm9
   
    !xor rax,rax
    !boucle:
   
    !vmovupd ymm4,ymm0
    !vmovupd ymm5,ymm1
    !vmulpd ymm4,ymm4,ymm4 ; aa = a * a
    !vmulpd ymm5,ymm5,ymm5 ; bb = b * b
   
    !vmovupd ymm6,ymm4
    !vaddpd ymm6,ymm6,ymm5
   
    ;If (aa + bb) > 4.0:Goto fin:EndIf
    !vcmpltpd ymm6,ymm6,ymm7
    !vmovmskpd rcx,ymm6
    !vandpd  ymm6,ymm6,ymm8 ; and 1
    !vpaddq ymm9,ymm9,ymm6  ; add 1
    !Or rcx,rcx
    !jz fin
   
    ;x1 = 2 * x0 * x1 + x3
    !vmulpd ymm1,ymm1,ymm0 ; b*a
    !vaddpd ymm1,ymm1,ymm1 ; (a*b)*2
    !vaddpd ymm1,ymm1,ymm3 ; (a*b*2)+d
   
    ;x0 = x4 - x5 + c
    !vmovupd ymm0,ymm4 ; a = aa
    !vsubpd ymm0,ymm0,ymm5 ; a = aa - bb
    !vaddpd ymm0,ymm0,ymm2 ;  a = aa - bb + c
   
    !inc rax
    !cmp rax,[p.v_cp]
    !jb boucle
   
    !fin:
    ;00c1 00c2  00c3  00c4
    !pxor xmm2,xmm2
    !Vpshufd ymm9,ymm9,216 ; 0000 c1c2  0000  c3c4
    !VEXTRACTI128 xmm0,ymm9,1 ; 0000 c1c2
    !packusdw xmm9,xmm2 ; 32 -> 16 bits
    !packusdw xmm0,xmm2
    !psllq xmm0,32
    !por xmm9,xmm0 ; c1c2c3c4  4*16
    !movq rax,xmm9
   
    !VZEROALL
   ProcedureReturn
EndProcedure




Procedure fractale(num)
    Protected i,j,jdeb,jfin  ,x.d,y.d
    jdeb=ey/ndt*num
    jfin=ey/ndt*(num+1)
   
    For j=jdeb To jfin-1
        y=py+j*scale
        x=px
        For i=0 To ex-1 Step 4
            var=mandelbrot(x, y,scale)
            bmp(j,i+0)=pal(var & 255)
            bmp(j,i+1)=pal(var >>16 & 255)
            bmp(j,i+2)=pal(var >>32 & 255)
            bmp(j,i+3)=pal(var >>48 & 255)
            x+scale*4
        Next
    Next
EndProcedure

OpenWindow(0, 0, 0, ex, ey, "Mandelbrot - Use mouse + wheel", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
CanvasGadget(0, 0, 0, ex, ey, #PB_Canvas_ClipMouse | #PB_Canvas_Keyboard)
CreateImage(0, ex, ey, 32)
InitPalette(8)
scale=2/ey
px=-ex/2*scale
py=-ey/2*scale
dx=1

Repeat 
   
    Event = WaitWindowEvent()
   
    If Event = #PB_Event_Gadget
        X = GetGadgetAttribute(0, #PB_Canvas_MouseX)
        Y = GetGadgetAttribute(0, #PB_Canvas_MouseY)
        z = GetGadgetAttribute(0, #PB_Canvas_WheelDelta)
        b = GetGadgetAttribute(0, #PB_Canvas_Buttons)
        If b:dx=x-ax:dy=y-ay:Else:dx=0:dy=0:EndIf         
        ax=x:ay=y
    EndIf
   
    ;If z Or dx Or dy
        ascale=scale
        scale=scale*(1-0.1*z)
        px-dx/2*scale
        py+dy/2*scale
        px+(-x)*(scale-ascale)
        py+(y-ey)*(scale-ascale)
        t=ElapsedMilliseconds()
        For i=0 To ndt-1:Thread(i)=CreateThread(@fractale(),i):Next
        For i=0 To ndt-1:If Thread(i) : WaitThread(thread(i)):EndIf:Next
        t=ElapsedMilliseconds()-t
        ;ndt=1:fractale(0)
        StartDrawing(ImageOutput(0))
        CopyMemory(@bmp(0,0),DrawingBuffer(),ex*ey*4)
        StopDrawing()
       
        StartDrawing(CanvasOutput(0))
        DrawImage(ImageID(0),0,0,ex,ey)
        DrawingMode(#PB_2DDrawing_Transparent )
        DrawText(10,10,"Scale:" + scale+ "  X:"+px +"  Y:"+py)
        DrawText(10,25,Str(t))
        StopDrawing()
    ;EndIf
   
Until Event = #PB_Event_CloseWindow
Avatar de l’utilisateur
Guillot
Messages : 532
Inscription : jeu. 25/juin/2015 16:18

Re: conversion de routine PB en ASM - Mandelbrot

Message par Guillot »

merci
mais comme je te disais plus haut, si on active la gestion des threads (options de compilation) ça fonctionne aussi en mode débogueur
manababel
Messages : 136
Inscription : jeu. 14/mai/2020 7:40

Re: conversion de routine PB en ASM - Mandelbrot

Message par manababel »

rajout de l'option "Julia" dans la procédure "Mandelbrot"


Code : Tout sélectionner

Global limit=255,ex=1280,ey=720,ascale.d,scale.d,     i,x.q,y.q,z.q,b,ax,ay,dx,dy
Global Event,px.d, py.d
Global option=0
Global Dim pal.l(limit+1)
Global Dim Bmp.l(ey - 1, ex - 1)

Global ndt=CountCPUs(#PB_System_ProcessCPUs )
Global Dim Thread(ndt+1)

Procedure ColorBlend(color1.l, color2.l, blend.f)
    Protected r.w,g.w,b.w,a.w
    r=  Red(color1) + (Red(color2)     - Red(color1)) * blend
    g=Green(color1) + (Green(color2) - Green(color1)) * blend
    b= Blue(color1) + (Blue(color2) -   Blue(color1)) * blend
    a=Alpha(color1) + (Alpha(color2) - Alpha(color1)) * blend
    ProcedureReturn  RGBA(r,g,b,a)
EndProcedure

Procedure InitPalette1(n=8)
    Protected i,j,c1,c2
    c2=Random($ffffff)
    For j = 0 To limit/n
        c1=c2
        c2=Random($ffffff)
        For i=0 To n-1
            pal(j*n+i) = ColorBlend(c1,c2, i/ n)
        Next
    Next
    pal(limit) = 0
EndProcedure

Procedure InitPalette2()
  For i=0 To limit-1
     r=Mod(i,32)*7;
     g=Mod(i,16)*14;
     b=Mod(i,128)*2;
    r=r<<16;
    g=g<<8;
    pal(i)=r+g+b
  Next
EndProcedure

Procedure InitPalette3()
  For i=0 To limit-1
    If (i<63):r=i: EndIf
    If (i>62 And i<128):g=i: EndIf
    If (i>127 And i<192):b=i: EndIf
    If (i>191) :r=i/4 :g=r :b=r : EndIf
    r=r<<16;
    g=g<<8;
    pal(i)=r+g+b
  Next
EndProcedure

Procedure InitPalette4(v=0)
  For i=0 To limit-1
    If (v=0):r=255*Cos(Radian(i*2)):g=255*Cos(Radian(i*3)):b=255*Cos(Radian(i*4)):EndIf
    If (v=1):r=255*Cos(i*3):g=255*Cos(i*4):b=255*Cos(i*2):EndIf
    If (v=2):r=255*Cos(i*4):g=255*Cos(i*2):b=255*Cos(i*3):EndIf
    If (v=3):r=255*Cos(i*2):g=255*Cos(i*4):b=255*Cos(i*3):EndIf
    r=r<<16;
    g=g<<8;
    pal(i)=r+g+b
  Next
EndProcedure

Procedure InitPalette5(v=0)
  For i=0 To limit-1
    c=i*8*$10000+i*4*$100+i
    pal(i)=c
  Next
EndProcedure

Procedure.q mandelbrot(a.d,b.d,scale.d,opt=0)

    cp.q=limit
    fo.d=4
    add.q=1
    
    juliacx.d=-0.7
    juliacy.d=0.27015

    !movsd xmm1,[p.v_scale] ; 00 sc
    !movsd xmm2,[p.v_a] ;     00 a0
    !addpd xmm2,xmm1 ; a1
    !addpd xmm2,xmm1 ; a2
    !addpd xmm2,xmm1 ; a3
    !movddup xmm0,xmm2 ; a3 a3
    !subpd xmm0,xmm1   ; a3 a2
    !VINSERTF128 ymm0,ymm0,xmm0,1 ; a3 a2 a3 a2
    !subpd xmm0,xmm1   ; a1
    !movddup xmm0,xmm0 ; a3 a2 a1 a1
    !subpd xmm0,xmm1   ; a3 a2 a1 a0
    
    !mov rax,1
    !movq xmm8,rax
    !VPBROADCASTQ ymm8,xmm8
    !VBROADCASTSD ymm7,[p.v_fo] ; 4  4  4  4
    !VBROADCASTSD ymm1,[p.v_b] ;  b  b  b  b
    !vmovupd ymm2,ymm0 ; c=a
    !vmovupd ymm3,ymm1 ; d=b
    !vpxor ymm9,ymm9,ymm9
    
    !xor rax,rax
    !cmp rax,[p.v_opt]
    !je NotJulia
      !VBROADCASTSD ymm2,[p.v_juliacx]
      !VBROADCASTSD ymm3,[p.v_juliacy]
    !NotJulia:
    !xor rax,rax
    !boucle:
   
    !vmovupd ymm4,ymm0
    !vmovupd ymm5,ymm1
    !vmulpd ymm4,ymm4,ymm4 ; aa = a * a
    !vmulpd ymm5,ymm5,ymm5 ; bb = b * b
   
    !vmovupd ymm6,ymm4
    !vaddpd ymm6,ymm6,ymm5
   
    ;If (aa + bb) > 4.0:Goto fin:EndIf
    !vcmpltpd ymm6,ymm6,ymm7
    !vmovmskpd rcx,ymm6
    !vandpd  ymm6,ymm6,ymm8 ; and 1
    !vpaddq ymm9,ymm9,ymm6  ; add 1
    !Or rcx,rcx
    !jz fin
   
    ;x1 = 2 * x0 * x1 + x3
    !vmulpd ymm1,ymm1,ymm0 ; zx*zy
    !vaddpd ymm1,ymm1,ymm1 ; (zx*zy)*2
    !vaddpd ymm1,ymm1,ymm3 ; (zx*zy*2)+cy
   
    ;x0 = x4 - x5 + c
    !vmovupd ymm0,ymm4 ; a = aa
    !vsubpd ymm0,ymm0,ymm5 ; a = aa - bb
    !vaddpd ymm0,ymm0,ymm2 ;  a = zx - zy + cx
   
    !inc rax
    !cmp rax,[p.v_cp]
    !jb boucle
   
    !fin:
    ;00c1 00c2  00c3  00c4
    !pxor xmm2,xmm2
    !Vpshufd ymm9,ymm9,216 ; 0000 c1c2  0000  c3c4
    !VEXTRACTI128 xmm0,ymm9,1 ; 0000 c1c2
    !packusdw xmm9,xmm2 ; 32 -> 16 bits
    !packusdw xmm0,xmm2
    !psllq xmm0,32
    !por xmm9,xmm0 ; c1c2c3c4  4*16
    !movq rax,xmm9
   
    !VZEROALL
   ProcedureReturn
 EndProcedure




Procedure fractale(num)
    Protected i,j,jdeb,jfin  ,x.d,y.d
    jdeb=ey/ndt*num
    jfin=ey/ndt*(num+1)
   
    For j=jdeb To jfin-1
        y=py+j*scale
        x=px
        For i=0 To ex-1 Step 4
            var=mandelbrot(x, y,scale,option)
            bmp(j,i+0)=pal(var & 255)
            bmp(j,i+1)=pal(var >>16 & 255)
            bmp(j,i+2)=pal(var >>32 & 255)
            bmp(j,i+3)=pal(var >>48 & 255)
            x+scale*4
        Next
    Next
EndProcedure

OpenWindow(0, 0, 0, ex, ey, "Mandelbrot - Use mouse + wheel", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
CanvasGadget(0, 0, 0, ex, ey, #PB_Canvas_ClipMouse | #PB_Canvas_Keyboard)
CreateImage(0, ex, ey, 32)
InitPalette1(8)
scale=2/ey
px=-ex/2*scale
py=-ey/2*scale
dx=1

If CreateMenu(0, WindowID(0))
  MenuTitle( "choix")
    MenuItem( 1, "Mandelbrot")
    MenuItem( 2, "Julia")
  
    MenuTitle( "color")
    MenuItem( 10, "v1")
    MenuItem( 11, "v2")
    MenuItem( 12, "v3")
    MenuItem( 13, "v4")
    MenuItem( 14, "v5")
  EndIf
Repeat 
   
    Event = WaitWindowEvent()
    
    Select event
      Case #PB_Event_Menu
        Select EventMenu()
          Case 1 : rnd=1
            option=0
          Case 2 : rnd=1
            option=1
          Case 10
            InitPalette1() : rnd=1
          Case 11
            InitPalette2() : rnd=1
          Case 12
            InitPalette3() : rnd=1
          Case 13
            InitPalette4() : rnd=1  
          Case 14
            InitPalette5() : rnd=1 
            
        EndSelect
    EndSelect
    
        
    If Event = #PB_Event_Gadget
        X = GetGadgetAttribute(0, #PB_Canvas_MouseX)
        Y = GetGadgetAttribute(0, #PB_Canvas_MouseY)
        z = GetGadgetAttribute(0, #PB_Canvas_WheelDelta)
        b = GetGadgetAttribute(0, #PB_Canvas_Buttons)
        If b:dx=x-ax:dy=y-ay:Else:dx=0:dy=0:EndIf         
        ax=x:ay=y
    EndIf
   
    If z Or dx Or dy Or rnd=1
       rnd=0
        ascale=scale
        scale=scale*(1-0.1*z)
        px-dx/2*scale
        py+dy/2*scale
        px+(-x)*(scale-ascale)
        py+(y-ey)*(scale-ascale)
        t=ElapsedMilliseconds()
        For i=0 To ndt-1:Thread(i)=CreateThread(@fractale(),i):Next
        For i=0 To ndt-1:If Thread(i) : WaitThread(thread(i)):EndIf:Next
        t=ElapsedMilliseconds()-t
        ;ndt=1:fractale(0)
        StartDrawing(ImageOutput(0))
        CopyMemory(@bmp(0,0),DrawingBuffer(),ex*ey*4)
        StopDrawing()
       
        StartDrawing(CanvasOutput(0))
        DrawImage(ImageID(0),0,0,ex,ey)
        DrawingMode(#PB_2DDrawing_Transparent )
        DrawText(10,10,"Scale:" + scale+ "  X:"+px +"  Y:"+py)
        DrawText(10,25,Str(t))
        StopDrawing()
    EndIf
   
Until Event = #PB_Event_CloseWindow
manababel
Messages : 136
Inscription : jeu. 14/mai/2020 7:40

Re: conversion de routine PB en ASM - Mandelbrot

Message par manababel »

j'ai modifié la première version du programme que vous avez posté.

1-le tableau 'BMP' a ete agrandit pour que chaque thread puisse calculer deux lignes de plus.
-le problème, c'est que les threads ne sont pas exécutés dans l'ordre
-donc par moments on voit toujours quelques problèmes de ligne

2-puis ajoutés la procédure "fractale test" à la fin du programme
-ici on attend que tous les threads soient finis, puis on calcule toutes les lignes qui posent problème avec les threads
-en améliorant cette partie, on doit pouvoir supprimer le calcul des deux lignes supplémentaire


Code : Tout sélectionner

Global imax=255,lg=1280,ht=800

Global Dim pal.l(4096)
Global Dim Bmp.l(ht + 2, lg - 1)
Global Dim t.w(ht + 1 , lg + 1)
Global Dim cx.d(lg - 1)
Global Dim cy.d(ht - 1)

Global ndt=CountCPUs(#PB_System_ProcessCPUs )
Global Dim Thread(ndt+1)
Global Event.i,z.d, fx.d, fy.d,zoom.d

Global Dim test(ndt+1)

Procedure ColorBlend(color1.l, color2.l, blend.f)
    Protected r.w,g.w,b.w,a.w
    r=  Red(color1) + (Red(color2)     - Red(color1)) * blend
    g=Green(color1) + (Green(color2) - Green(color1)) * blend
    b= Blue(color1) + (Blue(color2) -   Blue(color1)) * blend
    a=Alpha(color1) + (Alpha(color2) - Alpha(color1)) * blend
    ProcedureReturn  RGBA(r,g,b,a)
EndProcedure

Procedure InitPalette()
    r.b:g.b:b.b
    Dim co.l(65)
    ;co = Array(&HFF0000, &HFF00FF, &HFF&, &HFFFF&, &H4400&, &HFF0000, &HFF00FF, &HFF&, &HFFFF&, &H4400&, &HFF0000, &HFF00FF, &HFF&, &HFFFF&, &H4400&, &HFF0000, &HFF00FF, &HFF&, &HFFFF&, &H4400&)
    For i = 0 To 64: co(i) = Random( $FFFFFF): Next
    cl = 8
    For j = 0 To 63
        For i = 0 To cl - 1
          ci = j * cl + i
            pal(ci) = ColorBlend(co(j), co(j + 1), i/ cl)
        Next
    Next
    pal(imax) = 0
    ;For i = 0 To 255: CoRGB pal(i), b, g, r: Line (i * 2, 0)-(i * 2 + 1, 10), RGB(b, g, r), BF: next
EndProcedure

Procedure mandelbrotx1(a.d,b.d)

  cp.q=imax
  fo.d=4.0
  !movsd xmm7,[p.v_fo]
  !movsd xmm0,[p.v_a] ; a
  !movsd xmm1,[p.v_b] ; b
  !movupd xmm2,xmm0 ; c=a
  !movupd xmm3,xmm1 ; d=b
  !xor rax,rax
  !boucle1:
   
  !movupd xmm4,xmm0
  !movupd xmm5,xmm1
  !mulpd xmm4,xmm4 ; aa = a * a
  !mulpd xmm5,xmm5 ; bb = b * b

  ;If (aa + bb) > 4.0:Goto fin:EndIf
  !movupd xmm6,xmm4
  !addpd xmm6,xmm5
  !VCMPGTPD xmm6,xmm6,xmm7
  !movq rdx,xmm6
  !cmp rdx,0
  !jne fin1

  ;x1 = 2 * x0 * x1 + x3
  !mulpd xmm1,xmm0 ; b*a
  !addpd xmm1,xmm1 ; (a*b)*2
  !addpd xmm1,xmm3 ; (a*b*2)+d

  ;x0 = x4 - x5 + c
  !movupd xmm0,xmm4 ; a = aa
  !subpd xmm0,xmm5 ; a = aa - bb
  !addpd xmm0,xmm2 ;  a = aa - bb + c

  !inc rax
  !cmp rax,[p.v_cp]
  !jb boucle1

  !fin1:

  ProcedureReturn ; rax

EndProcedure

Procedure.q mandelbrotx4(a.d,b.d,scale.d,opt=0)

    cp.q=imax
    fo.d=4
    
    juliacx.d=-0.7
    juliacy.d=0.27015

    !movsd xmm1,[p.v_scale] ; 00 sc
    !movsd xmm2,[p.v_a] ;     00 a0
    !addpd xmm2,xmm1 ; a1
    !addpd xmm2,xmm1 ; a2
    !addpd xmm2,xmm1 ; a3
    !movddup xmm0,xmm2 ; a3 a3
    !subpd xmm0,xmm1   ; a3 a2
    !VINSERTF128 ymm0,ymm0,xmm0,1 ; a3 a2 a3 a2
    !subpd xmm0,xmm1   ; a1
    !movddup xmm0,xmm0 ; a3 a2 a1 a1
    !subpd xmm0,xmm1   ; a3 a2 a1 a0
    
    !mov rax,1
    !movq xmm8,rax
    !VPBROADCASTQ ymm8,xmm8
    !VBROADCASTSD ymm7,[p.v_fo] ; 4  4  4  4
    !VBROADCASTSD ymm1,[p.v_b] ;  b  b  b  b
    !vmovupd ymm2,ymm0 ; c=a
    !vmovupd ymm3,ymm1 ; d=b
    !vpxor ymm9,ymm9,ymm9
    
    !xor rax,rax
    !cmp rax,[p.v_opt]
    !je NotJulia
      !VBROADCASTSD ymm2,[p.v_juliacx]
      !VBROADCASTSD ymm3,[p.v_juliacy]
    !NotJulia:
    !xor rax,rax
    !boucle:
   
    !vmovupd ymm4,ymm0
    !vmovupd ymm5,ymm1
    !vmulpd ymm4,ymm4,ymm4 ; aa = a * a
    !vmulpd ymm5,ymm5,ymm5 ; bb = b * b
   
    !vmovupd ymm6,ymm4
    !vaddpd ymm6,ymm6,ymm5
   
    ;If (aa + bb) > 4.0:Goto fin:EndIf
    !vcmpltpd ymm6,ymm6,ymm7
    !vmovmskpd rcx,ymm6
    !vandpd  ymm6,ymm6,ymm8 ; and 1
    !vpaddq ymm9,ymm9,ymm6  ; add 1
    !Or rcx,rcx
    !jz fin
   
    ;x1 = 2 * x0 * x1 + x3
    !vmulpd ymm1,ymm1,ymm0 ; zx*zy
    !vaddpd ymm1,ymm1,ymm1 ; (zx*zy)*2
    !vaddpd ymm1,ymm1,ymm3 ; (zx*zy*2)+cy
   
    ;x0 = x4 - x5 + c
    !vmovupd ymm0,ymm4 ; a = aa
    !vsubpd ymm0,ymm0,ymm5 ; a = aa - bb
    !vaddpd ymm0,ymm0,ymm2 ;  a = zx - zy + cx
   
    !inc rax
    !cmp rax,[p.v_cp]
    !jb boucle
   
    !fin:
    ;00c1 00c2  00c3  00c4
    !pxor xmm2,xmm2
    !Vpshufd ymm9,ymm9,216 ; 0000 c1c2  0000  c3c4
    !VEXTRACTI128 xmm0,ymm9,1 ; 0000 c1c2
    !packusdw xmm9,xmm2 ; 32 -> 16 bits
    !packusdw xmm0,xmm2
    !psllq xmm0,32
    !por xmm9,xmm0 ; c1c2c3c4  4*16
    !movq rax,xmm9
   
    !VZEROALL
   ProcedureReturn
 EndProcedure

Procedure fractale(num)
    Protected i,j,jdeb,jfin  ,px,py , scale.d , var.q
    jdeb=ht/ndt*num
    jfin=ht/ndt*(num+1)+2
    test(num)=jfin
    
    If jfin>=ht : ProcedureReturn : EndIf
      
    For i = 0 To lg - 1   : cx(i) = zoom * ((2 * i) / lg -1) * lg / ht + fx: Next
    For j = jdeb To jfin-1: cy(j) = zoom * ((2 * j) / ht -1) + fy: Next
   scale = cx(2) - cx(0)
    dd=2
    For jj = jdeb/dd To jfin/dd:j=jj*dd
        For ii = 0 To lg/dd - 1 Step 4:i=ii*dd
          var = mandelbrotx4((cx(i)), (cy(j)),scale)
          t(j,i+0) = ( var >> 00 ) & 255
          t(j,i+2) = ( var >> 16 ) & 255
          t(j,i+4) = ( var >> 32 ) & 255
          t(j,i+6) = ( var >> 48 ) & 255
        Next
    Next
    While dd > 1
        d = dd / 2: c1 = 0: c2 = 0: c3 = 0
       
        For jj = jdeb/dd+2  To jfin/dd :j=jj*dd-d
            For ii = 2 To lg/dd - 1 :i=ii*dd-d
                pa = t(j - d,i - d)
                pb = t(j + d,i - d)
                pc = t( j + d,i + d)
                pd = t(j - d,i + d)
                If pa = pb And pa = pc And pa = pd
                    pe = pa: c1 = c1 + 1
                Else
                    pe = mandelbrotx1((cx(i)), (cy(j)))
                EndIf
                pf = t( j,i - dd)
                pg = t(j - dd,i)
                t(j,i) = pe
                If pe = pc And pe = pf And pe = pd
                    t( j,i + d) = pe: c2 = c2 + 1
                Else
                    t(j,i + d) = mandelbrotx1(cx(i + d), cy(j))
                EndIf
                If pe = pc And pe = pg And pe = pb
                    t(j + d,i) = pe: c3 = c3 + 1
                Else
                    t(j + d,i) = mandelbrotx1((cx(i)), (cy(j + d)))
                EndIf
            Next
        Next
        dd = dd / 2
        ;Debug lg * ht, c1; c2; c3
    Wend
   
    For j=jdeb To jfin-1
        For i=0 To lg-1
            bmp(j+1,i)=pal(t(j,i))
        Next
    Next

   
EndProcedure

Procedure fractale_test()
  
    For k=0 To ndt-1
      jdeb=test(k)
      jfin=jdeb
   
    If jfin>=ht : ProcedureReturn : EndIf
    ;For i = 0 To lg - 1   : cx(i) = zoom * ((2 * i) / lg -1) * lg / ht + fx: Next
    ;For j = jdeb To jfin-1: cy(j) = zoom * ((2 * j) / ht -1) + fy: Next
   
    
    dd=2
    ;While dd > 1
    d = dd / 2: c1 = 0: c2 = 0: c3 = 0
    j=jfin

        ;For jj = jdeb/dd+2  To jfin/dd :j=jj*dd-d
            For ii = 2 To lg/dd - 1 :i=ii*dd-d
                pa = t(j - 1,i - 1)
                pb = t(j + 1,i - 1)
                pc = t( j + 1,i + 1)
                pd = t(j - 1,i + 1)
                If pa = pb And pa = pc And pa = pd
                    pe = pa: c1 = c1 + 1
                Else
                    pe = mandelbrotx1((cx(i)), (cy(j)))
                EndIf
                pf = t( j,i - dd)
                pg = t(j - dd,i)
                t(j,i) = pe
                If pe = pc And pe = pf And pe = pd
                    t( j,i + 1) = pe: c2 = c2 + 1
                Else
                    t(j,i + 1) = mandelbrotx1(cx(i + 1), cy(j))
                EndIf
                If pe = pc And pe = pg And pe = pb
                    t(j + 1,i) = pe: c3 = c3 + 1
                Else
                    t(j + 1,i) = mandelbrotx1((cx(i)), (cy(j + 1)))
                EndIf
            Next
        ;Next
        dd = dd / 2
        ;Debug lg * ht, c1; c2; c3
    ;Wend
   
    ;For j=jdeb To jfin-1
        For i=0 To lg-1
            bmp(j+1,i)=pal(t(j,i))
        Next
      ;Next
      
    Next
    
    EndProcedure
    
    
    
OpenWindow(0, 0, 0, lg, ht, "Mandelbrot", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
CanvasGadget(0, 0, 0, lg, ht, #PB_Canvas_ClipMouse | #PB_Canvas_Keyboard)
CreateImage(0, lg, ht, 32)
InitPalette()
zoom = 1:dx=1

Repeat 
   
    Event = WindowEvent()
   
    If Event = #PB_Event_Gadget
        X = GetGadgetAttribute(0, #PB_Canvas_MouseX)
        Y = GetGadgetAttribute(0, #PB_Canvas_MouseY)
        z = GetGadgetAttribute(0, #PB_Canvas_WheelDelta)
        b = GetGadgetAttribute(0, #PB_Canvas_Buttons)
        If b:dx=x-ax:dy=y-ay:Else:dx=0:dy=0:EndIf         
        ax=x:ay=y
    EndIf
   
    If z Or dx Or dy
        zoom=zoom*(1-0.1*z)
        fx-dx*zoom/ht
        fy+dy*zoom/ht
        ;ti=ElapsedMilliseconds()
        For i=0 To ndt-1:Thread(i)=CreateThread(@fractale(),i):Next
        For i=0 To ndt-1:If Thread(i) : WaitThread(thread(i)):EndIf:Next
        fractale_test()
        StartDrawing(ImageOutput(0))
        CopyMemory(@bmp(0,0),DrawingBuffer(),lg*ht*4)
        StopDrawing()
        ;ti=ElapsedMilliseconds()-ti
        
        StartDrawing(CanvasOutput(0))
        DrawImage(ImageID(0),0,0,lg,ht)
        DrawingMode(#PB_2DDrawing_Transparent )
        DrawText(10,10,"Z:" + zoom+ "  X:"+x +"  Y:"+y)
        ;DrawText(10,25,Str(ti))
        StopDrawing()
    EndIf
       
Until Event = #PB_Event_CloseWindow
Répondre