I wrote a simple script for a Schroder fractal (a distant cousin of Newton
fractals) of a fourth-order polynomial:
function SchroderFourth(a0::Float64, a1::Float64, a2::Float64, a3::Float64,
a4::Float64, c::Complex{Float64})
x = c
x2 = c
M = 100
R = 10.0^-6
a = 0.1
for j in 1:M
x2 = x
f = a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4
fP = a1 + 2*a2*x + 3*a3*x^2 + 4*a4*x^3
fPP = 2*a2 + 6*a3*x + 12*a4*x^2;
x = x - f*fP/(fP^2 - fPP*f)
if abs(x - x2) < R
return exp(-a*(j - (log(abs(x - x2))-log(R))/log(R)) + angle(x)im)
end
end
return 0.0 + 0.0im
end
I tested it using @time as follows:
@time arr = [SchroderFourth(1.0,2.0,3.0,4.0,5.0,x+y*im) for x in
linspace(-2,2,1000), y in linspace(-2,2,1000)];
resulting in
elapsed time: 3.371773428 seconds (2239441280 bytes allocated, 48.23% gc
time)
So it's allocating 2.2GB throughout the computation, and half the time is
gc time. I have heard that sometimes this is caused by type instability,
but the SchroderFourth function is type-stable (always returns
Complex{Float64}). Am I committing any "great sins" in my implementation?
Given that the array has 1M computations, and each computation loops at
most M=100 loops, I suppose it's possible that gigabytes of allocation
could occur if a memory address was gobbled up after each iteration of the
loop. In comparison, running the same computation on the (simpler) "mandel"
fractal on the JuliaLang website gives much less memory allocation (8MB):
function mandel(z)
c = z
maxiter = 80
for n = 1:maxiter
if abs(z) > 2
return n-1
end
z = z^2 + c
end
return maxiter
end
@time arr = [mandel(x+y*im) for x in linspace(-2,2,1000), y in
linspace(-2,2,1000)];
giving
elapsed time: 0.123358763 seconds (8016176 bytes allocated)
Any clues as to the large disparity in memory allocation between the two
functions? I assume I'm doing something dumb that's causing memory to be
allocated and gc'ed after each iteration, but am having trouble finding
what is causing it.
