Author Topic: Mandelbrot  (Read 5962 times)

0 Members and 2 Guests are viewing this topic.

Peter

  • Guest
Mandelbrot
« on: October 14, 2015, 10:26:07 AM »
Hello,
unfortunately, Mike did not show his Mandelbrot collection!
I am interested in the speed!
Code: [Select]
include "asm.inc"
window 512,512,1

float bx,by, int bw=512,bh=512,col
float sx=-2.2,sy=-1.7,sw=3.4,sh=3.4
float t,gx,gy,zx,zy,nzx,r,v,b

t = ticks
for x=0 to bw
for y=0 to bh
    gx=x/bw*sw+sx
    gy=y/bh*sh+sy
    zx=gx
    zy=gy
for c=0 to 255
    col = c
    nzx = zx*zx - zy*zy + gx
    zy  = 2*zx*zy+gy
    zx  = nzx
    if zx*zx + zy*zy > 4 then
       col = c
       exit for
    end if
next
r = col     
v = col*32 
b = col*64 
color r,v,b
setpixel x,y
next
next

color 255,255,255
text 0,0,14, str((ticks()-t)/1000,4) + " secondes"

waitkey
winExit



[attachment deleted by admin]

JRS

  • Guest
Re: Mandelbrot
« Reply #1 on: October 14, 2015, 01:00:18 PM »
@Peter: Don't make me dig out the word count challenge again;D

Mike Lobanovsky

  • Guest
Re: Mandelbrot
« Reply #2 on: October 14, 2015, 01:50:55 PM »
Hi Peter,

My absolute best for a 400x400 px window was



Eros Olmi's thinBasic + O2 add-on for 400x400 pixels was a very, very close second-best with



which, for all intents and purposes and given GetTickCount() inaccuracy, would be practically the same. :)

The entire epopee can still be found starting with this message by Eros on the FBSL forum.



Can you please compile your O2 code to an exe for me, Peter?

Peter

  • Guest
Re: Mandelbrot
« Reply #3 on: October 14, 2015, 03:09:47 PM »
Quote
Can you please compile your O2 code to an exe for me, Peter?
but sure!


[attachment deleted by admin]

Mike Lobanovsky

  • Guest
Re: Mandelbrot
« Reply #4 on: October 15, 2015, 02:08:55 AM »
Thank you!

Not too bad, Peter, I'd say, very decent for a typical BASIC penis-metering contest using some form of PSet. :)

My 15/31 msec duel with Eros used direct writes into the DIB pixel array rather than GDI function calls. That is why it yielded much faster results, but unfortunately this method inverts the original MS Windows colors and therefore is rarely used in practice.

[attachment deleted by admin]
« Last Edit: October 15, 2015, 02:17:49 AM by Mike Lobanovsky »

Peter

  • Guest
Re: Mandelbrot
« Reply #5 on: October 15, 2015, 07:15:15 AM »
Hello,
here's another Mandel Julia test.
Code: [Select]
include "asm.inc"
window 1024,768,1  '0=fullscreen, be careful with windowsXP due to the desctop icons!

float r,i,x,y,z t
int c[32]

for i=0 to 32
    c[i] = RGB(i*8,i*14,i*8)
next

Sub Invers() 
d = x*x + y*y
if d=0 then
   r=100000
   x=r
   i=r
   y=r
else
   r=r/d
   x=r
   i=i/d
   y=i
end if
End Sub

Sub Inside(float d, fac, z, it)
i = Abs(Log(d))*fac
if z = it then
   r = GetR(c[i])
   g = GetG(c[i])
   b = GetB(c[i])
   color r,g,b
else
   color 0,0,0
end if
End Sub

Sub Julia_Init(float xss,yss,lt,rt,bm,tp,it,bb,inv,ins,fac,dx,dy,ci,cr)
dxpp = (rt-lt)/xss
dypp = (tp-bm)/yss
r=bm
for xr=0 to xss
    i = lt
