[R] Re: Mandelbrot set and C code
ucgamdo@ucl.ac.uk
ucgamdo at ucl.ac.uk
Wed Oct 1 15:09:36 CEST 2003
I decided to take on the 'proper' solution to calculate the
Mandelbrot set in R, i.e. to do the raw calculations in C and then link
that code with R. I thought it would be a hard task, but I was pleasantly
surprised when I saw how easily was to write the bit of C code (I am not a
C programmer myself!) and how even easier was to link it with R using .C()
and R CMD SHLIB. If you are interested (or anyone else) here is the (very
fast) code to calculate the mandelbrot set, I hope you enjoy playing with
it as much as I did.
# BEGINNING OF R CODE: ################################################
#######################################################################
# Function to calculate the Mandelbrot set. This function calls a #
# C routine in order to perform the calculations faster. #
# #
# Written by Mario dos Reis. September 2003 #
#######################################################################
mandelbrot <- function(x = c(-3, 1), # x limits
y = c(-1.8, 1.8), # y limits
nx = 600, # x resolution
ny = 600, # y resolution
iter = 20) # maximun number of iterations
{
xcoo <- seq(x[1], x[2], len = nx) # x coordinates
ycoo <- seq(y[1], y[2], len = ny) # y coordinates
set = numeric(nx*ny) # this will store the output of
# the C routine
# This is the call to the C function itself
the.set = .C("mandelbrot",
xcoo = as.double(xcoo),
ycoo = as.double(ycoo),
nx = as.integer(nx),
ny = as.integer(ny),
set = as.integer(set),
iter = as.integer(iter))$set
# Create a list with elements x, y and z,
# suitable for image(), persp(), etc. and return it.
return(list(x = xcoo, y = ycoo, z = matrix(the.set, ncol = ny, byrow = T)));
}
# END OF R CODE. #######################################################
/* BEGINNING OF C CODE: ***********************************/
#include <stdio.h>
/***********************************************************
* Function to evaluate a series of complex numbers that *
* form a subset of the Argand plane, to see if *
* they belong to the Mandelbrot set. *
* This function has been written with the intention of *
* linking it to some R code in order to create a *
* practical way to play with the set. *
* *
* Written by Mario dos Reis. September 2003 *
* *
***********************************************************/
void mandelbrot(double *xcoo, double *ycoo, int *nx,
int *ny, int *set, int *iter)
/* 'xcoo' and 'ycoo' are the x and y coordinates respectively
* of all the points to be evaluated.
* 'nx' and 'ny' number of divisions along the x and y axes
* 'set' will store the practical output of the function
* 'iter' is the maximun number of iterations
*/
{
int i, j, k;
double z[2], c[2], oldz[2];
for(i = 0; i < *nx; i++) {
for(j = 0; j < *ny; j++) {
/* initialise the complex point to be tested
* c[0] (or z[0]) is the real part
* c[1] (or z[1]) is the imaginary part
*/
c[0] = xcoo[i]; c[1] = ycoo[j];
z[0] = 0; z[1] = 0;
for(k = 1; k < *iter + 1; k++) {
oldz[0] = z[0]; oldz[1] = z[1];
// the mandelbrot mapping z -> z^2 + c
z[0] = oldz[0]*oldz[0] - oldz[1]*oldz[1] + c[0];
z[1] = 2 * oldz[0]*oldz[1] + c[1];
if((z[0]*z[0] + z[1]*z[1]) > 4) {
break;
}
}
/* fills the set vector
* notice the trick 'i * (*ny) + j' to find
* the apporpiate position of the output in the
* vector set. The R function will take care of
* transforming this vector into a suitable matrix
* for plotting, etc.
*/
set[i * (*ny) + j] = k;
}
}
return;
}
/* END OF C CODE. ******************************************
To anyone interested in trying out these functions, just save the C code
into a file, say 'mandelbrot.c', to be able to use the C code with R a
shared object needs to be created. This is done simple running the
following command from the directory where 'mandelbrot.c' is stored (note
this code was tried in Linux, for Win users I think the instructions are
similar but you need you check it out):
R CMD SHLIB mandelbrot.c
This will create the file 'mandelbrot.so' in that same directory.
Start R from this directory and paste the above R code into the console. Type
> dyn.load("mandelbrot.so")
in order to load the shared object. Now you are ready to appreaciate in
full detail the set.
> image(mandelbrot(), col = c(heat.colors(49), "black"))
This function is many orders of magnitude faster than the one I had
published previously
(see http://maths.newcastle.edu.au/~rking/R/help/03b/3390.html)
that function also had a slight error, it tested (z[1] + z[2])^2 > 4
instead of z[1]^2 + z[2]^2 > 4 which creates the asymmetric result you
pointed out in the outer regions of the set. The C code is also at
extremely fast with high zooms and iteratons, try:
xcoo = c(-1.1854735004165451, -1.1854647240344973)
ycoo = c(-0.3057632375298113, -0.3057544611477636)
frac <- mandelbrot(x = xcoo, y = ycoo, iter = 2000) # took ~ 12s in P4 2.2GHz
image(frac, col = c(heat.colors(49), "black"))
This is at least 19x faster than the code you suggested me. (same settings
took ~ 230s in the same machine).
Regards,
Mario.
At 10:21 18/09/03 +0200, you wrote:
>
> Mario> Well, I started playing with fractals in R, and wrote
> Mario> a function to generate de Mandelbrot set, which might
> Mario> be of interest to some people
>
> Mario>
###################################################################
> Mario> # Mandelbrot set
> Mario>
###################################################################
>
> Mario> mandelbrot <- function(x = c(-3.0, 1.0), # x coordinates
> Mario> y = c(-1.8, 1.8), # y coordinates
> Mario> b = .05, # by 'steps'
> Mario> iter = 20) # maximum number of
iterations
>
> <..................>
> <..................>
> <..................>
>
>Well, only a bit more than year ago I had posted my version of
>mandelbrot() and a drawing function, see
>
> http://finzi.psych.upenn.edu/R/Rhelp02a/archive/5898.html
>
> [ I have a slightly updated version of the R code, that is now
> available as ftp://stat.ethz.ch/U/maechler/R/Mandelbrot.R
> ]
>
>which is an order of magnitude more efficient than yours (22 x
>faster for your "b = 0.01", *1), *2)
>and with an iterImage() function for drawing its result in
>several (simple) color schemes.
>I agree that the proper solution would *definitely* use C code!
>
>Your idea of using persp() -- while seen frequently in fractal books,
>is new ``for R'' however, and nice! Thank you.
>
>
>Regards,
>Martin
>
>*1) mainly because I only have the iteration loop
> and do the rest vectorized, but also because I have a check
> for symmetry and make use of it when appropriate
>
>*2) For speed comparison, note that I use 10 x more iterations
> and a much higher resolution by default
>
>*3) When I do the comparison, I see that your function's result
> slightly differs from mine -- not for the Mandelbrot set
> itself, but for the very ``outer points'' (the dark orange
> ones in the image plot). Your result is
> asymmetric and hence can't be quite correct.
>
>
More information about the R-help
mailing list