multijulia for Paul

namespace ct_multi_julia_experimental
{
    typedef std::complex<float> ct_complex;

    void manifest(
        ct::plot::cairo::plot_2d& plot
    ) {
        std::cout << "\ct_multi_julia_experimental::manifest()\n";
        std::cout << "__________________________________________\n\n";

        {
            unsigned long n = 20000000;

            std::array<ct_complex, 3> julia_sets;

            float radius = 3.f;

            {

                float angle_base = CT_PI2 / julia_sets.size();
                float radius_base = radius / julia_sets.size();

                for (unsigned long i = 0; i < julia_sets.size(); ++i)
                {
                    float angle = angle_base * i;
                    float radius = radius_base * i;

                    ct_complex p0_normal = {
                        std::cos(angle + std::sin(angle * 1.f)),
                        std::sin(angle)
                    };

                    p0_normal *= radius;

                    std::cout << "[" << i << "] = " << p0_normal << "\n";

                    julia_sets[i] = p0_normal;
                }
            }

            ct_complex z = { 4, 4 };

            for (unsigned long i = 0; i < n; ++i)
            {
                // gain random's
                std::size_t rand_julia_set = std::rand() % julia_sets.size();
                float rand_root = std::rand() / (RAND_MAX - 0.f);

                // choose julia set
                ct_complex& julia = julia_sets[rand_julia_set];

                // root logic
                ct_complex dif = z - julia;

                ct_complex root = { 0, 0 };

                // select roots
                if (rand_root < 0.5f)
                {
                    root = std::sqrt(dif);
                }

                else
                {
                    root = -std::sqrt(dif);
                }


                // plot
                {
                    glm::vec2 point = { z.real(), z.imag() };
                    //plot.set_pixelf(point, CT_RGBF(1, 0, 0));


                    ct::plot::cairo::pixel* pixel = plot.get_pixelf(point);

                    if (pixel)
                    {
                        if (pixel->r < 255 - 16)
                        {
                            pixel->r += 16;
                        }

                        pixel->g += 1;
                        pixel->b += 5;
                    }

                }


                // update loop
                z = root;


                //std::cout << "[" << i << "] = " << julia << "\n";
                //std::cout << "[" << i << "] = " << root << "\n";
            }
        }
    }
}