for yr=0 to yss
    x = r
    y = i
    b = x
    m = y
    if inv=1 then Invers()
    xs = x*x
    ys = y*y
    while key(27)=0
    y = x*y
    Y = y+y-cr
    x = xs-ys-ci
    xs = x*x
    ys = y*y
    d = xs+ys
    if d < bb and z < it
       z +=1
    else
       exit while
    end if
    wend 
    if ins=0
       r = GetR(c[z])
       g = GetG(c[z])
       b = GetB(c[z])
       color r,g,b
    else
       Inside(d, fac, z, it)
    end if
    setpixel xr+dx, yr+dy
    z=0
    r=b
    i=m
    i += dypp
next
r += dxpp
redraw
next
End Sub

lt=-2 : rt=2
bm=-2 : tp=2
t = ticks
Julia_Init(1024,768,lt,rt,bm,tp,20,30.0,1,1,5,1,1,-0.32,-0.043)
wintext "TIME: " + (ticks()-t)/1000 + " SECONDS"
waitkey
WinExit


[attachment deleted by admin]

Mike Lobanovsky

  • Guest
Re: Mandelbrot
« Reply #6 on: October 15, 2015, 01:25:11 PM »
Hi Peter,

Thanks for that one! On my PC, it takes exactly 2.25 seconds to draw. My equivalent rewritten in FBSL draws in 0.296 seconds with 20 iterations as your original, and in exactly 0.5 seconds, with 100 iterations for much richer output. FBSL also uses my hand-written approximation of log() function that's about 7 times faster than the one from the standard C library.

[UPD] On revisiting the code in the morning, I was able to further reduce the 100-iteration version's exec time to 0.391 seconds by eliminating a good deal of redundant SetPixelV() calls with Julia's functionality intact. Code updated, snapshot and executable reuploaded.

Please see below for the FBSL sources, snapshot, and executable.

