Oxygen Basic
Information => Open Forum => Topic started by: RobbeK on May 08, 2015, 07:18:19 AM
-
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
.
-
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. :)
-
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 8)
best Rob
.
-
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... :)
-
i am wondering is possible to build this without this complex library directly in o2
-
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 (http://bbs.vbstreets.ru/viewtopic.php?f=28&t=45329) 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. :)
-
thanks MIke
i see : kompljeksne čisla :D
i will try... ;)
-
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 ..
-
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
.
-
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?
-
hello RobbeK
did you modify fbcomplex ?
if so, care to share the code ?
-
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 : 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
-
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
.
-
Final version
.
-
RobbeK
I experimented with your fractal program a little, I copied fbcomplex.bas and fbcomplex.bi to fbcomplex-s.bas and fbcomplex-s.bi and edited them to use single instead of double
then compiled using: fbc -c -gen gcc -Wc -Ofast fbcomplex-s.bas
and compiled your fractal program using the same switches, it's roughly twice as fast, but there is a visible dot in the yellow disk
.
-
the included library has issues, when using factor 5 it's distortion is severe.
-
Hi Jack,
a lot of thanks
don't worry about that black hole on 0,0 -- it should be there , it's a (near) by zero division ...
"distortion is severe" -- maybe you have to change the loop escape mechanism (I wrote this with double floats in mind)
(residu > 0.0000001) replace-> (residu > 0.00001). it's possible with such a small value , the single floats become "bouncing"
The idea of scaling the coefficients is not so bright , I'll include a coordinate system scaling/transformation.
best, Rob
-
RobbeK
I changed the resedue as suggested and it fixed the strange result when using scaling factor 5 :)
-
hi Rob
Your example of Mandelbrot is interesting because fractal is created from
outside to inside ..nice effect.. :)
-
Thanks Aurel
Hi Jack, yes I suspected that value was too small for single floats. I added the viewport width now ....
now embedding Newlisp into freebasic , could be vary handy sometimes
best , Rob
.