Mandelbrot Render Speed With Arbritary Precision vs Double Precision

Question

The Problem:

My C# code that renders the mandelbrot set does so in ~450ms, while my c++ code using MPFR and MPIR takes ~10,800ms at 5 iterations. Both of these times are fairly slow, but 10 seconds is glacial for rendering the set. I am aware that this question, this question, and a slew of other questions exist regarding the speed of mandelbrot rendering, but my question focuses on the differences in speed between C# and C++ render times. I also know that Ultra Fractal can render the mandelbrot set using MPFR and MPIR at iterations above 5000 in less than a second (atleast on my computer), meaning that it is possible to speed up fractal rendering. I know that Ultra Fractal uses many different algorithms for optimization, so I'm not expecting that kind of speed, but I find it hard to believe that their base code (their code without things like guessing pixel values) would take anywhere near as long as mine at 5 iterations. I want to know why my C++ code is so slow and how to optimize it. I don't need to optimize the C# code, because it won't be in the program when I'm done. It's just for testing. Some help speeding up the image saving would be nice since it takes ~600ms, but it is not included in the timings.

Timings

Note that these timings vary a little between tests. Also, the tests were conducted at 1440p

C#

5 iterations: ~450ms

250 iterations: ~1,100ms

C++

5 iterations: ~10,800ms

250 iterations: ~102,400ms (102 seconds)

As you can tell, the time it takes to render also does not go up at the same rate between the C++ and C# code. The time it takes the C# code to render at 250 iterations is ~3 times longer, while the C++ code takes ~10 times longer. I am not sure why this is the case.

Code

C#

using System;
using System.Threading;
using System.Threading.Tasks;
using System.Drawing;

namespace FractalToDesktop
{

    public class calculateFrame
    {
        private static int WIDTH, HEIGHT;
        private static int MAX_ITERATION;
        private static double xZoomPoint, yZoomPoint;
        private static double ZOOM_RATE;

        private static double realFactor;
        private static double imaginaryFactor;
        private static double ZCoordinateRealSquared, ZCoordinateImaginarySquared;
        private static int iteration;
        private static int maxIterationModifier;

        private Bitmap finalImg;
        LockBitmap img;
        private Color[] color;

        //These are the starting dimensions of the coordinate plane. 
        private static double viewportMinRealX = -2.5;
        private static double viewportMaxRealX = 2.5;
        private static double viewportMinImaginaryY = -1.75;
        private static double viewportMaxImaginaryY = 0; //This is calculated to prevent the screen from being streched. See calculation in constructor.

        public calculateFrame(double pixelation, int WIDTHin, int HEIGHTin, int MAX_ITERATIONin, double xZoomPointin, double yZoomPointin, double ZOOM_RATEin, string palette)
        {
            WIDTH = Convert.ToInt32(WIDTHin / pixelation);
            HEIGHT = Convert.ToInt32(HEIGHTin / pixelation);
            MAX_ITERATION = MAX_ITERATIONin;
            xZoomPoint = xZoomPointin;
            yZoomPoint = yZoomPointin;
            ZOOM_RATE = ZOOM_RATEin;

            maxIterationModifier = Convert.ToInt32(ZOOM_RATE * 50);

            finalImg = new Bitmap(WIDTH, HEIGHT, System.Drawing.Imaging.PixelFormat.Format32bppPArgb);
            img = new LockBitmap(finalImg);
            color = ColorPalette.coolerChooser(palette);

            //These are just constants defined for speeding up the calculation. Look at where they're called for the whole calculation. 
            viewportMaxImaginaryY = viewportMinImaginaryY + (viewportMaxRealX - viewportMinRealX) * HEIGHT / WIDTH; //This is calculated here because the height and width are now defined. See declaration for more information.
            realFactor = (viewportMaxRealX - viewportMinRealX) / (WIDTH - 1);
            imaginaryFactor = (viewportMaxImaginaryY - viewportMinImaginaryY) / (HEIGHT - 1);
        }