Code: OxygenBasic
  1. #Define width 1024
  2. #Define height 768
  3.  
  4. Type RECT
  5.   Left As Long
  6.   Top As Long
  7.   Right As Long
  8.   Bottom As Long
  9. End Type
  10.  
  11. Type BITMAPINFOHEADER Align 2
  12.   biSize As Long
  13.   biWidth As Long
  14.   biHeight As Long
  15.   biPlanes As Integer * 16
  16.   biBitCount As Integer * 16
  17.   biCompression As Long
  18.   biSizeImage As Long
  19.   biXPelsPerMeter As Long
  20.   biYPelsPerMeter As Long
  21.   biClrUsed As Long
  22.   biClrImportant As Long
  23. End Type
  24.  
  25. Dim rc As RECT, bi As BITMAPINFOHEADER
  26.  
  27. With bi
  28.   .biSize = SizeOf(BITMAPINFOHEADER)
  29.   .biWidth = width
  30.   .biHeight = height
  31.   .biPlanes = 1
  32.   .biBitCount = 24
  33. End With
  34.  
  35. Dim %c[32]
  36. For Dim i = 0 To 32
  37.   c[i] = RGB(i * 8, i * 14, i * 8)
  38. Next
  39.  
  40. Dim lt = -2, rt = 2
  41. Dim bm = -2, tp = 2
  42. Dim hBmp = CreateDIBSection(NULL, @bi, 0, NULL, NULL, 0)
  43. Dim hDC = GetDC(ME), hMemDC = CreateCompatibleDC(NULL)
  44. Dim hOldBmp = SelectObject(hMemDC, hBmp): SetStretchBltMode(hDC, 4) ' HALFTONE
  45. Dim t = GetTickCount()
  46.  
  47. Julia(hMemDC, lt, rt, bm, tp, 100, 30, TRUE, TRUE, 5, 1, 1, -0.32, -0.043)
  48.  
  49. FbslSetText(ME, "-= GDI Julia =-  TIME: " & (GetTickCount() - t) / 1000 & " SECONDS")
  50. Resize(ME, 0, 0, width, height)
  51. Center(ME)
  52. Show(ME)
  53.  
  54. Begin Events
  55.   If CBMSG = &H111 AndAlso CBWPARAM = 2 Then ' WM_COMMAND && Esc
  56.    PostMessage(ME, &H10, 0, 0) ' WM_CLOSE
  57.  ElseIf CBMSG = &HF OrElse CBMSG = &H5 Then ' WM_PAINT || WM_SIZE
  58.    GetClientRect(ME, @rc)
  59.     ValidateRect(ME, NULL)
  60.     StretchBlt(hDC, 0, 0, rc.Right, rc.Bottom, hMemDC, 0, 0, width, height, &HCC0020) ' SRCCOPY
  61.    Return 0
  62.   ElseIf CBMSG = &H10 Then ' WM_CLOSE
  63.    DeleteObject(SelectObject(hMemDC, hOldBmp))
  64.     DeleteDC(hMemDC)
  65.     ReleaseDC(ME, hDC)
  66.   End If
  67. End Events
  68.  
  69. DynC Julia(%dc, %lt, %rt, !bm, !tp, %it, !bb, %inv, %ins, !fac, !dx, !dy, _
  70.           !ci, !cr, %clr = @c, !xss = width, !yss = height)
  71.  
  72.   double fabs(double);
  73.   void __attribute__((stdcall)) SetPixelV(int, int, int, int);
  74.  
  75.   float mylog(float x) { // 7 times faster than log(), accurate to two decimal places
  76.     union {float f; int i;} vx = {x};
  77.     float y = vx.i;
  78.     y *= 1.0 / (1 << 23);
  79.     return 0.69314718f * (y - 126.94269504f);
  80.   }
  81.  
  82.   void main(int hdc, int lt, int rt, float bm, float tp, int it, float bb, int inv, int ins,
  83.             float fac, float dx, float dy, float ci, float cr, int c[], float xss, float yss)
  84.   {
  85.     int xr, yr, z, ix, iy;
  86.     float xs, ys, b, m;
  87.     float dxpp = (rt - lt) / xss;
  88.     float dypp = (tp - bm) / yss;
  89.     float r = bm, i, x, y, xs, ys, d, b, m;
  90.    
  91.     for (xr = 0; xr < xss; xr++) {
  92.       i = lt;
  93.       for (yr = 0; yr < yss; yr++) {
  94.         b = x = r; m = y = i;
  95.         if (inv) { // inverse()
  96.           d = x * x + y * y;
  97.           if (d == 0.0f)
  98.             x = y = i = r = 100000.0f;
  99.           else {
  100.             x = r /= d; y = i /= d;
  101.           }
  102.         }
  103.         xs = x * x; ys = y * y;
  104.         while (1) {
  105.           y += (y *= x) - cr;
  106.           x = xs - ys - ci;
  107.           d = (xs = x * x) + (ys = y * y);
  108.           if (d < bb && z < it)
  109.             z++;
  110.           else
  111.             break;
  112.         }
  113.         if (!ins) // inside()
  114.           SetPixelV(hdc, xr + dx, yr + dy, c[z]);
  115.         else if (z == it)
  116.           SetPixelV(hdc, xr + dx, yr + dy, c[(int)(fabs(mylog(d)) * fac)]);
  117.         z = 0; r = b;
  118.         i = m + dypp;
  119.       }
  120.       r += dxpp;
  121.     }
  122.   }
  123. End DynC

[attachment deleted by admin]
« Last Edit: October 16, 2015, 12:36:04 AM by Mike Lobanovsky »

Peter

  • Guest
Re: Mandelbrot
« Reply #7 on: October 16, 2015, 01:08:00 AM »
Hi Mike,

indeed. very fast!  :D
SetDIBSection is a tad unfair, we don't have the same conditions.   ;D

[attachment deleted by admin]

Peter

  • Guest
Re: Mandelbrot
« Reply #8 on: October 16, 2015, 01:24:36 AM »
Hi Mike,

found this at basicprogramming.org
Code: [Select]
'vonoroi.bas for SmallBASIC 2015-10-14 MGA/B+
'modified bluatigro code from JB 2015-09-07 vonoroi diagram

include "asm.inc"
window 800,600,1

int sq=800,s2=sq/2,points=240
int x[240],y[240],kl[240]

single t = ticks
for i=0 to points
    x[i] = rand(32,sq)
    y[i] = rand(32,sq)
    g=127-127*(abs(s2-x[i])/s2)+127-127*(abs(s2-y[i])/s2)
    kl[i] = RGB(255-x[i]/sq*255,g,y[i]/sq*255)
