The FBcomplex lib from FreeBasic

[RobbeK]

Very easy & smooth integration of complex numbers in a FB source

----------------------

#include "fbcomplex.bi"

screenres 800,800,32

screencontrol 100,10,10
screencontrol 101,"Newton Raphson power 7"

 
dim shared as single x,y
dim shared as integer i,j,res

function newton (z as complex) as integer
    dim i as integer
    dim zz as complex
    dim residu as single = 100
    do while (residu > 0.0000001) and (i<40)
        zz=z-(z^7-1)/(7*z^6)
        residu = abs(cabs(z)-cabs(zz))
        z=zz
        i+=1
    loop
    return i
end function

sub main()
    'screenlock()
    i=0 : j=0
    for x=-2.0017 to 2.0 step 4/800
        j=0 : i+=1
        for y=-2.0017 to 2.0 step 4/800
          j+=1
          res=6*newton (Cmplx (x,y))
           pset (i,j),rgb(res-255,res,255*abs((400-i)*(400-j))/160000)
          next
    next
    'screenunlock()
    sleep 100000
end sub


main()

---------------------------------------------------

a 7th power upto 40 it deep, you need an arm length of numbers/symbols to write this as a f(x,y)

best Rob  source and exec. included in the 7zip file

addendum : that 2.0017 to avoid the 0 pole - the function is not defined there   because  x/0

.

[Mike Lobanovsky]

Hi Rob,

The fractal is spectacular and the code runs pretty fast though I think it can be made even faster if PSet() is avoided in favor of modifying the main window's DIB pixel data array directly and then bitblitting the modified DIB from memory onto the screent window canvas in one swoop.

Still I may be having some infantile complex which alerts my psyche to everything with 5 extremities and more that can be potentially alive and harmful while anything less can't. This somehow overlooks worms and serpents but my only reasonable guess is that I grew up in the south and consequently am not afraid of either.

[RobbeK]

Hi Mike, 

