The Computer Language
23.03 Benchmarks Game

mandelbrot C++ g++ program

source code

// The Computer Language Benchmarks Game
// https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
//
// Contributed by Markus Flad.
// Unlike the original version, this version uses std::array to make the
// compiler vectorize.
//
// Compile with following g++ flags
//  -std=c++17 -O3 -Wall -march=native -mno-fma

#include <string>
#include <iostream>
#include <vector>
#include <complex>
#include <algorithm>
#include <thread>
#include <climits>
#include <version>
#include <array>

const auto numberOfCpuCores = std::thread::hardware_concurrency();

// The PortableBinaryBitmap manages access to the pbm output file and provides
// interlaced canvases that allow threads to write to the bitmap in parallel.
class PortableBinaryBitmap {
public:
    using Size=std::size_t;
    PortableBinaryBitmap(std::ostream& ostr, Size width, Size height)
    : _ostr (ostr)
    , _width (roundToMultiple(width, CHAR_BIT))
    , _height (roundToMultiple(height, numberOfCpuCores))
    , _data ((_width * _height) / CHAR_BIT) {
        _ostr << "P4" << '\n';
        _ostr << _width << ' ' << _height << '\n';
    }
    ~PortableBinaryBitmap() {
        _ostr.write(_data.data(), _data.size());
    }
    Size width() const noexcept {
        return _width;
    }
    Size height() const noexcept {
        return _height;
    }
    Size widthInBytes() const noexcept {
        return _width / CHAR_BIT;
    }
    struct Line {
        constexpr static Size pixelsPerWrite() noexcept {
            return CHAR_BIT;
        }
        Size y;
        Size width;
        char* data;
    };
    // The InterlacedCanvas provides interlaced access to the bitmap data. Each
    // thread must use its own InterlacedCanvas to write to the bitmap.
    class InterlacedCanvas {
    public:
        class Iterator {
        public:
            Iterator(Size y, Size _width, char* data,
                    Size interlaceIncrement, Size dataPointerIncrement) noexcept
            : _il {y, _width, data}
            , _interlaceIncrement (interlaceIncrement)
            , _dataPointerIncrement (dataPointerIncrement) {
            }
            Line& operator*() noexcept {
                return _il;
            }
            bool operator!=(const Iterator& other) const noexcept {
                return _il.data != other._il.data;
            }
            Iterator& operator++() noexcept {
                _il.y += _interlaceIncrement;
                _il.data += _dataPointerIncrement;
                return *this;
            }
        private:
            Line _il;
            Size _interlaceIncrement;
            Size _dataPointerIncrement;
        };
        InterlacedCanvas(PortableBinaryBitmap& pbm, Size yStart, Size increment)
                noexcept
        : _pbm (pbm)
        , _yStart (yStart)
        , _increment (increment)
        , _dataStart (yStart * pbm.widthInBytes())
        , _dataPointerIncrement (increment * pbm.widthInBytes()) {
        }
        Size width() const noexcept {
            return _pbm.width();
        }
        Size height() const noexcept {
            return _pbm.height();
        }
        Iterator begin() noexcept {
            return Iterator(_yStart, _pbm.width(),
                    _pbm._data.data() + _dataStart,
                    _increment, _dataPointerIncrement);
        }
        Iterator end() noexcept {
            return Iterator(_yStart + _pbm.height(), _pbm.width(),
                    _pbm._data.data() + _pbm._data.size() + _dataStart,
                    _increment, _dataPointerIncrement);
        }
    private:
        PortableBinaryBitmap& _pbm;
        Size _yStart;
        Size _increment;
        Size _dataStart;
        Size _dataPointerIncrement;
    };
    std::vector<InterlacedCanvas> provideInterlacedCanvas(Size increment)
            noexcept {
        std::vector<InterlacedCanvas> interlacedCanvasVector;
        for (Size yStart=0; yStart<increment; ++yStart) {
            interlacedCanvasVector.emplace_back(*this, yStart, increment);
        }
        return interlacedCanvasVector;
    }
    static Size roundToMultiple (Size number, Size base) noexcept {
        return number + ((number % base) ? (base - number % base) : 0);
    }
private:
    std::ostream& _ostr;
    Size _width;
    Size _height;
    std::vector<char> _data;
};

template<typename NUMBER_TYPE>
class VectorizedNumber
{
public:
    constexpr static std::size_t SIZE = 8;
    using NumericArray = std::array<NUMBER_TYPE, SIZE>;