        public Image calculateOneFrame(double frame)
        {
            System.Diagnostics.Stopwatch stopwatch = System.Diagnostics.Stopwatch.StartNew();
            img.setImage(finalImg);

            for(int y = 0; y < HEIGHT; y ++)
            {
                double coordinateImaginary = yZoomPoint + (viewportMaxImaginaryY - (y * imaginaryFactor)) / Math.Pow(ZOOM_RATE, frame); //Convert the y pixel coordinate to a coordinate on the Imaginary plane

                for (int x = 0; x < WIDTH; x++)
                {
                    double coordinateReal = xZoomPoint + (viewportMinRealX + (x * realFactor)) / Math.Pow(ZOOM_RATE, frame); //Convert the x pixel coordinate to a coordinate on the Real plane

                    //Allows us to keep track of Z independent of starting point.
                    double ZCoordinateReal = coordinateReal;
                    double ZCoordinateImaginary = coordinateImaginary;

                    //This is done so that you don't have to square multiple times (optimization). 
                    ZCoordinateRealSquared = ZCoordinateReal * ZCoordinateReal;
                    ZCoordinateImaginarySquared = ZCoordinateImaginary * ZCoordinateImaginary;
                    iteration = 0;

                    while (ZCoordinateRealSquared + ZCoordinateImaginarySquared <= 4 && iteration < MAX_ITERATION)
                    {
                        // Calculate the sequence
                        ZCoordinateImaginary = (2 * ZCoordinateReal * ZCoordinateImaginary) + coordinateImaginary; //Calculates Z^2 + C for Imaginaries (2ab)
                        ZCoordinateReal = (ZCoordinateRealSquared - ZCoordinateImaginarySquared) + coordinateReal; //Calculates Z^2 + C for Reals (a^2 + b^2)

                        //Needs to be done inside the loop to update while condition. 
                        ZCoordinateRealSquared = ZCoordinateReal * ZCoordinateReal;
                        ZCoordinateImaginarySquared = ZCoordinateImaginary * ZCoordinateImaginary;
                        iteration++;
                    }

                    //If the number escapes, draw a black pixel. This the main body of the set.
                    if (iteration == MAX_ITERATION)
                    {
                        img.SetPixel(x, y, Color.FromArgb(0, 0, 0));
                    }
                    //This colors the edges of the set, which are the most interesting parts. Gradients can be chosen at the top of the runner. They are all stored in the ColorPalette class. You can add your own gradients easily there. 
                    else
                    {
                        img.SetPixel(x, y, color[iteration % color.Length]);
                    }
                }
            }
            img.UnlockBits();

            stopwatch.Stop();
            Console.WriteLine("Render Time: " + stopwatch.ElapsedMilliseconds);

            //MAX_ITERATION += maxIterationModifier;
            return finalImg;
        }
    }
}

C++

FrameCalculator.cpp

#include "FrameCalculator.h"

using mpfr::mpreal;
using namespace cimg_library;
using namespace std;

namespace FrameCalculator {

    ColorPalette::ColorPalette colorPaletteCreator;
    vector<vector<int>> colors;

    void calculateFrame::initializeCalculateFrame(double pixelation, int WIDTHin, int HEIGHTin, int MAX_ITERATIONin, double xZoomPointin, double yZoomPointin, double ZOOM_RATEin, const char* palette) {
        WIDTH = WIDTHin / pixelation;
        HEIGHT = HEIGHTin / pixelation;
        MAX_ITERATION = MAX_ITERATIONin;
        ZOOM_RATE = ZOOM_RATEin;
        mpfr_init_set_d(xZoomPoint, xZoomPointin, MPFR_RNDN);
        mpfr_init_set_d(yZoomPoint, yZoomPointin, MPFR_RNDN);

        mpfr_init(coordinateImaginary);
        mpfr_init(coordinateReal);
        mpfr_init(ZCoordinateReal);
        mpfr_init(ZCoordinateImaginary);
        mpfr_init(ZCoordinateRealSquared);
        mpfr_init(ZCoordinateImaginarySquared);

        mpfr_init(ZCoordinateRealSquaredPlusZCoordinateImaginarySquared);
        mpfr_init(ZCoordinateImaginaryTemp1);
        mpfr_init(ZCoordinateImaginaryTemp2);
        mpfr_init(ZCoordinateRealTemp1);

        mpfr_init(coordinateImaginaryTemp1);
        mpfr_init(coordinateImaginaryTemp2);
        mpfr_init(coordinateImaginaryTemp3);
        mpfr_init(coordinateRealTemp1);
        mpfr_init(coordinateRealTemp2);
        mpfr_init(coordinateRealTemp3);

        //This scales the viewport to the user's resolution to ensure the set isn't being streched. This is calculated here because the height and width are now defined. See declaration for more information.
        viewportMaxImaginaryY = viewportMinImaginaryY + (viewportMaxRealX - viewportMinRealX) * HEIGHT / WIDTH; 

        maxIterationModifier = ZOOM_RATEin * 50;

        colorPaletteCreator.ColorChooser(palette, colors);

        //These are just constants defined for speeding up the calculation. Look at where they're called for the whole calculation. 
        mpfr_init_set_d(realFactor, (double) (viewportMaxRealX - viewportMinRealX) / (WIDTH - 1), MPFR_RNDN);
        mpfr_init_set_d(imaginaryFactor, (double)(viewportMaxImaginaryY - viewportMinImaginaryY) / (HEIGHT - 1), MPFR_RNDN);
    }

