Using IEEE 754 double precision values in the Mandelbrot Set Program developed in a previous post, we reach the limits of precision at about radius \(10^{-14}\). Beyond that images become blocky because we can't differentiate between near by pixels because they are closer than \(1 \cdot ULP(x)\).
We can use higher-precision floating point types but usually at a cost of performance. x86-64 based CPUs support 80-bit extended floating point, with about 19 decimal digits of accuracy, in hardware. Some platforms support 128-bit quad precision floating point, but that is usually a software implementation. MacOS on Apple Silicon ARM64 CPUs do not have any hardware support for extended floating point. Using the Gnu MP multi-precision library we can implement software arbitrary precision floating point for any architecture. The Gnu MP library is a low-level library, but built on top of the Gnu MP library is a high-level C++ wrapper called mpreal.h that we will use to make programming easier and cleaner.
By adding a few defines and typedefs, we can refactor our Mandelbrot Set program to use different floating point data types including the Gnu MP-based MPREAL.
#define USE_EXTENDED_FP // Define to use extended precision floating point
#define FORCE_MPREAL // Force use of GMP MPREALs
#if MAX_PIXEL_LEVELS > 256
#define PIXEL_TYPE uint16_t
#else
#define PIXEL_TYPE uint8_t
#endif
int DIGITS = 16; // Decimal digits of accuracy
#ifdef USE_EXTENDED_FP
#if defined __FLOAT128__ && !defined FORCE_MPREAL
typedef __float128 double_type;
#elif defined(__SIZEOF_LONG_DOUBLE__) && (__SIZEOF_LONG_DOUBLE__ > 8) && !defined FORCE_MPREAL
#define LONG_DOUBLE
typedef long double double_type;
#else
#define HAS_MPREAL
#include "mpreal.h" // MPFR C++ Wrapper
using namespace mpfr;
#define PRECISION int(ceil(DIGITS * 3.3)) // Precision in bits (~3.3 bits per decimal digit)
#define DOUBLE_DIGITS DIGITS
typedef mpreal double_type;
inline mpreal ULP(mpreal x) {
mpreal next = nextabove(x); // Smallest representable number greater than x
mpreal ulp = next - x;
return ulp;
}
#endif
#else
typedef double double_type;
#endif
#ifndef HAS_MPREAL
#define DOUBLE_DIGITS numeric_limits<double_type>::digits10
#define PRECISION numeric_limits<double_type>::digits
/**
* @brief Computes the Unit in the Last Place (ULP) of a double-precision value.
*
* This function calculates the distance between the given value @p x and the next
* representable double greater than @p x in the direction of positive infinity.
* The result gives a measure of the precision at the value @p x.
*
* @param x The input double value.
* @return The ULP of @p x (i.e., the difference to the next larger representable double).
*
* @note If @p x is NaN or infinity, the result may be undefined.
*/
inline double_type ULP(double_type x) {
double_type next = nextafter(x, numeric_limits<double_type>::infinity());
double_type ulp = next - x;
return ulp;
}
#endif
For the MPREAL type we add a new command line option --decimal_digits that we use to specify the precision that will be used by the MPREALs in our program. A good rule of thumb is to use about five more digits than the absolute value of the radius exponent. So if we want to zoom in to a radius of \(10^{-20}\) we should use about 25 decimal digits of accuracy.
After refactoring our Mandelbrot Set program to use the new double_type typedef instead of double we can compile our mandelbrot_threaded_color_mp.cpp program to use more precise datatypes if we define use_extended_fp.
Complete mandelbrot_threaded_color_mp.cpp
/**
* @file mandelbrot_threaded_color_mp.cpp
* @brief Program to compute the Mandelbrot Set.
*
* This program can compute images of the Mandelbrot Set in its
* entirety or close ups of boundary areas. It will use as many
* CPUs as are available to speed up the rendering. This program
* will use extended precision floating point. The Gnu Multi-Precision
* library is required if hardware extended precision FP is not
* available.
*
* It also allows two RGB colors to be specified to define the
* coloring of the set based on the escape value of a pixel.
* Colors can be specified using #xxx or #xxxxxx hexadecimal.
*
* Requires:
* Linux<br>
* get cli11.hpp at https://github.com/CLIUtils/CLI11/releases
* get mpreal.h at https://github.com/advanpix/mpreal.git
* sudo apt install libgmp-dev<br>
* sudo apt install libmpfr-dev<br>
*
* macOS<br>
* brew install cli11<br>
* brew install gmp<br>
* brew install mpfr<br>
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Build:
* Linux<br>
* put CLI.hpp in /usr/local/include <br>
* put mpreal.h in /usr/local/include <br>
*
* macOS<br>
* put mpreal.h in /opt/homebrew/include
* Add to .zshrc<br>
* export CPPFLAGS="-I/opt/homebrew/include"<br>
* export LDFLAGS="-L/opt/homebrew/lib"
*
* Compile with<br>
* Linux<br>
* If 128-bit doubles are available on your platform <br>
* g++ --std=gnu++20 -o mandelbrot_threaded_color_mp mandelbrot_threaded_color_mp.cpp -lquadmath
* to use GMP<br>
* g++ --std=gnu++20 -o mandelbrot_threaded mandelbrot_threaded_color.cpp -lmpfr -lgmp
* macOS<br>
* g++ $CPPFLAGS $LDFLAGS --std=gnu++20 -o mandelbrot_threaded_color_mp mandelbrot_threaded_color_mp.cpp -lmpfr -lgmp
*
* Execute:
* mandelbrot_threaded real imaginary radius [--filename=basename]
* [--image_size=pixels] [--max_iterations=count]
[--low_color=#xxx] [--high_color=#xxx] [--rainbow]
[--decimal_digits=n]
* Test multi-precision with command <br>
* ./mandelbrot_threaded_color_mp 0.3663629834227643413289327308 0.59153377326144522757749349375 1e-27 -m 100000 -d 30
*
* @author Richard Lesh
* @date 2025-03-12
* @version 1.0
*/
#include <CLI/CLI.hpp>
#include <complex>
#include <fstream>
#include <string>
#include <thread>
#include <vector>
#include "Color.hpp"
#define MAX_PIXEL_LEVELS 256 // 256 for 8-bit
bool USE_COLOR = true; // True to allow colors to be used
bool USE_RAINBOW = false; // True to use rainbow colors
#define USE_EXTENDED_FP // Define to use extended precision floating point
#define FORCE_MPREAL // Force use of GMP MPREALs
#if MAX_PIXEL_LEVELS > 256
#define PIXEL_TYPE uint16_t
#else
#define PIXEL_TYPE uint8_t
#endif
int DIGITS = 16; // Decimal digits of accuracy
#ifdef USE_EXTENDED_FP
#if defined __FLOAT128__ && !defined FORCE_MPREAL
typedef __float128 double_type;
#elif defined(__SIZEOF_LONG_DOUBLE__) && (__SIZEOF_LONG_DOUBLE__ > 8) && !defined FORCE_MPREAL
#define LONG_DOUBLE
typedef long double double_type;
#else
#define HAS_MPREAL
#include "mpreal.h" // MPFR C++ Wrapper
using namespace mpfr;
#define PRECISION int(ceil(DIGITS * 3.3)) // Precision in bits (~3.3 bits per decimal digit)
#define DOUBLE_DIGITS DIGITS
typedef mpreal double_type;
inline mpreal ULP(mpreal x) {
mpreal next = nextabove(x); // Smallest representable number greater than x
mpreal ulp = next - x;
return ulp;
}
#endif
#else
typedef double double_type;
#endif
#ifndef HAS_MPREAL
#define DOUBLE_DIGITS numeric_limits<double_type>::digits10
#define PRECISION numeric_limits<double_type>::digits
/**
* @brief Computes the Unit in the Last Place (ULP) of a double-precision value.
*
* This function calculates the distance between the given value @p x and the next
* representable double greater than @p x in the direction of positive infinity.
* The result gives a measure of the precision at the value @p x.
*
* @param x The input double value.
* @return The ULP of @p x (i.e., the difference to the next larger representable double).
*
* @note If @p x is NaN or infinity, the result may be undefined.
*/
inline double_type ULP(double_type x) {
double_type next = nextafter(x, numeric_limits<double_type>::infinity());
double_type ulp = next - x;
return ulp;
}
#endif
using namespace std;
typedef Color<PIXEL_TYPE> RGB;
// Array used to store the image in a thread safe way as an array of arrays
// Each element of the array will be a pointer to an array containing the
// pixels for a row of the image. Rows will be computed by threads.
unique_ptr<unique_ptr<int32_t[]>[]> image = nullptr;
// Atomic int used to assign tasks (row to compute) to the threads
atomic_int counter{0};
/**
* @brief Computes the squared magnitude (norm) of a complex number without
* taking the square root.
*
* This function returns the sum of the squares of the real and imaginary parts
* of the complex number @p z. It avoids the overhead of computing a square
* root, making it useful for performance-sensitive contexts like Mandelbrot
* set calculations where only relative magnitude comparisons are needed.
*
* @tparam T The numeric type (e.g., float, double, mpreal).
* @param z The complex number whose squared magnitude is to be computed.
* @return The squared magnitude of @p z, equivalent to |z|^2.
*
* @see norm may be used as a standard alternative, but may perform
* additional checks.
*/
template<typename T>
inline T norm_squared(complex<T> z) {
return z.real() * z.real() + z.imag() * z.imag();
}
template<typename T>
inline T magnitude(const complex<T>& z) {
return sqrt(z.real() * z.real() + z.imag() * z.imag());
}
// Number of pixels that escape
atomic_int num_escaped{0};
// Number of cycles detected
atomic_int num_cycles{0};
// Number of pixels that reach max. iterations
atomic_int num_max_iter{0};
/**
* @brief Compute the number of iterations to escape.
*
* Function computes the number of iterations of the formula
* z(n+1) = z(n)^2 + c that can be executed before the value of |z|
* exceeds the escape limit radius of 2.0.
*
* The function uses Floyd’s Cycle Detection (Tortoise
* and Hare) algorithm to detect cycles.
* @param c Complex point used to perform the calculation.
* @param max_iteration Number of iterations, if exceeded, when we will declare
* as not escaping and therefore in the Mandelbrot set.
* @return Value [0, max_iteration) representing the escape iterations.
* 0 is considered as not escaping.
*/
int compute_escape(complex<double_type> c, int max_iteration) {
complex<double_type> tortoise{0., 0.};
complex<double_type> hare{0., 0.};
int i = 0;
while (i < max_iteration && norm_squared(tortoise) < 4.0) {
tortoise = tortoise * tortoise + c;
hare = hare * hare + c;
hare = hare * hare + c; // hare moves twice
double_type ulps = ULP(magnitude(tortoise));
ulps = 9 * ulps * ulps;
if (norm_squared(tortoise - hare) < ulps) {
++num_cycles;
return 0L; // cycle detected
}
++i;
}
if (i == max_iteration) { // Return 0 (black) if we are in the set
++num_max_iter;
return 0L;
} else {
++num_escaped;
return i;
}
}
/**
* @brief Compute a portion of the Mandelbrot Set.
*
* Function computes a portion of the Mandelbrot Set row by row.
* The atomic integer counter is used to get the row to compute.
* The row is stored in the image array.
*
* @param x_min Minimum real axis bound to plot.
* @param x_max Maximum real axis bound to plot.
* @param y_min Minimum imaginary axis bound to plot.
* @param y_max Maximum imaginary axis bound to plot.
* @param max_iteration Limit to the number of iterations we will compute per point.
* @param image_size Size, in pixels, of the image to generate.
*/
void compute_mandelbrot(double_type x_min, double_type x_max, double_type y_min, double_type y_max, int max_iteration, int image_size) {
#ifdef HAS_MPREAL
mpreal::set_default_prec(PRECISION); // Set default precision in bits for each thread
#endif
double_type width_per_pixel = (x_max - x_min) / image_size;
double_type height_per_pixel = (y_max - y_min) / image_size;
int modulus = ceil(image_size / 40.); // Used to control progress output
while (true) {
int j = counter.fetch_add(1);
if (j < image_size) {
image[j] = make_unique<int32_t[]>(image_size);
// Calculate the imaginary coordinate of the middle of the pixel row
double_type y0 = y_max - height_per_pixel * (j + 0.5);
for (int i = 0; i < image_size; ++i) {
// Calculate the real coordinate of the middle of the pixel
double_type x0 = width_per_pixel * (i + 0.5) + x_min;
// Compute the escape number of iterations [0, max_iteration)
int level = compute_escape(complex(x0, y0), max_iteration);
// Store the escape iterations in the image array.
image[j][i] = level;
}
// Output periodic progress
if (j % modulus == 0) cout << '.' << flush;
} else {
break;
}
}
}
/**
* @brief Write out a PGM/PPM file.
*
* Function writes out image data to a PGM or PPM format file. Can be an
* 8-bit/16-bit PGM or an 8-bit/16-bit PPM. This function normalizes the
* level values using a histogram.
*
* @param filename Pathspec to the file that will be created.
* @param max_iteration Limit to the number of iterations we used to compute each point.
* @param image_size Size, in pixels, of the image to save.
* @param color1 First color to use for near 1 iteration before escaping.
* @param color2 Second color to use for near max_iterations iterations before escaping.
* @param use_rainbow true to use rainbow coloring.
*/
void output_file(const string &filename, int max_iteration, int image_size, RGB color1, RGB color2, bool use_rainbow) {
string ext = ".pgm";
if (USE_COLOR) ext = ".ppm";
cout << endl << "Creating " << (filename + ext) << "..." << endl;
ofstream out_file((filename + ext).c_str(), ios::binary);
if (!out_file) { // Check if file opened successfully
cerr << "Error opening file " << (filename + ext) << " for writing!" << endl;
exit(1);
}
if (USE_COLOR) {
// File header is P6 for binary pix map, x_size y_size, max pixel value per RGB component
out_file << "P6\n" << image_size << " " << image_size << "\n" << (MAX_PIXEL_LEVELS - 1) << "\n";
} else {
// File header is P5 for binary gray map, x_size y_size, max pixel value
out_file << "P5\n" << image_size << " " << image_size << "\n" << (MAX_PIXEL_LEVELS - 1) << "\n";
}
// Find the maximum iteration count
int max = 0;
for (int j = 0; j < image_size; ++j) {
for (int i = 0; i < image_size; ++i) {
if (image[j][i] > max) max = image[j][i];
}
}
max_iteration = max + 1;
// Compute bucket size for histogram
int max_iteration_bits = static_cast<int>(ceil(log(max_iteration)/log(2.0)));
int bucket_shift_bits = max_iteration_bits > 20 ? max_iteration_bits - 20 : 0;
vector<int>::size_type histogram_size = (max_iteration >> bucket_shift_bits) + 1;
// Compute histogram
vector<int> histogram(histogram_size, 0);
int totalCount = 0;
for (int j = 0; j < image_size; ++j) {
for (int i = 0; i < image_size; ++i) {
if (image[j][i] > 0) {
++histogram[vector<int>::size_type(image[j][i] >> bucket_shift_bits)];
++totalCount;
}
}
}
// Compute cumulative histogram
int count = 0;
for (vector<int>::size_type i = 1; i < histogram_size; ++i) {
count += histogram[i];
histogram[i] = count;
}
// Normalize using histogram
for (int j = 0; j < image_size; ++j) {
for (int i = 0; i < image_size; ++i) {
if (USE_COLOR) {
RGB rgb(0,0,0);
if (image[j][i] != 0) {
if (use_rainbow)
rgb = level_2_RGB<PIXEL_TYPE>(
histogram[vector<int>::size_type(image[j][i] >> bucket_shift_bits)] /
static_cast<double>(totalCount + 1));
else
rgb = blend_colors(
histogram[vector<int>::size_type(image[j][i] >> bucket_shift_bits)] /
static_cast<double>(totalCount + 1), color1, color2);
}
if constexpr (sizeof(PIXEL_TYPE) == 2) {
out_file.write(reinterpret_cast<const char*>(&rgb.R) + 1, 1);
out_file.write(reinterpret_cast<const char*>(&rgb.R), 1);
out_file.write(reinterpret_cast<const char*>(&rgb.G) + 1, 1);
out_file.write(reinterpret_cast<const char*>(&rgb.G), 1);
out_file.write(reinterpret_cast<const char*>(&rgb.B) + 1, 1);
out_file.write(reinterpret_cast<const char*>(&rgb.B), 1);
} else {
out_file.write(reinterpret_cast<const char*>(&rgb.R), 1);
out_file.write(reinterpret_cast<const char*>(&rgb.G), 1);
out_file.write(reinterpret_cast<const char*>(&rgb.B), 1);
}
} else {
if constexpr (sizeof(PIXEL_TYPE) == 2) {
PIXEL_TYPE level = 0;
if (image[j][i] != 0)
level = static_cast<uint8_t>(color1.R +
histogram[vector<int>::size_type(image[j][i] >> bucket_shift_bits)] /
static_cast<double>(totalCount + 1) * (color2.R - color1.R + 1));
out_file.write(reinterpret_cast<const char*>(&level) + 1, 1);
out_file.write(reinterpret_cast<const char*>(&level), 1);
} else {
PIXEL_TYPE level = 0;
if (image[j][i] != 0)
level = static_cast<uint8_t>(color1.R +
histogram[vector<int>::size_type(image[j][i] >> bucket_shift_bits)] /
static_cast<double>(totalCount + 1) * (color2.R - color1.R + 1));
out_file.write(reinterpret_cast<const char*>(&level), 1);
}
}
}
}
out_file.close();
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: mandelbrot_threaded_color_mp real imaginary radius [--filename=basename]
* [--image_size=n] [--max_iterations=n] [--high_color=#xxx] [--low_color=#xxx]
* [--rainbow] [--decimal_digits=n]
*
* real & imaginary specify the coordinates of the Mandelbrot Set to display at
* the center of the image.
*
* radius is the radius in Mandelbrot Set coordinates out to the extreme edges
* of the image.
*
* max_iterations is the maximum number of iterations per pixel allowed to
* determine if it is in the set. The max_iterations will generally need to get
* larger as radius gets smaller.
*
* image_size is the size in pixels for the width and height of the square image.
*
* high_color and low_color are hexadecimal color values.
*
* filename is the basename of the image file. It defaults to "image".
*
* decimal_digits sets the accuracy of floating point computations when using GMP.
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char** argv) {
CLI::App app{"Program to generate Mandelbrot Set images."};
app.set_help_flag("", "");
argv = app.ensure_utf8(argv);
// use CLI11 to process command line options.
vector<string> coords;
app.add_option("coords", coords, "Location to plot (x y radius)")
->expected(3);
// Optional arguments
string filename = "image";
int max_iteration = 10000;
int image_size = 512;
string high_color = "#FFF";
string low_color = "#000";
bool use_rainbow = false;
int decimal_digits = 16;
bool show_help = false;
app.add_flag("--help", show_help, "Show help");
app.add_option("-f,--filename", filename, "Base filename (no extension)");
app.add_option("-m,--max_iterations", max_iteration, "Max. iterations");
app.add_option("-s,--image_size", image_size, "Image size in pixels");
app.add_option("-h,--high_color", high_color, "High color hex code (#xxx or #xxxxxx)");
app.add_option("-l,--low_color", low_color, "Low color hex code (#xxx or #xxxxxx)");
app.add_flag("-r,--rainbow", use_rainbow, "Use rainbow colors");
app.add_option("-d,--decimal_digits", decimal_digits, "Decimal digits of accuracy to use");
CLI11_PARSE(app, argc, argv);
if (show_help) {
cout << app.help() << endl;
return 0;
}
// Get the three required coordinates to plot
double_type x, y, radius;
if (coords.size() == 3) {
#ifdef HAS_MPREAL
DIGITS = decimal_digits;
mpreal::set_default_prec(PRECISION); // Set default precision in bits
x = double_type(coords[vector<string>::size_type(0)]);
y = double_type(coords[vector<string>::size_type(1)]);
radius = double_type(coords[vector<string>::size_type(2)]);
#elif defined __FLOAT128__
x = strtoflt128(coords[0], nullptr);
y = strtoflt128(coords[1], nullptr);
radius = strtoflt128(coords[2], nullptr);
#elif defined LONG_DOUBLE
x = stold(coords[0]);
y = stold(coords[1]);
radius = stold(coords[2]);
#else
x = stod(coords[0]);
y = stod(coords[1]);
radius = stod(coords[2]);
#endif
} else {
cout << "Arguments x, y and radius required but not found!" << endl;
return 1;
}
// Check the image size
if (image_size > sqrt(numeric_limits<int>::max())) {
cout << "Image size can not exceed " << floor(sqrt(numeric_limits<int>::max())) << endl;
return 1;
}
// Check the max number of iterations
if (max_iteration > (numeric_limits<int32_t>::max() / 2)) {
cout << "Max iterations can not exceed " << numeric_limits<int32_t>::max() / 2 << endl;
return 1;
}
// Get the low and high colors
RGB high_RGB = RGB(high_color);
RGB low_RGB = RGB(low_color);
if (high_RGB.isGray() && low_RGB.isGray()) USE_COLOR = false;
if (use_rainbow) USE_COLOR = true;
// Compute limits of image
const double_type x_min = x - radius;
const double_type x_max = x + radius;
const double_type y_min = y - radius;
const double_type y_max = y + radius;
// Output settings
cout << "Size of floating point: " << sizeof(double_type) << endl;
cout << "Floating point decimal digits: " << DOUBLE_DIGITS << endl;
cout << "Floating point precision: " << PRECISION << " bits" << endl;
cout << "Maximum Iterations: " << max_iteration << endl;
cout << (USE_COLOR ? "Using color" : "Using grayscale") << endl;
if (use_rainbow) {
cout << "Rainbow colors" << endl;
} else {
cout << "High escape color: " << high_color << endl;
cout << "Low escape color: " << low_color << endl;
}
cout << "Image Bounds: " << setprecision(DOUBLE_DIGITS) << x_min << ", " << x_max << ", " << y_min << ", " << y_max << endl;
cout << "Image Size: " << image_size << " x " << image_size << endl;
// Allocate an array for the row pointers.
image = make_unique<unique_ptr<int32_t[]>[]>(image_size);
// Do a simple check to see if our doubles have enough precision
if ((x_max - x_min) / image_size < ULP(x) ||
(y_max - y_min) / image_size < ULP(y))
cout << "Doubles may not have enough precision!" << endl;
// Get the number of CPU cores.
unsigned int num_cpus = thread::hardware_concurrency();
if (num_cpus == 0) {
cout << "Could not determine number of CPUs. Using 4.\n";
num_cpus = 4;
} else {
cout << "Number of CPU cores: " << num_cpus << endl;
}
// Create the threads
vector<thread> threads;
threads.reserve(num_cpus);
for (int i = 0; i < num_cpus; ++i) {
threads.emplace_back(compute_mandelbrot, x_min, x_max, y_min, y_max,
max_iteration, image_size);
}
// Wait for all threads to complete
for (auto& t : threads) {
t.join();
}
cout << endl;
cout << "Number of escaped pixels: " << num_escaped << endl;
cout << "Number of cycles detected: " << num_cycles << endl;
cout << "Number of max iteration pixels: " << num_max_iter << endl;
// Generate the output file
output_file(filename, max_iteration, image_size, low_RGB, high_RGB, use_rainbow);
return 0;
}
Perturbation Theory
The big disadvantage to using software multi-precision is that it is much slower than hardware floating point. Even with cycle detection, it can take a long time to compute images with high zoom factors. As we zoom in, we also typically need to set the maximum iterations higher which compounds the performance problem.
We can take a page from perturbation theory and develop an algorithm that is faster by using multi-precision to compute the orbit for just one pixel. Then because other pixels nearby are very close at high magnification factors, we can use perturbations to compute just the differences from the reference orbit. These differences can be computed using hardware IEEE754 floats so these computations can be done very quickly compared to multi-precision computations. Zhuoran posted just such an algorithm on the Fractal Forums in 2021. The program mandelbrot_threads_color_mp_fast.cpp implements the Zhuoran algorithm with a mixture of MPREALs and doubles to achieve much faster rendering at large magnifications.
Pseudocode Algorithm by Zhuoran
complex Reference[]; // Reference orbit (MUST START WITH ZERO)
int MaxRefIteration; // The last valid iteration of the reference (the iteration just before it escapes)
complex dz = 0, dc; // Delta z and Delta c
int Iteration = 0, RefIteration = 0;
while (Iteration < MaxIteration) {
dz = 2 * dz * Reference[RefIteration] + dz * dz + dc;
RefIteration++;
complex z = Reference[RefIteration] + dz;
if (|z| > BailoutRadius) break;
if (|z| < |dz| || RefIteration == MaxRefIteration) {
dz = z;
RefIteration = 0;
}
Iteration++;
}
Complete mandelbrot_threads_color_mp_fast.cpp
/**
* @file mandelbrot_threaded_color_mp_fast.cpp
* @brief Program to compute the Mandelbrot Set.
*
* This program can compute images of the Mandelbrot Set in its
* entirety or close ups of boundary areas. It will use as many
* CPUs as are available to speed up the rendering. This program
* will use extended precision floating point. The Gnu Multi-Precision
* library is required if hardware extended precision FP is not
* available. The algorithm used is based on perturbation theory
* so that most computations can be done with doubles.
*
* It also allows two RGB colors to be specified to define the
* coloring of the set based on the escape value of a pixel.
* Colors can be specified using #xxx or #xxxxxx hexadecimal.
*
* Requires:
* Linux<br>
* get cli11.hpp at https://github.com/CLIUtils/CLI11/releases
* get mpreal.h at https://github.com/advanpix/mpreal.git
* sudo apt install libgmp-dev<br>
* sudo apt install libmpfr-dev<br>
*
* macOS<br>
* brew install cli11<br>
* brew install gmp<br>
* brew install mpfr<br>
* get mpreal.h at https://github.com/advanpix/mpreal.git
*
* Build:
* Linux<br>
* put CLI.hpp in /usr/local/include <br>
* put mpreal.h in /usr/local/include <br>
*
* macOS<br>
* put mpreal.h in /opt/homebrew/include
* Add to .zshrc<br>
* export CPPFLAGS="-I/opt/homebrew/include"<br>
* export LDFLAGS="-L/opt/homebrew/lib"
*
* Compile with<br>
* Linux<br>
* If 128-bit doubles are available on your platform <br>
* g++ --std=gnu++20 -o mandelbrot_threaded_color_mp_fast mandelbrot_threaded_color_mp_fast.cpp -lquadmath
* to use GMP<br>
* g++ --std=gnu++20 -o mandelbrot_threaded_color_mp_fast mandelbrot_threaded_color_mp_fast.cpp -lmpfr -lgmp
* macOS<br>
* g++ $CPPFLAGS $LDFLAGS --std=gnu++20 -o mandelbrot_threaded_color_mp_fast mandelbrot_threaded_color_mp_fast.cpp -lmpfr -lgmp
*
* Execute:
* mandelbrot_threaded_color_mp_fast real imaginary radius [--filename=basename]
* [--image_size=pixels] [--max_iterations=count]
* [--low_color=#xxx] [--high_color=#xxx] [--rainbow]
* [--decimal_digits=n] [--cycle_detection]
*
* Test multi-precision with command <br>
* ./mandelbrot_threaded_color_mp_fast 0.3663629834227643413289327308 0.59153377326144522757749349375 1e-27 -m 100000 -d 30
*
* @author Richard Lesh
* @date 2025-03-12
* @version 1.0
*/
#include <bit>
#include <CLI/CLI.hpp>
#include <cassert>
#include <complex>
#include <fstream>
#include <string>
#include <thread>
#include <unordered_set>
#include <vector>
#include "Color.hpp"
#define MAX_PIXEL_LEVELS 256 // 256 for 8-bit
bool USE_COLOR = true; // True to allow colors to be used
bool USE_RAINBOW = false; // True to use rainbow colors
#define USE_EXTENDED_FP // Define to use extended precision floating point
#define FORCE_MPREAL // Force use of GMP MPREALs
#if MAX_PIXEL_LEVELS > 256
#define PIXEL_TYPE uint16_t
#else
#define PIXEL_TYPE uint8_t
#endif
int DIGITS = 16; // Decimal digits of accuracy
#ifdef USE_EXTENDED_FP
#if defined __FLOAT128__ && !defined FORCE_MPREAL
typedef __float128 double_type;
#elif defined(__SIZEOF_LONG_DOUBLE__) && (__SIZEOF_LONG_DOUBLE__ > 8) && !defined FORCE_MPREAL
#define LONG_DOUBLE
typedef long double double_type;
#else
#define HAS_MPREAL
#include "mpreal.h" // MPFR C++ Wrapper
using namespace mpfr;
#define PRECISION int(ceil(DIGITS * 3.3)) // Precision in bits (~3.3 bits per decimal digit)
#define DOUBLE_DIGITS DIGITS
typedef mpreal double_type;
inline mpreal ULP(mpreal x) {
mpreal next = nextabove(x); // Smallest representable number greater than x
mpreal ulp = next - x;
return ulp;
}
#endif
#else
typedef double double_type;
#endif
#ifndef HAS_MPREAL
#define DOUBLE_DIGITS numeric_limits<double_type>::digits10
#define PRECISION numeric_limits<double_type>::digits
/**
* @brief Computes the Unit in the Last Place (ULP) of a double-precision value.
*
* This function calculates the distance between the given value @p x and the next
* representable double greater than @p x in the direction of positive infinity.
* The result gives a measure of the precision at the value @p x.
*
* @param x The input double value.
* @return The ULP of @p x (i.e., the difference to the next larger representable double).
*
* @note If @p x is NaN or infinity, the result may be undefined.
*/
inline double_type ULP(double_type x) {
double_type next = nextafter(x, numeric_limits<double_type>::infinity());
double_type ulp = next - x;
return ulp;
}
#endif
using namespace std;
typedef Color<PIXEL_TYPE> RGB;
constexpr std::string::size_type operator"" _sz(unsigned long long n) {
return static_cast<size_t>(n);
}
// Array used to store the image in a thread safe way as an array of arrays
// Each element of the array will be a pointer to an array containing the
// pixels for a row of the image. Rows will be computed by threads.
unique_ptr<unique_ptr<int32_t[]>[]> image = nullptr;
// Atomic int used to assign tasks (row to compute) to the threads
atomic_int counter{0};
/**
* @brief Computes the squared magnitude (norm) of a complex number without
* taking the square root.
*
* This function returns the sum of the squares of the real and imaginary parts
* of the complex number @p z. It avoids the overhead of computing a square
* root, making it useful for performance-sensitive contexts like Mandelbrot
* set calculations where only relative magnitude comparisons are needed.
*
* @tparam T The numeric type (e.g., float, double, mpreal).
* @param z The complex number whose squared magnitude is to be computed.
* @return The squared magnitude of @p z, equivalent to |z|^2.
*
* @see norm may be used as a standard alternative, but may perform
* additional checks.
*/
template<typename T>
inline T norm_squared(complex<T> z) {
return z.real() * z.real() + z.imag() * z.imag();
}
// Number of pixels that escape
atomic_int num_escaped{0};
// Number of cycles detected
atomic_int num_cycles{0};
// Number of pixels that reach max. iterations
atomic_int num_max_iter{0};
/**
* @brief Compute the number of iterations to escape.
*
* Function computes the number of iterations of the formula
* z(n+1) = z(n)^2 + c that can be executed before the value of |z|
* exceeds the escape limit radius of 2.0.
*
* @param c Complex point used to perform the calculation.
* @param max_iteration Number of iterations, if exceeded, when we will declare
* as not escaping and therefore in the Mandelbrot set.
* @return vector of complex points representing the orbit of c.
*/
vector<complex<double_type>>
compute_escape(complex<double_type> c, int max_iteration) {
int iteration = 0; // Number of iterations before escaping
complex<double_type> z{0., 0.};
vector<complex<double_type>> reference;
reference.push_back(z);
for (int i = 0; i < max_iteration && norm_squared(z) < 1e20; ++i) {
z = z * z + c;
reference.push_back(z);
++iteration;
}
return reference;
}
// Convert double to its bitwise representation as int64_t
inline int64_t toULP(double x) {
auto i = std::bit_cast<int64_t>(x);
// adjust negative ones so that the resulting int64_t sort the same as the doubles.
return i < 0 ? numeric_limits<int64_t>::min() - i : i;
}
struct QuantizedComplex {
int64_t ulp_real;
int64_t ulp_imag;
int64_t d_ulp_real;
int64_t d_ulp_imag;
QuantizedComplex(complex<double> x, complex<double> dx, int ulp_bin_shift = 2) {
// Quantize by grouping into bins of ulp_bin ULPs
ulp_real = toULP(x.real()) >> ulp_bin_shift;
ulp_imag = toULP(x.imag()) >> ulp_bin_shift;
d_ulp_real = toULP(dx.real()) >> ulp_bin_shift;
d_ulp_imag = toULP(dx.imag()) >> ulp_bin_shift;
}
bool operator==(const QuantizedComplex& other) const {
return ulp_real == other.ulp_real && ulp_imag == other.ulp_imag &&
d_ulp_real == other.d_ulp_real && d_ulp_imag == other.d_ulp_imag;
}
static bool ulp_diff(double a, double b, int max_ulps) {
return std::abs(toULP(a) - toULP(b)) <= max_ulps;
}
};
// Hash based on quantized ULPs
template <>
struct std::hash<QuantizedComplex> {
size_t operator()(const QuantizedComplex& q) const noexcept {
return q.ulp_real ^ q.ulp_imag ^ q.d_ulp_real ^ q.d_ulp_imag;
}
};
/**
* @brief Compute a portion of the Mandelbrot Set.
*
* Function computes a portion of the Mandelbrot Set row by row.
* The atomic integer counter is used to get the row to compute.
* The row is stored in the image array.
*
* @param cell_size size of the computation cell in pixels
* @param x_min Minimum real axis bound to plot.
* @param x_max Maximum real axis bound to plot.
* @param y_min Minimum imaginary axis bound to plot.
* @param y_max Maximum imaginary axis bound to plot.
* @param max_iteration Limit to the number of iterations we will compute per point.
* @param image_size Size, in pixels, of the image to generate.
*/
void compute_mandelbrot(int cell_size, const double_type& x_min, const double_type& x_max, const double_type& y_min, const double_type& y_max, int max_iteration, int image_size, bool use_cycle_detection = false) {
#ifdef HAS_MPREAL
mpreal::set_default_prec(PRECISION); // Set default precision in bits for each thread
#endif
int num_cells = image_size / cell_size;
double_type width = x_max - x_min;
double_type height = y_max - y_min;
int modulus = ceil(num_cells * num_cells / 40.);
int cycle_detection_limit = 10000;
double_type cell_width = (x_max - x_min) / num_cells;
double_type cell_height = (y_max - y_min) / num_cells;
while (true) {
int cell_num = counter.fetch_add(1);
if (cell_num < num_cells * num_cells) {
int i = cell_num % num_cells;
int j = cell_num / num_cells;
double_type cell_x_min = x_min + i * cell_width;
double_type cell_x_max = x_min + i * cell_width + cell_width;
double_type cell_y_min = y_max - j * cell_width - cell_width;
double_type cell_y_max = y_max - j * cell_width;
complex<double_type> center{(cell_x_max + cell_x_min) / 2.,
(cell_y_max + cell_y_min) / 2.};
vector<complex<double_type>> reference = compute_escape(center, max_iteration);
vector<double_type>::size_type max_ref_iteration = reference.size() - 1;
assert(max_ref_iteration > 1);
bool enable_cycle_detection = false;
for (int i0 = 0; i0 < cell_size; ++i0) {
for (int j0 = 0; j0 < cell_size; ++j0) {
double_type y0 = cell_y_max - (cell_height / cell_size) * (j0 + 0.5);
double_type x0 = (cell_width / cell_size) * (i0 + 0.5) + cell_x_min;
complex<double_type> point(x0, y0);
complex<double_type> dc_full = point - center;
complex<double> dc(static_cast<double>(dc_full.real()), static_cast<double>(dc_full.imag()));
int iteration = 0;
vector<double_type>::size_type ref_iteration = 0;
complex<double> dz(0., 0.);
complex<double_type> ref = reference[ref_iteration];
complex<double> ref_z(
static_cast<double>(ref.real()),
static_cast<double>(ref.imag()));
std::unordered_set<QuantizedComplex> seen;
bool cycle_found = false;
while (iteration < max_iteration) {
dz = 2. * dz * ref_z + dz * dz + dc;
++ref_iteration;
ref = reference[ref_iteration];
ref_z = complex<double>(
static_cast<double>(ref.real()),
static_cast<double>(ref.imag()));
complex<double> z = ref_z + dz;
double norm = norm_squared(z);
if (norm > 4.0) break;
if (enable_cycle_detection && iteration > cycle_detection_limit) {
QuantizedComplex qz(ref_z, dz);
if (seen.count(qz)) {
cycle_found = true;
break;
} else {
seen.insert(qz);
}
}
if (norm < norm_squared(dz) || ref_iteration == max_ref_iteration) {
dz = z;
ref_iteration = 0;
ref = reference[ref_iteration];
ref_z = complex<double>(
static_cast<double>(ref.real()),
static_cast<double>(ref.imag()));
}
++iteration;
}
if (cycle_found) {
++num_cycles;
cycle_detection_limit = static_cast<int>(0.9 * iteration);
iteration = 0;
} else if (iteration == max_iteration) {
++num_max_iter;
enable_cycle_detection = use_cycle_detection;
cycle_detection_limit = static_cast<int>(0.9 * max_iteration);
} else {
++num_escaped;
if (cycle_detection_limit < iteration) {
cycle_detection_limit = static_cast<int>(1.1 * iteration);
}
}
image[j * cell_size + j0][i * cell_size + i0] = iteration == max_iteration ? 0 : iteration;
}
}
if (cell_num % modulus == 0) cout << '.' << flush;
} else {
break;
}
}
}
/**
* @brief Write out a PGM/PPM file.
*
* Function writes out image data to a PGM or PPM format file. Can be an
* 8-bit/16-bit PGM or an 8-bit/16-bit PPM. This function normalizes the
* level values using a histogram.
*
* @param filename Pathspec to the file that will be created.
* @param max_iteration Limit to the number of iterations we used to compute each point.
* @param image_size Size, in pixels, of the image to save.
* @param color1 First color to use for near 1 iteration before escaping.
* @param color2 Second color to use for near max_iterations iterations before escaping.
* @param use_rainbow true to use rainbow coloring.
*/
void output_file(const string &filename, int max_iteration, int image_size, RGB color1, RGB color2, bool use_rainbow) {
string ext = ".pgm";
if (USE_COLOR) ext = ".ppm";
cout << endl << "Creating " << (filename + ext) << "..." << endl;
ofstream out_file((filename + ext).c_str(), ios::binary);
if (!out_file) { // Check if file opened successfully
cerr << "Error opening file " << (filename + ext) << " for writing!" << endl;
exit(1);
}
if (USE_COLOR) {
// File header is P6 for binary pix map, x_size y_size, max pixel value per RGB component
out_file << "P6\n" << image_size << " " << image_size << "\n" << (MAX_PIXEL_LEVELS - 1) << "\n";
} else {
// File header is P5 for binary gray map, x_size y_size, max pixel value
out_file << "P5\n" << image_size << " " << image_size << "\n" << (MAX_PIXEL_LEVELS - 1) << "\n";
}
// Find the maximum iteration count
int max = 0;
for (int j = 0; j < image_size; ++j) {
for (int i = 0; i < image_size; ++i) {
if (image[j][i] > max) max = image[j][i];
}
}
max_iteration = max + 1;
// Compute bucket size for histogram
int max_iteration_bits = static_cast<int>(ceil(log(max_iteration)/log(2.0)));
int bucket_shift_bits = max_iteration_bits > 20 ? max_iteration_bits - 20 : 0;
vector<int>::size_type histogram_size = (max_iteration >> bucket_shift_bits) + 1;
// Compute histogram
vector<int> histogram(histogram_size, 0);
int totalCount = 0;
for (int j = 0; j < image_size; ++j) {
for (int i = 0; i < image_size; ++i) {
if (image[j][i] > 0) {
++histogram[vector<int>::size_type(image[j][i] >> bucket_shift_bits)];
++totalCount;
}
}
}
// Compute cumulative histogram
int count = 0;
for (vector<int>::size_type i = 1; i < histogram_size; ++i) {
count += histogram[i];
histogram[i] = count;
}
// Normalize using histogram
for (int j = 0; j < image_size; ++j) {
for (int i = 0; i < image_size; ++i) {
if (USE_COLOR) {
RGB rgb(0,0,0);
if (image[j][i] != 0) {
if (use_rainbow)
rgb = level_2_RGB<PIXEL_TYPE>(
histogram[vector<int>::size_type(image[j][i] >> bucket_shift_bits)] /
static_cast<double>(totalCount + 1));
else
rgb = blend_colors(
histogram[vector<int>::size_type(image[j][i] >> bucket_shift_bits)] /
static_cast<double>(totalCount + 1), color1, color2);
}
if constexpr (sizeof(PIXEL_TYPE) == 2) {
out_file.write(reinterpret_cast<const char*>(&rgb.R) + 1, 1);
out_file.write(reinterpret_cast<const char*>(&rgb.R), 1);
out_file.write(reinterpret_cast<const char*>(&rgb.G) + 1, 1);
out_file.write(reinterpret_cast<const char*>(&rgb.G), 1);
out_file.write(reinterpret_cast<const char*>(&rgb.B) + 1, 1);
out_file.write(reinterpret_cast<const char*>(&rgb.B), 1);
} else {
out_file.write(reinterpret_cast<const char*>(&rgb.R), 1);
out_file.write(reinterpret_cast<const char*>(&rgb.G), 1);
out_file.write(reinterpret_cast<const char*>(&rgb.B), 1);
}
} else {
if constexpr (sizeof(PIXEL_TYPE) == 2) {
PIXEL_TYPE level = 0;
if (image[j][i] != 0)
level = static_cast<uint8_t>(color1.R +
histogram[vector<int>::size_type(image[j][i] >> bucket_shift_bits)] /
static_cast<double>(totalCount + 1) * (color2.R - color1.R + 1));
out_file.write(reinterpret_cast<const char*>(&level) + 1, 1);
out_file.write(reinterpret_cast<const char*>(&level), 1);
} else {
PIXEL_TYPE level = 0;
if (image[j][i] != 0)
level = static_cast<uint8_t>(color1.R +
histogram[vector<int>::size_type(image[j][i] >> bucket_shift_bits)] /
static_cast<double>(totalCount + 1) * (color2.R - color1.R + 1));
out_file.write(reinterpret_cast<const char*>(&level), 1);
}
}
}
}
out_file.close();
}
/**
* @brief Main function.
*
* The entry point of the program.
*
* Syntax: mandelbrot_threaded_color_mp_fast real imaginary radius [--filename=basename]
* [--image_size=n] [--max_iterations=n] [--high_color=#xxx] [--low_color=#xxx]
* [--rainbow] [--decimal_digits=n] [--cycle_detection]
*
* real & imaginary specify the coordinates of the Mandelbrot Set to display at
* the center of the image.
*
* radius is the radius in Mandelbrot Set coordinates out to the extreme edges
* of the image.
*
* max_iterations is the maximum number of iterations per pixel allowed to
* determine if it is in the set. The max_iterations will generally need to get
* larger as radius gets smaller.
*
* image_size is the size in pixels for the width and height of the square image.
*
* high_color and low_color are hexadecimal color values.
*
* filename is the basename of the image file. It defaults to "image".
*
* decimal_digits sets the accuracy of floating point computations when using GMP.
*
* @return int Returns 0 on successful execution.
*/
int main(int argc, char** argv) {
CLI::App app{"Program to generate Mandelbrot Set images."};
app.set_help_flag("", "");
argv = app.ensure_utf8(argv);
// use CLI11 to process command line options.
vector<string> coords;
app.add_option("coords", coords, "Location to plot (x y radius)")
->expected(3);
// Optional arguments
string filename = "image";
int max_iteration = 10000;
int image_size = 512;
string high_color = "#FFF";
string low_color = "#000";
bool use_rainbow = false;
int decimal_digits = 16;
int cell_size = 0;
bool use_cycle_detection = false;
bool show_help = false;
app.add_flag("--help", show_help, "Show help");
app.add_option("-f,--filename", filename, "Base filename (no extension)");
app.add_option("-m,--max_iterations", max_iteration, "Max. iterations");
app.add_option("-s,--image_size", image_size, "Image size in pixels");
app.add_option("-h,--high_color", high_color, "High color hex code (#xxx or #xxxxxx)");
app.add_option("-l,--low_color", low_color, "Low color hex code (#xxx or #xxxxxx)");
app.add_flag("-r,--rainbow", use_rainbow, "Use rainbow colors");
app.add_option("-d,--decimal_digits", decimal_digits, "Decimal digits of accuracy to use");
app.add_option("-c,--cell_size", cell_size, "Size in pixels of computation cell");
app.add_flag("-C,--cycle_detection", use_cycle_detection, "Use cycle detection algorithm");
CLI11_PARSE(app, argc, argv);
if (show_help) {
cout << app.help() << endl;
return 0;
}
// Get the three required coordinates to plot
double_type x, y, radius;
if (coords.size() == 3) {
#ifdef HAS_MPREAL
DIGITS = decimal_digits;
mpreal::set_default_prec(PRECISION); // Set default precision in bits
x = double_type(coords[vector<string>::size_type(0)]);
y = double_type(coords[vector<string>::size_type(1)]);
radius = double_type(coords[vector<string>::size_type(2)]);
#elif defined __FLOAT128__
x = strtoflt128(coords[0], nullptr);
y = strtoflt128(coords[1], nullptr);
radius = strtoflt128(coords[2], nullptr);
#elif defined LONG_DOUBLE
x = stold(coords[0]);
y = stold(coords[1]);
radius = stold(coords[2]);
#else
x = stod(coords[0]);
y = stod(coords[1]);
radius = stod(coords[2]);
#endif
} else {
cout << "Arguments x, y and radius required but not found!" << endl;
return 1;
}
// Check the image size
if (image_size > sqrt(numeric_limits<int>::max())) {
cout << "Image size can not exceed " << floor(sqrt(numeric_limits<int>::max())) << endl;
return 1;
}
// Check the max number of iterations
if (max_iteration > (numeric_limits<int32_t>::max() / 2)) {
cout << "Max iterations can not exceed " << numeric_limits<int32_t>::max() / 2 << endl;
return 1;
}
// Get the low and high colors
RGB high_RGB = RGB(high_color);
RGB low_RGB = RGB(low_color);
if (high_RGB.isGray() && low_RGB.isGray()) USE_COLOR = false;
if (use_rainbow) USE_COLOR = true;
// Compute limits of image
const double_type x_min = x - radius;
const double_type x_max = x + radius;
const double_type y_min = y - radius;
const double_type y_max = y + radius;
// Get the number of CPU cores.
unsigned int num_cpus = thread::hardware_concurrency();
if (num_cpus == 0) {
cout << "Could not determine number of CPUs. Using 4.\n";
num_cpus = 4;
} else {
cout << "Number of CPU cores: " << num_cpus << endl;
}
if (cell_size <= 0) {
cell_size = static_cast<int>(image_size / std::ceil(std::sqrt(2 * num_cpus)));
if (cell_size <= 1) cell_size = 1;
else if (cell_size < 4) cell_size = 2;
else if (cell_size < 8) cell_size = 4;
else if (cell_size < 16) cell_size = 8;
else if (cell_size < 32) cell_size = 16;
else if (cell_size < 64) cell_size = 32;
else if (cell_size < 128) cell_size = 64;
else if (cell_size < 256) cell_size = 128;
else cell_size = 256;
}
image_size = ((image_size + cell_size - 1) / cell_size) * cell_size;
// Output settings
cout << "Size of floating point: " << sizeof(double_type) << endl;
cout << "Floating point decimal digits: " << DOUBLE_DIGITS << endl;
cout << "Floating point precision: " << PRECISION << " bits" << endl;
cout << "Maximum Iterations: " << max_iteration << endl;
cout << (USE_COLOR ? "Using color" : "Using grayscale") << endl;
if (use_rainbow) {
cout << "Rainbow colors" << endl;
} else {
cout << "High escape color: " << high_color << endl;
cout << "Low escape color: " << low_color << endl;
}
cout << "Image Bounds: " << setprecision(DOUBLE_DIGITS) << x_min << ", " << x_max << ", " << y_min << ", " << y_max << endl;
cout << "Image Size: " << image_size << " x " << image_size << endl;
cout << "Cell Size: " << cell_size << " x " << cell_size << endl;
cout << "Cycle Detection: " << (use_cycle_detection ? "enabled" : "disabled") << endl;
// Do a simple check to see if our doubles have enough precision
if ((x_max - x_min) / image_size < ULP(x) ||
(y_max - y_min) / image_size < ULP(y))
cout << "Doubles may not have enough precision!" << endl;
// Allocate an array for the row pointers.
image = make_unique<unique_ptr<int32_t[]>[]>(image_size);
for (int i = 0; i < image_size; ++i)
image[i] = make_unique<int32_t[]>(image_size);
vector<thread> threads;
threads.reserve(num_cpus);
// Create the threads
for (int i = 0; i < num_cpus; ++i) {
threads.emplace_back(compute_mandelbrot, cell_size, x_min, x_max, y_min,
y_max, max_iteration, image_size, use_cycle_detection);
}
// Wait for all threads to complete
for (auto& t : threads) {
t.join();
}
cout << endl;
cout << "Number of escaped pixels: " << num_escaped << endl;
if (use_cycle_detection)
cout << "Number of cycles detected: " << num_cycles << endl;
cout << "Number of max iteration pixels: " << num_max_iter << endl;
// Generate the output file
output_file(filename, max_iteration, image_size, low_RGB, high_RGB, use_rainbow);
return 0;
}
You can find an exploration of one point on the Mandelbrot Set at ever increasing magnifications that illustrates the use of the multi-precision version of the program at https://www.lesh.cloud/mandelbrot/deep_mandelbrot/