Reverse Julia^2 Iteration Mutation...

complex_t
ct_reversex(
    complex_t o, // center
    xfloat_t r, // radius
    complex_t z, // start point
    complex_t c, // const point
    xfloat_t p, // prob factor wrt roots
    glm::vec4 const& color, // color
    unsigned int imax // iterations
){
    xfloat_t smax = 0.0;
    complex_t cx = { -c.real(), c.imag() };
    xfloat_t P = 0.2;

    for (unsigned int i = 0; i < imax; ++i)
    {
        complex_t d = z - c; // the master difference

        xfloat_t rn0 = RAND_FLOAT(); // [0...1] prob range
        xfloat_t rn1 = RAND_FLOAT(); // [0...1] prob range

        xfloat_t l = std::abs(d); // distance
        xfloat_t s = std::sqrt(l); // sqrt of distance

        // Buddha like mutation
        if (rn0 > 0.97) // slim prob factor... ;^)
        {
            // The log of distance extends the fractal
            // internal and external structures! Wow.
            s = std::log(l);
        }

        xfloat_t a = std::atan2(d.imag(), d.real()) / 2.0; // angle
        xfloat_t sgn = -1.0f;

        smax = std::max(smax, s); // gain largest radius


        // notice the (a + PI) here... Simple, but powerful indeed ;^)
        // gain both roots
        complex_t nz_0 = { s * std::cos(a), s * std::sin(a + PI) };
        complex_t nz_1 = { -1.0 * s * std::cos(a), -1.0 * s * std::sin(a) };

        // project roots onto the canvas (unsigned integer plane)
        complex_t nz_p0 = {
            std::floor(nz_0.real() * r + o.real()),
            std::floor(nz_0.imag() * r + o.imag())
        };

        complex_t nz_p1 = {
            std::floor(nz_1.real() * r + o.real()),
            std::floor(nz_1.imag() * r + o.imag())
        };

        // plot the damn thing in the hit list
        if (nz_p0.real() > -1 && nz_p0.real() < m_gc.m_width &&
            nz_p0.imag() > -1 && nz_p0.imag() < m_gc.m_height)
        {
            unsigned int ihit = nz_p0.real() + m_gc.m_width * nz_p0.imag();
            if (ihit < m_hits.size()) m_hits[ihit] += 1.0f;
        }

        if (nz_p1.real() > -1 && nz_p1.real() < m_gc.m_width &&
            nz_p1.imag() > -1 && nz_p1.imag() < m_gc.m_height)
        {
            unsigned int ihit = nz_p1.real() + m_gc.m_width * nz_p1.imag();
            if (ihit < m_hits.size()) m_hits[ihit] += 1.0f;
        }

        // choose random root with probability factor p
        z = (rn1 > p) ? nz_0 : nz_1;

        if (!(i % 1000000))
        {
            std::printf("iterated: %u of %u\r", i, imax);
        }
    }

    std::printf("\r\nComplete!\r\n");

    return z;
}