    void calculateFrame::calculateOneFrame(int frame, const char* saveLocation) {
        chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();

        CImg<unsigned char> image(WIDTH, HEIGHT, 1, 3, 0);

        for (int y = 0; y < HEIGHT; y++) {
            mpfr_mul_si(coordinateImaginaryTemp1, imaginaryFactor, y, MPFR_RNDN);
            mpfr_d_sub(coordinateImaginaryTemp2, viewportMaxImaginaryY, coordinateImaginaryTemp1, MPFR_RNDN);
            mpfr_div_d(coordinateImaginaryTemp3, coordinateImaginaryTemp2, pow(ZOOM_RATE, frame), MPFR_RNDN);
            mpfr_add(coordinateImaginary, coordinateImaginaryTemp3, yZoomPoint, MPFR_RNDN);

            for (int x = 0; x < WIDTH; x++) {
                mpfr_mul_si(coordinateRealTemp1, realFactor, x, MPFR_RNDN);
                mpfr_add_d(coordinateRealTemp2, coordinateRealTemp1, viewportMinRealX, MPFR_RNDN);
                mpfr_div_d(coordinateRealTemp3, coordinateRealTemp2, pow(ZOOM_RATE, frame), MPFR_RNDN);
                mpfr_add(coordinateReal, coordinateRealTemp3, xZoomPoint, MPFR_RNDN);

                //Allows us to keep track of Z independent of starting point.
                mpfr_set(ZCoordinateReal, coordinateReal, MPFR_RNDN);
                mpfr_set(ZCoordinateImaginary, coordinateImaginary, MPFR_RNDN);

                //This is done so that you don't have to square multiple times (optimization). 
                mpfr_mul(ZCoordinateRealSquared, ZCoordinateReal, ZCoordinateReal, MPFR_RNDN);
                mpfr_mul(ZCoordinateImaginarySquared, ZCoordinateImaginary, ZCoordinateImaginary, MPFR_RNDN);

                iteration = 0;

                mpfr_add(ZCoordinateRealSquaredPlusZCoordinateImaginarySquared, ZCoordinateRealSquared, ZCoordinateImaginarySquared, MPFR_RNDN);

                while (mpfr_cmp_d(ZCoordinateRealSquaredPlusZCoordinateImaginarySquared, 4) <= 0 && iteration < MAX_ITERATION) {
                    //Calculate the sequence
                    //Calculates Z^2 + C for Imaginaries (2ab + C)
                    mpfr_mul(ZCoordinateImaginaryTemp1, ZCoordinateReal, ZCoordinateImaginary, MPFR_RNDN);
                    mpfr_mul_si(ZCoordinateImaginaryTemp2, ZCoordinateImaginaryTemp1, 2, MPFR_RNDN);
                    mpfr_add(ZCoordinateImaginary, ZCoordinateImaginaryTemp2, coordinateImaginary, MPFR_RNDN);

                    //Calculates Z^2 + C for Reals (a^2 + b^2 + C)
                    mpfr_sub(ZCoordinateRealTemp1, ZCoordinateRealSquared, ZCoordinateImaginarySquared, MPFR_RNDN);
                    mpfr_add(ZCoordinateReal, ZCoordinateRealTemp1, coordinateReal, MPFR_RNDN);

                    //Squares the two terms to update the while condition.
                    mpfr_mul(ZCoordinateRealSquared, ZCoordinateReal, ZCoordinateReal, MPFR_RNDN);
                    mpfr_mul(ZCoordinateImaginarySquared, ZCoordinateImaginary, ZCoordinateImaginary, MPFR_RNDN);
                    mpfr_add(ZCoordinateRealSquaredPlusZCoordinateImaginarySquared, ZCoordinateRealSquared, ZCoordinateImaginarySquared, MPFR_RNDN);
                    iteration++;
                }

                if (iteration == MAX_ITERATION) {
                    unsigned char black[3] = { 0,0,0 };
                    image.draw_point(x, y, black, 255);
                }
                else {
                    unsigned char color[3] = { colors[iteration % colors.size()][0], colors[iteration % colors.size()][1] , colors[iteration % colors.size()][2] };
                    image.draw_point(x, y, color, 255);
                }
            }
        }

        chrono::high_resolution_clock::time_point t2 = chrono::high_resolution_clock::now();
        auto duration = chrono::duration_cast<chrono::milliseconds>(t2 - t1).count();
        cout << "C++ Render Time: ";
        cout << duration << endl;

        image.save(saveLocation);
        //MAX_ITERATION += maxIterationModifier;
    }
}