    VectorizedNumber() noexcept
    : _x(_values.data()) {
    }
    explicit VectorizedNumber(NUMBER_TYPE value) noexcept
    : _x(_values.data()) {
        std::fill(_values.begin(), _values.end(), value);
    }
    VectorizedNumber(const VectorizedNumber& other) noexcept
    : _x(_values.data()) {
        *this = other;
    }
    VectorizedNumber& operator=(const VectorizedNumber& other) noexcept {
        // In GCC (unlike Clang), using the std::array assignment operator and
        // also its copy constructor slows down SIMD performance. Therefore we
        // use a raw loop here.
        for (size_t i=0; i<SIZE; ++i) {
            _values[i] = other._values[i];
        }
        return *this;
    }
    NUMBER_TYPE operator[](std::size_t i) const noexcept{
        return _values[i];
    }
    NUMBER_TYPE& operator[](std::size_t i) noexcept {
        return _values[i];
    }
    typename NumericArray::const_iterator begin() const noexcept {
        return _values.begin();
    }
    typename NumericArray::const_iterator end() const noexcept {
        return _values.end();
    }
    bool operator>(NUMBER_TYPE value) const noexcept {
        return (std::all_of(_values.begin(), _values.end(),
                [&value](NUMBER_TYPE v) {return v > value;}));
    }
    char lteToPixels(NUMBER_TYPE threshold) const noexcept {
        char result = 0;
        if (_values[0] <= threshold) result |= 0b10000000;
        if (_values[1] <= threshold) result |= 0b01000000;
        if (_values[2] <= threshold) result |= 0b00100000;
        if (_values[3] <= threshold) result |= 0b00010000;
        if (_values[4] <= threshold) result |= 0b00001000;
        if (_values[5] <= threshold) result |= 0b00000100;
        if (_values[6] <= threshold) result |= 0b00000010;
        if (_values[7] <= threshold) result |= 0b00000001;
        return result;
    }
private:
    NumericArray _values;
    // A pointer to the data, which is not used outside, but helps the compiler
    // with the SIMD optimization (especially GCC 11).
    NUMBER_TYPE* _x;
};

// VectorizedComplex provides a convenient interface to deal with complex
// numbers and uses the power of SIMD for high execution speed.
template <typename NUMBER_TYPE>
class VectorizedComplex {
public:
    using Size = std::size_t;

    VectorizedComplex() = default;
    VectorizedComplex(const VectorizedComplex&) = default;
    VectorizedComplex& operator=(const VectorizedComplex&) = default;
    VectorizedComplex(const VectorizedNumber<NUMBER_TYPE>& reals,
            NUMBER_TYPE commonImagValue) noexcept
    : _reals(reals)
    , _imags(commonImagValue){
    }
    VectorizedComplex& squareAndAdd(const VectorizedComplex& c,
            VectorizedNumber<NUMBER_TYPE>& squaredAbs) noexcept {
        for (Size i=0; i<VectorizedNumber<NUMBER_TYPE>::SIZE; ++i) {
            auto realSquared = _reals[i] * _reals[i];
            auto imagSquared = _imags[i] * _imags[i];
            auto realTimesImag = _reals[i] * _imags[i];
            _reals[i] = realSquared - imagSquared + c._reals[i];
            _imags[i] = realTimesImag + realTimesImag + c._imags[i];
            squaredAbs[i] = realSquared + imagSquared;
        }
        return *this;
    }
private:
    VectorizedNumber<NUMBER_TYPE> _reals;
    VectorizedNumber<NUMBER_TYPE> _imags;
};

// The ComplexPlaneCalculator performs function f(c), with c as a
// VectorizedComplex and a byte as the return value. Due to its eightfold
// vectorization, each returned bit can return a Boolean value from the
// calculation f(c). The full byte is then written to the canvas. This is done
// until the whole bitmap is filled.
template <typename NUMBER_TYPE, class Functor>
class ComplexPlaneCalculator {
public:
    using VComplex = VectorizedComplex<NUMBER_TYPE>;
    using Line = typename PortableBinaryBitmap::Line;
    using Size = std::size_t;

    ComplexPlaneCalculator(const std::complex<NUMBER_TYPE>& cFirst,
            const std::complex<NUMBER_TYPE>& cLast,
            PortableBinaryBitmap::InterlacedCanvas& canvas, Functor f) noexcept
    : _cFirst(cFirst)
    , _cLast(cLast)
    , _canvas(canvas)
    , _f(f) {
    }
    void operator()() noexcept {
        const NUMBER_TYPE realRange = _cLast.real() - _cFirst.real();
        const NUMBER_TYPE imagRange = _cLast.imag() - _cFirst.imag();
        const NUMBER_TYPE rasterReal = realRange / _canvas.width();
        const NUMBER_TYPE rasterImag = imagRange / _canvas.height();
        std::vector<VectorizedNumber<NUMBER_TYPE>> cRealValues;
        cRealValues.reserve(_canvas.width() / Line::pixelsPerWrite());
        for (Size x=0; x<_canvas.width(); x+=Line::pixelsPerWrite()) {
            VectorizedNumber<NUMBER_TYPE> cReals;
            for (Size i=0; i<Line::pixelsPerWrite(); ++i) {
                cReals[i] = _cFirst.real() + (x+i)*rasterReal;
            }
            cRealValues.push_back(cReals);
        }
        for (Line& line : _canvas) {
            char* nextPixels = line.data;
            char lastPixels = 0x00;
            const NUMBER_TYPE cImagValue = _cFirst.imag() + line.y*rasterImag;
            for (const VectorizedNumber<NUMBER_TYPE>& cReals : cRealValues) {
                const VComplex c(cReals, cImagValue);
                *nextPixels = _f(c, lastPixels);
                lastPixels = *nextPixels;
                nextPixels++;
            }
        }
    }
private:
    std::complex<NUMBER_TYPE> _cFirst;
    std::complex<NUMBER_TYPE> _cLast;
    PortableBinaryBitmap::InterlacedCanvas _canvas;
    Functor _f;
};

