Template Meta-Mandelbrot
Typically, we think of software as taking in data, transforming that data, and then presenting it somehow–say, printing to standard out.
But often times, some aspect of the calculation occurs exactly the same way, with the same result, no matter what data you supply to your program.
If you’ve been keeping up with the c++ template metaprogramming literature, then you know that it is an incredibly common requirement that you multiply some number by four factorial.
Obviously, the multiplication must be done anew each time, but there is no good reason to recalculate each and every time you run your program.
One way to deal with this would be to perform the calculation beforehand, and define a constant.
const unsigned int Factorial4 = 24;
Then along comes some wicked smaht sunofagun who figures out that you can use the c++ template system to do this for you at compile time.
template< int i >
struct factorial
{
enum { result = i * factorial<i - 1>::result };
};
template<>
struct factorial<0>
{
enum { result = 1 };
};
constexpr unsigned int Factorial4 = factorial<4>::result;
But why should we stop at calculating individual constants? In scientific computing (and I’ll bet quite a few other areas as well), it is common not only to have constants that are known at compile time, but also containers, say a vector or matrix.
As a case study, in this post we’ll walk through how you might go about calculating the Mandelbrot Set and storing it in an std::array at compile time using c++14.
Because the Mandelbrot Set, like template metaprogramming, is just so damn practical.
The Mandelbrot Set
For those of us who haven’t committed this most essential information to memory, a quick primer on the Mandelbrot Set.
The Mandelbrot Set is a set of complex numbers. Every complex number is either in the set, or it isn’t. If I give you a complex number, , you determine whether it is a member of the Mandelbrot Set by iteratively calculating , where in the first iteration (the tuple represents the real and imaginary components). If the absolute value of remains bounded as the number of iterations approaches infinity, then is a member of the Mandelbrot Set.
Let’s take two examples, and , and calculate a few iterations.
| Iteration | ||
|---|---|---|
| 0 | 1.0 | 1.0 |
| 1 | 0.0 | 2.0 |
| 2 | 1.0 | 5.0 |
| 3 | 0.0 | 26.0 |
| 4 | 1.0 | 677.0 |
| 5 | 0.0 | 458330.0 |
| 6 | 1.0 | 210066388901.0 |
We would say that is a member of the Mandelbrot Set, wherease is not.
Setting Up the Problem
We’re going to define a subset of the complex plane, and then for each member of that subset, we will calculate whether or not the pixel is a member of the Mandelbrot set.
Since the complex plane is , and an std::array is not, we’ll need to do a little bookkeeping to go from the array index to the position in space.
Let’s restrict our set to in both the real and complex directions, and lets sample pixels per unit.
We can calculate (at compile time, of course!) the stride (i.e., the number by which the array index increases upon advancing one row), and the total number of pixels in the image.
struct Mbrot
{
// Number of pixels per unit space
static constexpr size_t res = 16;
// Stride (Number of pixels across the entire image in either direction)
// 4*res+1, since we will span [-2, 2] in both the real (x) and complex (y) directions
static constexpr size_t str = 4*Mbrot::res+1;
// Total number of pixels in the image
static constexpr size_t pix = Mbrot::str*Mbrot::str;
};
Using the stride, we can calculate the row and column from the array index:
struct Mbrot
{
// ...
// Return the column number given an index into the flat array.
static constexpr size_t
Col(const size_t &i)
{
if (i >= Mbrot::pix) { throw std::out_of_range("Index out of range."); }
return i % Mbrot::str;
}
// Return the row number given an index into the flat array.
static constexpr size_t
Row(const size_t &i)
{
if (i >= Mbrot::pix) { throw std::out_of_range("Index out of range."); }
return i / Mbrot::str;
}
};
Note that if either of these functions are called at compile time with an out of range index, it will fail to compile, which is just what you want.
For example, if I try constexpr auto p = Mbrot::Col(Mbrot::pix);, I get a “subexpression not valid in a constant expression” error, referring to the throw.
Now, using these utilities, we can get directly from an index to the position in the complex plane:
struct Mbrot
{
// ...
// Return the position in the complex plane given an index into the flat array.
static constexpr std::complex<double>
Pos(const size_t &i)
{
if (i >= Mbrot::pix) { throw std::out_of_range("Index out of range."); }
return std::complex<double>(
(double(Mbrot::Col(i))-double(Mbrot::str/2))/double(Mbrot::res),
(double(Mbrot::Row(i))-double(Mbrot::str/2))/double(Mbrot::res)
);
}
};
Determining Membership in the Set
Having converted from index to complex position, we can use some straightforward constexpr recursion to determine whether a given complex number is a member of the Mandelbrot Set.
Note that, since we initialize with , the first iteration always returns .
Therefore, in our implementation we skip this and initialize to be .
To prevent going into an infinite loop, we’ll set a maximum number of iterations.
Also, you might find it odd that in is_inside_impl, I’ve written my own functions for adding, multiplying, and taking the absolute value of std::complex.
It turns out that this is necessary because, perplexingly, none of operator+, operator*, or std::abs are declared constexpr.
Luckily, all these functions are straightforward to write (see the source code if you’re interested in the definitions).
struct Mbrot
{
// ...
// Maximum number of iterations to determine convergence
static constexpr size_t its = 8;
// Determine whether a given position in the complex plane is in the Mandelbrot set.
static constexpr bool
IsInSet(std::complex<double> p)
{
return Mbrot::is_inside_impl(p, p, Mbrot::its);
}
private:
// Private implementation of the recursive function determining membership in the Mandelbrot set.
static constexpr bool
is_inside_impl(std::complex<double> initial,
std::complex<double> final,
size_t it)
{
return (0 == it || dv::abs(final) > 2) ?
dv::abs(final) < 2 : is_inside_impl(initial, dv::add(dv::mul(final,final), initial), it - 1);
}
};
Initializing the Container
We’re almost there. Given an index, we can calculate the corresponding position in the complex plane, and then determine whether that position is a member of the Mandelbrot Set. Intuitively, what we’d like to do is loop over the indices, assigning to the array as we go. Something like:
struct Mbrot
{
// ...
// Return the array.
static constexpr Array
GetArray()
{
auto a = Array();
for (size_t i = 0; i < Mbrot::pix; ++i) a[i] = Mbrot::IsInSet(Mbrot::Pos(i));
return a;
}
};
Unfortunately, operator[] isn’t declared constexpr in c++14 (though it looks like this will change).
non-constexpr function 'operator[]' cannot be used in a constant expression
for (size_t i = 0; i < Mbrot::pix; ++i) a[i] = Mbrot::IsInSet(Mbrot::Pos(i));
To get ourselves out of this mess with our dignity, we’ll need to introduce a new c++14 class template representing a compile-time sequence of integers:
template< class T, T... Ints >
class integer_sequence;
As seen in the definition, an std::integer_sequence is templated over the integer type, and the integers in the sequence. For example:
using sequence = std::integer_sequence<size_t, 0, 5, 7>;
Additionally, the standard library provides a few useful helper functions, for the cases where the integer type is size_t and where the sequence goes from .
using TSeq = std::integer_sequence<size_t, 0, 1, 2, 3>;
static_assert(std::is_same<TSeq, std::index_sequence<0, 1, 2, 3>>::value, "");
static_assert(std::is_same<TSeq, std::make_integer_sequence<size_t, 4>>::value, "");
static_assert(std::is_same<TSeq, std::make_index_sequence<4>>::value, "");
With that background, let’s get back to our Mandelbrot code:
struct Mbrot
{
// ...
// Type of the array to be returned.
using TArray = std::array<bool, Mbrot::pix>;
// std::integer_sequence instantiation, necessary.
using TSequence = std::make_index_sequence<Mbrot::pix>;
};
Using the c++11 syntax, we’re defining a boolean array type with length equal to the number of pixels and type. Easy enough. Next, we’re defining an integer sequence that is as if we had manually typed:
using TSequence = std::index_sequence<0, 1, 2, ..., Mbrot::pix>;
Here, the template parameter size_t... ints is 0, 1, 2, ..., Mbrot::pix. We exploit this fact when implementing the function that instantiates the array:
struct Mbrot
{
// ...
// Return the array.
template<size_t ...i>
static constexpr TArray
get_array_impl(std::index_sequence<i...> s)
{
if (!std::is_same<std::index_sequence<i...>, TSequence>::value)
{ throw std::out_of_range("Incorrect sequence passed."); }
return TArray{ { Mbrot::IsInSet(Mbrot::Pos(i))... } };
}
};
Let’s unpack that a bit.
This member function is templated over a parameter pack of integers, which it deduces from its argument, s.
To ensure that we haven’t passed the wrong integer sequence, we define a new index sequence type over the deduced parameter pack, and throw an exception if its type doesn’t match TSequence.
We then apply our member functions over the initializer list, creating a list of boolean values indicating membership in the set.
Note the position of the ellipse operator outside the function.
This causes a function to be applied to each member of a parameter pack, rather than passing the entire parameter pack to a single function call, as would happen if the ellipse were within the parentheses.
Finally, we use aggregate initialization to initialize the array.
Party Like It’s 1978
We don’t want our users to need to know what an std::integer_sequence is in order to use our class, so let’s define a nicer public interface:
struct Mbrot
{
// ...
// Return the array.
static constexpr TArray
GetArray()
{
constexpr auto s = Mbrot::TSequence();
return Mbrot::get_array_impl(s);
}
};
And let’s give our class a whirl!
int main()
{
constexpr auto a = Mbrot::GetArray();
for (unsigned int i = 0; i < Mbrot::pix; ++i)
{
std::cout << (a.at(i) ? " *" : " ");
if (Mbrot::str - 1 == i % Mbrot::str) std::cout << std::endl;
}
return EXIT_SUCCESS;
}
Finally, we can sit back and enjoy some hard earned, mid century modern graphics.
:-)
Almost as good as the original, huh?
As always, check out the repo for source.
* * *
* *
* * * *
* * * * *
* * * * *
* * * * * * * * * * *
* * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * *
* * * * * * * * * * *
* * * * *
* * * * *
* * * *
* *
* * *