Complex Number Library
This is a work in progress collection of operators and functions for complex numbers. Mostly in the form of SUBs and FUNCTIONs at the end of the code. The demo shows the library in use in several complex number topics with domain coloring plots as well as a simple example of how to plot the Mandelbrot set at the end.
This is a work in progress collection of operators and functions for complex numbers. Mostly in the form of SUBs and FUNCTIONs at the end of the code. The demo shows the library in use in several complex number topics with domain coloring plots as well as a simple example of how to plot the Mandelbrot set at the end.
Code: (Select All)
defdbl a-z
const sw = 800
const sh = 600
dim shared pi
pi = 4*atn(1)
zoom = 140
dim as long i, j, k, xx, yy
screen _newimage(sw, sh, 32)
for i=0 to 9
'''plots
for yy=0 to sh
for xx=0 to sw
x = (xx - sw/2)/zoom
y = (sh/2 - yy)/zoom
select case i
case 0
u = x
v = y
pset (xx, yy), hrgb(u, v)
case 1
u = x
v = y
pset (xx, yy), checker(u, v)
case 2
'w = 1/z
cdiv u, v, 1, 0, x, y
'pset (xx, yy), hrgb(u, v)
pset (xx, yy), checker(u, v)
case 3
'w = sin(1/z)
'extra zoom
d = 0.35*(x*x + y*y)
if d<>0 then
u = sin(x/d)*cosh(-y/d)
v = cos(x/d)*sinh(-y/d)
else
u = 0
v = 0
end if
pset (xx, yy), hrgb(u, v)
'pset (xx, yy), checker(u, v)
case 4
'extra zoom
u = 0.56*x
v = 0.56*y
for j=0 to 14
uu = u*u - v*v + 0.35
vv = 2*u*v + 0.0
u = uu
v = vv
next
pset (xx, yy), hrgb(u, v)
'pset (xx, yy), checker(u, v)
case 5
cmul u, v, 1, 0, x - cos(2*pi/3), y + sin(2*pi/3)
cmul u, v, u, v, x - cos(2*pi/3), y - sin(2*pi/3)
cmul u, v, u, v, x - 1, y
'cdiv u, v, u, v, x - 1, y
pset (xx, yy), hrgb(u, v)
'pset (xx, yy), checker(u, v)
case 6
'CIF numerical integration
'f(z_0) = (2 pi i)^-1 int_C f(z)/(z - z0) dz
n = 35
uu = 0
vv = 0
for j=0 to n - 1
'C: z(t)
p = 1.5*cos(j*2*pi/n)
q = 1.5*sin(j*2*pi/n)
'f(z(t)):
cmul u, v, 1, 0, p - cos(2*pi/3), q + sin(2*pi/3)
cmul u, v, u, v, p - cos(2*pi/3), q - sin(2*pi/3)
cmul u, v, u, v, p - 1, q
'f(z)/(z - z0)
cdiv u, v, u, v, p - x, q - y
'dz/dt
cmul u, v, u, v, -1.5*sin(j*2*pi/n), 1.5*cos(j*2*pi/n)
if j = 0 or j = n - 1 then
uu = uu + 0.5*u
vv = vv + 0.5*v
else
uu = uu + u
vv = vv + v
end if
next
'dt
u = uu*2*pi/n
v = vv*2*pi/n
'1/(2 pi i)
cmul u, v, u, v, 0, -1/(2*pi)
pset (xx, yy), hrgb(u, v)
'pset (xx, yy), checker(u, v)
case 7
'extra zoom
x = x*0.5
y = y*0.5
p = 1
q = 0
for j=0 to 5
cmul uu, vv, 1, 0, -0.4, -0.18*(j - 2.1)
cmul p, q, p, q, x - uu - 0.2, y - vv
cdiv p, q, p, q, x - uu + 0.2, y - vv + 0.1
next
for j=0 to 2
cmul uu, vv, 1, 0, 0.4, -0.18*(j - 2.1) - 0.18*2.1/2
cdiv p, q, p, q, x - uu - 0.2, y - vv
cmul p, q, p, q, x - uu + 0.2, y - vv + 0.1
next
u = p
v = q
pset (xx, yy), grey(u, v)
case 8
'extra zoom
u = 0.66*x - 0.5
v = 0.66*y
x0 = u
y0 = v
for j=0 to 3
uu = u*u - v*v + x0
vv = 2*u*v + y0
u = uu
v = vv
next
'pset (xx, yy), hrgb(u, v)
pset (xx, yy), checker(u, v)
case 9
'extra zoom
u = 0.66*x - 0.5
v = 0.66*y
x0 = u
y0 = v
for j=0 to 70
uu = u*u - v*v + x0
vv = 2*u*v + y0
u = uu
v = vv
next
'pset (xx, yy), hrgb(u, v)
pset (xx, yy), checker(u, v)
end select
next
next
'''diagrams
select case i
case 0
_title "w = z polar contouring"
case 1
_title "w = z checkerboard"
case 2
_title "w = 1/z singularity"
case 3
_title "w = sin(1/z) essential singularity"
case 4
_title "Julia fractal"
case 5
_title "tri mass"
case 6
_title "Cauchy integral formula"
a = 0
x = 1.5*cos(a)
y = 1.5*sin(a)
circle (x*zoom + sw/2, sh/2 - y*zoom), 3, _rgb(255,255,0)
for a=0 to 2*pi step 2*pi/n
x = 1.5*cos(a)
y = 1.5*sin(a)
line -(x*zoom + sw/2, sh/2 - y*zoom), _rgb(255,255,0)
circle step(0,0), 3, _rgb(255,255,0)
next
case 7
_title "mutual inductance"
sleep
'extra zoom
m = zoom/0.5
'this diagram is a JB original
'left
a = -pi/2
x = sw/2/m + 0.2*cos(a) - 0.4
y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5/pi) - 2.1*0.18
for t=x to 0 step -0.001
circlef (x - t)*m, y*m, 1, _rgb(0,0,0)
next
for a = -pi/2 to 5*2*pi + 2*pi + pi/2 step 0.01
x = sw/2/m + 0.2*cos(a) - 0.4
y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5/pi) - 2.1*0.18
circlef x*m, y*m, 1, _rgb(0,0,0)
next
a = 5*2*pi + 2*pi + pi/2
x = sw/2/m + 0.2*cos(a) - 0.4
y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5/pi) - 2.1*0.18
for t=x to 0 step -0.001
circlef (x - t)*m, y*m, 1, _rgb(0,0,0)
next
'right
a = -pi/2
x = sw/2/m - 0.2*cos(a) + 0.4
y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5)/pi - 2.1*0.18 + 0.18*2.1/4
for t=0 to 1.5 step 0.001
circlef (x + t)*m, y*m, 1, _rgb(0,0,0)
next
for a = -pi/2 to 2*2*pi + 2*pi + pi/2 step 0.01
x = sw/2/m - 0.2*cos(a) + 0.4
y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5)/pi - 2.1*0.18 + 0.18*2.1/4
circlef x*m, y*m, 1, _rgb(0,0,0)
next
a = 2*2*pi + 2*pi + pi/2
x = sw/2/m - 0.2*cos(a) + 0.4
y = sh/2/m + 0.18*sin(a) + (a*0.18*0.5)/pi - 2.1*0.18 + 0.18*2.1/4
for t=0 to 1.5 step 0.001
circlef (x + t)*m, y*m, 1, _rgb(0,0,0)
next
case 8
_title "checkerboard Mandelbrot"
end select
sleep
next
_title "escape time example with nth order Mandelbrot"
zoom = 210
for n=2 to 6
'line (0, 0)-(sw, sh), _rgb(255,255,255), bf
for yy=0 to sh
for xx=0 to sw
x = (xx - sw/2)/zoom - 0.5
y = (sh/2 - yy)/zoom
u = 0
v = 0
for i=0 to 140
'f(z) = z^n + c
cexp u, v, u, v, n, 0
u = u + x
v = v + y
if u*u + v*v > 4 then exit for
next
if i>=140 then
pset (xx, yy), _rgb(0,0,0)
else
pset (xx, yy), _rgb(255,255,255)
end if
next
next
sleep
next
system
sub circlef(x as long, y as long, r as long, c as long)
dim as long x0, y0, e
x0 = r
y0 = 0
e = -r
do while y0 < x0
if e <=0 then
y0 = y0 + 1
line (x - x0, y + y0)-(x + x0, y + y0), c, bf
line (x - x0, y - y0)-(x + x0, y - y0), c, bf
e = e + 2*y0
else
line (x - y0, y - x0)-(x + y0, y - x0), c, bf
line (x - y0, y + x0)-(x + y0, y + x0), c, bf
x0 = x0 - 1
e = e - 2*x0
end if
loop
line (x - r, y)-(x + r, y), c, bf
end sub
function grey~&(x, y)
m = sqr(x*x + y*y)
a = (pi + _atan2(y, x))/(2*pi)
m = log(1 + 100*m)
'polar contouring
n = 16
mm = m*5000 mod 500
p = abs(a*n - int(a*n))
g = 1 - 0.0007*mm - 0.21*p
grey = _rgb(255*g, 255*g, 255*g)
end function
function checker~&(xx, yy)
if 1 then
x = xx
y = yy
else 'polar checkerboard
x = _atan2(yy, xx)/(pi/4)
y = sqr(xx*xx + yy*yy)
y = log(1 + 1000*y)
end if
z = abs(x - int(x)) xor abs(y - int(y))
if z then checker = _rgb(0,0,0) else checker = _rgb(255,255,255)
end function
function hrgb~&(x, y)
m = sqr(x*x + y*y)
a = (pi + _atan2(y, x))/(2*pi)
'm = log(1 + 1000*m)
r = 0.5 - 0.5*sin(2*pi*a - pi/2)
g = (0.5 + 0.5*sin(2*pi*a*1.5 - pi/2)) * -(a < 0.66)
b = (0.5 + 0.5*sin(2*pi*a*1.5 + pi/2)) * -(a > 0.33)
'polar contouring
n = 16
mm = m*500 mod 500
p = abs(a*n - int(a*n))
r = r - 0.0005*mm - 0.14*p
g = g - 0.0005*mm - 0.14*p
b = b - 0.0005*mm - 0.14*p
'cartesian shading
if 0 then
t = 0.03 'thickness
xx = abs(x - int(x)) < t or abs(-x - int(-x)) < t
yy = abs(y - int(y)) < t or abs(-y - int(-y)) < t
if xx or yy then
'if m > 1 then 'dont shade origin
r = r - 0.5
g = g - 0.5
b = b - 0.5
'end if
end if
end if
hrgb = _rgb(255*r, 255*g, 255*b)
end function
sub cmul(u, v, xx, yy, aa, bb)
x = xx
y = yy
a = aa
b = bb
u = x*a - y*b
v = x*b + y*a
end sub
sub cdiv(u, v, xx, yy, aa, bb)
x = xx
y = yy
a = aa
b = bb
d = a*a + b*b
u = (x*a + y*b)/d
v = (y*a - x*b)/d
end sub
sub cexp(u, v, xx, yy, aa, bb)
x = xx
y = yy
a = aa
b = bb
lnz = x*x + y*y
if lnz = 0 then
u = 0
v = 0
else
lnz = 0.5*log(lnz)
argz = _atan2(y, x)
m = exp(a*lnz - b*argz)
a = a*argz + b*lnz
u = m*cos(a)
v = m*sin(a)
end if
end sub
sub clog(u, v, xx, yy)
x = xx
y = yy
lnz = x*x + y*y
if lnz=0 then
u = 0
v = 0
else
u = 0.5*log(lnz)
v = _atan2(y, x)
end if
end sub
function cosh(x)
cosh = 0.5*(exp(x) + exp(-x))
end function
function sinh(x)
sinh = 0.5*(exp(x) - exp(-x))
end function
sub csin(u, v, xx, yy)
x = xx
y = yy
u = sin(x)*cosh(y)
v = cos(x)*sinh(y)
end sub
sub ccos(u, v, xx, yy)
x = xx
y = yy
u = cos(x)*cosh(y)
v =-sin(x)*sinh(y)
end sub
function factorial~&(n)
if n = 0 then
factorial = 1
else
factorial = n*factorial(n - 1)
end if
end function