FrameCalculator.h

#pragma once
#include <vector>
#include <iostream>
#include <gmp_util.h>
#include <cmath>
#include <chrono>
#include "CImg\CImg.h"
#include "ColorPalette.h"

using mpfr::mpreal;

namespace FrameCalculator {
    class calculateFrame {
    public:
        void initializeCalculateFrame(double pixelation, int WIDTHin, int HEIGHTin, int MAX_ITERATIONin, double xZoomPointin, double yZoomPointin, double ZOOM_RATEin, const char* palette);

        void calculateOneFrame(int frame, const char* saveLocation);
    private:
        int WIDTH, HEIGHT;
        int MAX_ITERATION;
        mpfr_t xZoomPoint, yZoomPoint;
        double ZOOM_RATE;

        mpfr_t realFactor, imaginaryFactor;
        mpfr_t ZCoordinateRealSquared, ZCoordinateImaginarySquared;
        mpfr_t coordinateImaginary, coordinateReal;
        mpfr_t ZCoordinateReal, ZCoordinateImaginary;

        //These are variables used in calculations because mpfr only supports calculations one at a time.
        mpfr_t ZCoordinateRealSquaredPlusZCoordinateImaginarySquared;
        mpfr_t ZCoordinateImaginaryTemp1;
        mpfr_t ZCoordinateImaginaryTemp2;
        mpfr_t ZCoordinateRealTemp1;

        mpfr_t coordinateImaginaryTemp1;
        mpfr_t coordinateImaginaryTemp2;
        mpfr_t coordinateImaginaryTemp3;
        mpfr_t coordinateRealTemp1;
        mpfr_t coordinateRealTemp2;
        mpfr_t coordinateRealTemp3;

        int iteration;
        int maxIterationModifier;

        //These are the starting dimensions of the coordinate plane. 
        double viewportMinRealX = -2.5;
        double viewportMaxRealX = 2.5;
        double viewportMinImaginaryY = -1.75;
        double viewportMaxImaginaryY; //This is calculated to prevent the screen from being streched. See calculation in constructor.
    };
}

Clarifications

I'm using the CImg library for drawing and saving in C++. I've left my method of timing the program to make sure that everyone knows what is being timed and what isn't. I've also left my comments in to improve readability. That being said, I'm not very satisfied with my C++ code's readability, so suggestions for readability that don't affect performance are appreciated. Also, the temp values in the c++ code are there because, as far as I can tell, MPFR doesn't support doing multiple lines of calculation at once while mixing doubles and mpfr_t types. You can look at the C# code to see the expression in one line (it's a lot easier to read there). A solution to having to use multiple variables for intermediate calculations is appreciated.

Interop

The C++ code isn't actually called from C++. I don't know whether this is relevant, but I thought I'd mention it, because I'm fairly new to C++/C# interop and I don't know how that affects performance. Basically, the C# code calls my C++ /clr project. Then that calls the code below, which has /clr turned off. The whole process takes 20 ms, and I've taken care to not include that in my timings. You can see my method of timing my code above.

Optimizations That I already Know

1. Multi-threading: I am going to handle multi-threading through a different method in C#. I just want to get my base code sorted out before that.
2. Using double precision at the beginning and switching to arbitrary precision: I plan to do this, but I want to make sure my arbitrary precision code works fine before I do that.
3. xZoomPoint and yZoomPoint can be doubles: I am aware of this, but in the future, I am going to implement an autopilot function, which will need more precision than double.

Answer 1

I realized that the answer probably does not lie within my code. There was a new method of calculating fractals using something called perturbation theory to significantly speed up calculations using arbritary precision. This allowed one to calculate images that used to take days in just a few days. Here is a MathSE question that goes into a little more detail on it and links to the paper itself. Once I'm done with my code, I'm going to comment it and post it on GitHub for everyone to use, but it the meantime, if you'd like a simple example, there is a program called Antelbrot on GitHub, which has some working code. I'm not 100% sure that Ultra Fractal uses perturbation theory, but it matches the results of the code.

Anyways, the speed up from perturbation theory allowed me to render the C++ code at 250 iterations in ~200ms or about 512 times faster. That speed increase only becomes more apparent as you zoom into the mandelbrot set or increase the iteration count.