[Please do not mail me a copy of your followup] I'm interested in talking with other folks about the best way to generalize various fractal algorithms in the style of the standard library: iterators, containers and algorithms. If there are any libraries that already exist that apply to this sort of code, I'd love to hear about them. Most fractals consist of iterating a sequence for a point in the complex plane. The iterating sequence is parameterized for each point in the complex plane and some characteristic of the iterating sequence is used to color that point. The result is an image. For instance, the classic Mandelbrot fractal is: z <- z^2 + c z0 := 0 The complex number 'c' changes for each point in the complex plane. This sequence is iterated until norm(z) > 4, or some maximum iteration count is reached. Pixels in the image correspond to each value of 'c' and are typically colored according to the number of iterations needed until norm(z) exceeded 4. There are a couple things that make for interesting variation in the algorithms used in rendering fractal images: - choice of formula to iterate in the complex plane - parameters for deciding the "bailout" criteria: - iteration counts - predicate function on z, i.e. [](complex<float> z){ return norm(z) > 4; } Iterating until the bailout test passesis the simplest mechanism, but there are gains to be made by applying more sophisticated mechanisms. A significant optimization is to look for cycles in the sequences. With very large iteration counts, detecting cycles in the iteration can save significant amounts of time. The cycle detection algorithm is basically from Knuth and isn't specific to fractals but ends up being mixed together with the Mandelbrot set computation in most programs. This seems to be another opportunity to use generic composition in order to apply this orthogonal feature to any iterated formula in the complex plane. I did a little experiment casting the Mandelbrot formula as an iterator and got the code at the end of this message. Ideally I'd like to explore approaches that result in code that is ignorant of the underlying representation of a complex number in order to have algorithms that can be applied to arbitrary precision numeric types. My little experiment works, but even just looking at this little experiment makes me think that the mandelbrot iterator is doing too much: it is performing the z <- z^2 + c business, but it is also doing the bailout test and these seem to be separate responsibilities. #include <complex> #include <iostream> #include <iterator> #include <vector> class mandelbrot : public std::iterator< std::input_iterator_tag, std::complex<float>, unsigned> { public: mandelbrot() : c_{}, iter_{}, max_iter_{}, bailout_{}, z_{}, bailed_out_{true} { } mandelbrot(std::complex<float> const& c, int max_iter = 256, float bailout = 2.0f) : c_{c}, max_iter_{max_iter}, bailout_{bailout*bailout}, iter_{0}, z_{0.f, 0.f}, bailed_out_{false} { } mandelbrot& operator++() { z_ *= z_; z_ += c_; bailed_out_ = ++iter_ > max_iter_ || std::norm(z_) >= bailout_; return *this; } bool operator==(mandelbrot const& rhs) const { if (bailed_out_ && rhs.bailed_out_) { return true; } if (bailed_out_ != rhs.bailed_out_) { return false; } return z_ == rhs.z_; } std::complex<float> const& operator*() const { return z_; } private: std::complex<float> c_; int const max_iter_; float const bailout_; int iter_; std::complex<float> z_; bool bailed_out_; }; bool operator!=(mandelbrot const &lhs, mandelbrot const &rhs) { return !(lhs == rhs); } int main() { float const r_min = -3.0f; float const r_max = 1.0f; float const i_min = -2.0f; float const i_max = 2.0f; std::vector<float> imag; imag.resize(24); for (int y = 0; y < 24; ++y) { imag[y] = i_min + (i_max - i_min)*(23 - y)/24.0f; } std::vector<float> real; real.resize(80); for (int x = 0; x < 80; ++x) { real[x] = r_min + (r_max - r_min)*x/80.0f; } std::string printables; for (int i = 32; i < 128; ++i) { printables.push_back(char(i)); } for (int y = 0; y < 24; ++y) { float const im = imag[y]; for (int x = 0; x < 80; ++x) { float const re = real[x]; mandelbrot end; mandelbrot it{std::complex<float>{re, im}}; unsigned count = std::distance(it, end); std::cout << printables[count % printables.size()]; } std::cout << '\n'; } } -- "The Direct3D Graphics Pipeline" free book http://tinyurl.com/d3d-pipeline The Computer Graphics Museum http://computergraphicsmuseum.org The Terminals Wiki http://terminals.classiccmp.org Legalize Adulthood! (my blog) http://legalizeadulthood.wordpress.com