Thanks  , yes it will run faster , but (as you can see in the code - but while not bitblitting ) i changed back from a back- , frontbuffer & flipping (which is also faster) to direct plotting. With either of these 2 methods the screen will be empty for some seconds which looks rather frustrating 
(the gain is very low against the huge processing time -- FB does this about the double speed comperated with Clozure Common Lisp (though I did no type declarations which speed up things somewhat )

attached the 7z+5z buffered graphics.

I will write the code with a Graph. interface for any power combination upto 7   (don't worry , during my quarantaine I'm 100% cured from JAVA   

best Rob


.

[Mike Lobanovsky]

Very fine-looking, thanks a lot, Rob!

Actually, what takes a few seconds on your PC might take only a fraction of a second on others' ones...

[Aurel]

i am wondering is possible to build this without this complex library directly in o2

[Mike Lobanovsky]

Hi Aurel,

Sure it can be done, especially if you only need a couple of functions that operate on complex numbers. Here is the link to the Russian Visual Basic site where this piece of VB code is an extended library of math functions including complex numbers. The functions that emulate complex numbers have a cx prefix. You can see that a complex number is actually a structure with two double-precision numbers, and all you have to do is just perform two math operations on two double values instead of one.

[Aurel]

thanks MIke
i see : kompljeksne čisla
i will try...

[RobbeK]

Hi Aurel,

No computer (even no human) can do even something simple as adding two complex numbers ... it has to be split in a real and imaginary part. 
The basics are --  complex numbers written as Z1=(a,b) and Z2=(c,d)

(a,b)+(c,d) = (a+c , b+d)
(a,b)-(c,d) = (a-c , b-d)   because m*(c,d) = (mc,md) consequently -(c,d) = -1*(c,d) = (-c,-d)

(a,b)*(c,d) = (ac - bd , bc + ad) consequently    (a,b)² = (a²-b², 2*a*b)

(a,b)/(c,d) = m*(ac+bd , bc-ad) where m=c²+d²  a real number.  (it is done by multiplying with the conjugate complex number )

you will have to do Z*Z*Z*Z*Z*Z for Z^6   ,  the FBmath module does A=Z*Z*Z en then Z^6 = A*A as an example. 


best to declare :

type complex
  r as single
  i as single
end type

z as complex    gives z.r (real part)  and z.i (imag part )

with this you can do this in any language.

oops , forgot the magnitude (modulus)
that's the length of the vector

by Pythagoras rule   cabs(a,b) = sqrt(a*a + b*b)

best Rob
 
addendum : in many math programs , the sqrt is just skipped (because expensive)  , you better square your criterium .

x < sqrt 5  is the same as x² < 5 if x > 0   this is done in all those Mandelbrot progs ..

[RobbeK]

Hi all ,

good news -- reduced the if then's to 7 per cycle (by using simul operations on the nominator and denominator ),  speed about upto 5x faster.

bitblitting now ?  -- or OGL and going 3D ??

best, Rob

also attached a version that takes complex coefficients (more fun) -- and same with buffered grahics

.

[Aurel]

best to declare :

type complex
  r as single
  i as single
end type

hi rob..
so if I wish to do this in o2 i must use UDT for complex numbers and try to create such a fractal
right?

[jack]

hello RobbeK
did you modify fbcomplex ?
if so, care to share the code ?

[RobbeK]

Hi,

@ Aurel :  I'll write the O² functions for you , or do you wish a DLL with these basics ?

@ Jack : no , but the nominator and denominater are updated at the same time, something as :

for every power of the polynomial :
  if the modulus of the coefficient is non zero then update the nominator and the denominator

the optimum would be a pre processor as the one I wrote in Scheme

(define optimized-formula (simplify the-formula))      removing the cabs(z)=0 therms (not sure this is the correct english - in maths i'm more used to german)

(eval (cons operator optimized-formula))  execute the data to code switch   (Lisp is such powerful,  but I only can get it (SBCL / Clozure CL) half the speed of O² and lacking something superb as the Firefly RAD )

There is one further optimization because        f' (z^n) = nz^(n-1) we don't have to expand the f'(z) as a power but convert it into a more simple division from f(z)  ;  in the next and final version I'll include a scaling factor for a closer or wider look around the origin (0,0) 

best , Rob

[Aurel]

@ Aurel :  I'll write the O² functions for you , or do you wish a DLL with these basics ?

hi rob...
thanks 
show us functions better than dll..or maybe someone want dll ...its up to you

[RobbeK]

Hi,

Quickly coded in FB (should be identical with O² sofar ?)   it's included in the attached + a demo (doing a Mandelbrot with complex numbers --  much easier when using complex calculus )

'--------------------------------------------------------------------------------------------------------

type complex
    r as single
    i as single
end type

 


function cmplx (x as single , y as single ) as complex
    dim z as complex
    z.r=x : z.i=y
    return z
end function

sub cprint(z as complex)
   print str$(z.r)+ " + i " +str$(z.i)
end sub


function Cadd (x as complex , y as complex) as complex
   dim z as complex
   z.r=x.r+y.r : z.i=x.i+y.i
   return z
end function

function Csub (x as complex , y as complex) as complex
    dim z as complex
    z.r=x.r-y.r : z.i=x.i-y.i
    return z
end function

function Cmul (x as complex , y as complex) as complex
    dim z as complex
    z.r=x.r*y.r-x.i*y.i : z.i=x.r*y.i+y.r*x.i
    return z
end function

function Csqr (x as complex) as complex
    return cmul(x,x)
end function

function modulus(x as complex) as single
    return sqr (x.r*x.r+x.i*x.i)
end function

function conjugate(x as complex) as complex
    dim z as complex
    z.r=x.r : z.i=-x.i
    return z
end function

function scalar(m as single, x as complex) as complex
    dim z as complex
    z.r=m*x.r : z.i=m*x.i
    return z
end function   

function cinvers(x as complex) as complex
    return scalar(1/(x.r*x.r + x.i*x.i),conjugate(x))
end function

'------------------------------------------------------------------------------------------------------

there is no complex division , because in fact it does not exist -- you have do a multiplication with the inverse complex number

A / B  is identical with A * cinvers(B)


good luck !

Rob

.

[RobbeK]

Final version




.