// Functor calculating a Mandelbrot iteration for a VectorizedComplex. This
// means that for eight complex numbers the Mandelbrot calculation is
// (potentially) executed in parallel. The result is a byte that contains a 1
// for each bit if the corresponding complex number is in the Mandelbrot set,
// and a 0 if it is not.
template <typename NUMBER_TYPE>
class MandelbrotFunction {
public:
    using VComplex = VectorizedComplex<NUMBER_TYPE>;
    using Size = std::size_t;
    constexpr static Size ITERATIONS_WITHOUT_CHECK = 5;
    constexpr static char NONE_IN_MANDELBROT_SET = 0x00;

    MandelbrotFunction(Size maxIterations, NUMBER_TYPE pointOfNoReturn = 2.0)
            noexcept
    : _maxOuterIterations(maxIterations / ITERATIONS_WITHOUT_CHECK - 2)
    , _squaredPointOfNoReturn(pointOfNoReturn * pointOfNoReturn) {
    }
    static void doMandelbrotIterations(VComplex& z, const VComplex& c,
            VectorizedNumber<NUMBER_TYPE>& squaredAbs) noexcept {
        for (Size j=0; j<ITERATIONS_WITHOUT_CHECK; ++j) {
            z.squareAndAdd(c, squaredAbs);
        }
    }
    char operator()(const VComplex& c, char lastPixels) const noexcept {
        VComplex z = c;
        VectorizedNumber<NUMBER_TYPE> squaredAbs;
        if (lastPixels == NONE_IN_MANDELBROT_SET) {
            for (Size i=0; i<_maxOuterIterations; ++i) {
                doMandelbrotIterations(z, c, squaredAbs);
                if (squaredAbs > _squaredPointOfNoReturn) {
                    return NONE_IN_MANDELBROT_SET;
                }
            }
        } else {
            for (Size i=0; i<_maxOuterIterations; ++i) {
                doMandelbrotIterations(z, c, squaredAbs);
            }
        }
        doMandelbrotIterations(z, c, squaredAbs);
        doMandelbrotIterations(z, c, squaredAbs);
        return squaredAbs.lteToPixels(_squaredPointOfNoReturn);
    }
private:
    Size _maxOuterIterations;
    NUMBER_TYPE _squaredPointOfNoReturn;
};

int main(int argc, char** argv) {
    using NumberType = double;
    using ComplexNumber = std::complex<NumberType>;
    using MandelbrotCalculator = ComplexPlaneCalculator<NumberType,
            MandelbrotFunction<NumberType>>;
    std::size_t n = 16000;
    if (argc>=2) {
        n = atoi(argv[1]);
    }
    const std::size_t maxIterations = 50;
    PortableBinaryBitmap pbm(std::cout, n, n);
    auto canvasVector = pbm.provideInterlacedCanvas(numberOfCpuCores);
    std::vector<std::thread> threads;
    for (auto& canvas : canvasVector) {
        threads.emplace_back(MandelbrotCalculator (ComplexNumber(-1.5, -1.0),
                ComplexNumber(0.5, 1.0), canvas,
                MandelbrotFunction<NumberType> (maxIterations)));
    }
    for (auto& t : threads) {
        t.join();
    }
    return 0;
}
    

notes, command-line, and program output

NOTES:
64-bit Ubuntu quad core
Ubuntu 12.2.0-3ubuntu1


Wed, 25 Jan 2023 02:28:05 GMT

MAKE:
/usr/bin/g++ -c -pipe -O3 -fomit-frame-pointer -march=ivybridge  --std=c++17 -I ./Include/gpp/vcl mandelbrot.gpp-0.c++ -o mandelbrot.gpp-0.c++.o &&  \
        /usr/bin/g++ mandelbrot.gpp-0.c++.o -o mandelbrot.gpp-0.gpp_run -pthread 
rm mandelbrot.gpp-0.c++

4.72s to complete and log all make actions

COMMAND LINE:
./mandelbrot.gpp-0.gpp_run 16000

(BINARY) PROGRAM OUTPUT NOT SHOWN