next

for xx=0 to sq
for yy=0 to sq
    d = 307201
for i=0 to points
    a=x[i]-xx: b=y[i]-yy
    q=a*a+b*b
    if q<d then
       d=q
       kkl = i
    end if   
next
PutPixel xx, yy, kl[kkl]
next
next
WinText "TIME: " + (ticks()-t)/1000 + " SECONDS"

waitkey
winExit



[attachment deleted by admin]

Peter

  • Guest
Re: Mandelbrot
« Reply #9 on: October 16, 2015, 01:39:53 AM »
Here's a SDL version.
Code: [Select]
'vonoroi.bas for SmallBASIC 2015-10-14 MGA/B+
'modified bluatigro code from JB 2015-09-07 vonoroi diagram

actually, his name is Voronoy!

include "sdl.inc"
window 800,600,1

int sq=800,s2=sq/2,points=240
int x[240],y[240],kl[240]

single t = ticks
for i=0 to points
    x[i] = rand(32,sq)
    y[i] = rand(32,sq)
    g=127-127*(abs(s2-x[i])/s2)+127-127*(abs(s2-y[i])/s2)
    kl[i] = RGBA(255-x[i]/sq*255,g,y[i]/sq*255,255)
next

for xx=0 to sq
for yy=0 to sq
    d = 307201
for i=0 to points
    a=x[i]-xx: b=y[i]-yy
    q=a*a+b*b
    if q<d then
       d=q
       kkl = i
    end if   
next
PutPixel xx, yy, kl[kkl]
next
next
SetCaption "TIME: " + (ticks()-t)/1000 + " SECONDS"
waitkey
winExit

[attachment deleted by admin]
« Last Edit: October 16, 2015, 09:52:40 AM by Peter »

Mike Lobanovsky

  • Guest
Re: Mandelbrot
« Reply #10 on: October 16, 2015, 03:25:13 AM »
Hi Peter,

SetDIBSection is a tad unfair, we don't have the same conditions.   ;D

DIB section isn't critical in this context as long as we're using anything SetPixelV-/PSet-/setpixel-compatible. We could as well have used CreateCompatibleBitmap instead and it wouldn't change the speed in any way. SetPixelV is reasonably fast for BASIC when drawing static images to an invisible memory DC for picture persistency but becomes dramatically slow for per-pixel draws directly on screen. Subsequent bit-blitting or stretch-blitting of the entire DC in response to the WM_PAINT/WM_RESIZE messages is very fast if not artificially looped, because it is harware-assisted on all modern GPUs.

However, if we would try and store colors directly into the bitmap's pixel array instead of PSet-ing pixels in the in-memory canvas like we did, we would have to use the DIB section only. Direct writes to the DIB could speed up the program by yet another 30% but they would also require a different [x, y] order of addressing pixels and the loop code would have to be changed.

Quote from: Peter's new Big_Julia
TIME: 0.516 SECONDS

That's nice! That's what I had been getting before I eliminated the SitPixelV calls that would draw black pixels over the initially black background in the newly created memory context. :)

I'm also going to give your Voronoi bananas a try in the evening, but now it's time for a short nap.


Aurel

  • Guest
Re: Mandelbrot
« Reply #11 on: October 16, 2015, 10:02:39 AM »
Hey Mike
Your julia100 get 2.068 sec on my old XP 1.7Ghz machine...is that ok?
 ;D

Mike Lobanovsky

  • Guest
Re: Mandelbrot
« Reply #12 on: October 16, 2015, 10:17:49 AM »
Hi Aurel,

How long does Peter's latest big_julia take on your PC?

Peter

  • Guest
Re: Mandelbrot
« Reply #13 on: October 16, 2015, 10:26:50 AM »
Perhaps but none of us will live to see that, so why bother advertizing in a technical thread?

advertising
Sorry Mike, but I am English learner!  ;D

Mike Lobanovsky

  • Guest
Re: Mandelbrot
« Reply #14 on: October 16, 2015, 10:54:41 AM »
« Last Edit: October 16, 2015, 11:08:49 AM by Mike Lobanovsky »