Gfrktyl

#include <vector>

template<typename T, typename SOURCE, typename DEST, std::size_t SIZE=8>

void calculate(const T &_dLimit, const unsigned int _iMax, SOURCE &_rSource, DEST &_rDest)

{       // starting point -- does not change during iteration

    float adX0[SIZE], adY0[SIZE];

        // coordinates of point -- for reporting back number of iterations

    unsigned int aiX[SIZE], aiY[SIZE];

        // current value of iteration

    float adX[SIZE], adY[SIZE];

        // number of iterations

    unsigned int aiIt[SIZE];

        // is the point a valid point

    bool ab[SIZE];

        // number of valid points currently being calculated

    unsigned int iValid = 0;

        // initial -- get SIZE points

    for (std::size_t i = 0; i < SIZE; ++i)

    {   const auto s = _rSource.getPoint();

        if (!(ab[i] = s.b))

            if (iValid == 0)

                return;

            else

                ;

        else

            ++iValid;

            // the coordinates of the point

        aiX[i] = s.ix;

        aiY[i] = s.iy;

            // the complex value of the point

        adX[i] = adX0[i] = s.dx;

        adY[i] = adY0[i] = s.dy;

            // number of iterations

        aiIt[i] = 0;

    }

    while (1)

    {       // square values

        float ad2X[SIZE], ad2Y[SIZE];

            // stop condition fulfilled?

        bool abl[SIZE];

            // x square

        for (std::size_t i = 0; i < SIZE; ++i)

            ad2X[i] = adX[i]*adX[i];

            // y square

        for (std::size_t i = 0; i < SIZE; ++i)

            ad2Y[i] = adY[i]*adY[i];

            // stop condition met

        for (std::size_t i = 0; i < SIZE; ++i)

            if (aiIt[i] > _iMax || ad2Y[i] + ad2X[i] > _dLimit)

                abl[i] = false;

            else

                abl[i] = true;

            // new x value

        for (std::size_t i = 0; i < SIZE; ++i)

            ad2X[i] = ad2X[i] - ad2Y[i] + adX0[i];

            // new y value

        for (std::size_t i = 0; i < SIZE; ++i)

            ad2Y[i] = 2*adX[i]*adY[i] + adY0[i];

            // copy over to x

        for (std::size_t i = 0; i < SIZE; ++i)

            adX[i] = ad2X[i];

            // copy over to y

        for (std::size_t i = 0; i < SIZE; ++i)

            adY[i] = ad2Y[i];

            // if condition is met, report results and request new value

        for (std::size_t i = 0; i < SIZE; ++i)

            if (!abl[i] && ab[i])

            {   _rDest.report(aiX[i], aiY[i], aiIt[i]);

                const auto s = _rSource.getPoint();

                if (!(ab[i] = s.b))

                    if (!--iValid)

                        return;

                aiX[i] = s.ix;

                aiY[i] = s.iy;

                adX[i] = adX0[i] = s.dx;

                adY[i] = adY0[i] = s.dy;

                aiIt[i] = 0;

            }

            else

                ++aiIt[i];

    }

}

struct point

{   float dx, dy;

    unsigned int ix, iy;

    bool b;

};

struct source

{   point getPoint(void);

    void report(unsigned int, unsigned int, unsigned int);

};

int main()

{

    source sSource;

    calculate(10.0, 100, sSource, sSource);

}