mandelbrot/mandelbrot.cpp

/*
The GPLv3 License (GPLv3)

Copyright (c) 2022 Akos Horvath

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "mandelbrot.h"
#include "mpreal.h"
#include "utils.h"
#include <cstdint>
#include <iostream>
#include <limits>

using namespace mpfr;

Mandelbrotc::Mandelbrotc(Vec2mp const &f, Vec2mp const &t, Vec2i const &s, uint32_t mi)
                       : f(f), t(t), s(s), max_iter(mi)
{
    mpreal::set_default_prec(digits2bits(100));
    this->threads = std::vector<std::thread>();
    this->screen = new double*[s.x];
    for(int i = 0; i < s.x; ++i)
        this->screen[i] = new double[s.y];

    this->done = false;
    this->diff = mpreal(0.00001);
}

void Mandelbrotc::start_threads(bool first)
{
    if (!this->done && !first)
        return;
    Mandelbrotc const &m = *this;
    for (int i = 0; i < this->thread_count; i++) {
        std::thread t = std::thread(&Mandelbrotc::calc, this, i);
        this->threads.push_back(std::move(t));
        //std::cout << "thread created" << std::endl;
    }
    this->done = false;
}

void Mandelbrotc::stop_threads()
{
    for (int i = 0; i < this->thread_count; i++) {
        while (!this->threads[i].joinable())
            this->threads[i].join();
    }
}
void Mandelbrotc::calc(const uint8_t tid)
{
    mpreal x, y;
    std::vector<uint8_t> v = std::vector<uint8_t>();

    for(double i = tid; i < this->s.x; i+=this->thread_count) {
        for(double j = 0; j < this->s.y; j++) {
            x = (mpfr::abs(this->f.x - this->t.x) * (i / this->s.x)) + this->f.x;
            y = (mpfr::abs(this->f.y - this->t.y) * (j / this->s.y)) + this->f.y;

            this->screen[(int)i][(int)j] = this->mandelbrot(Vec2mp(x, y));
        }
    }

    this->updatePrec();

    //std::cout << "thread " << tid << " done" << std::endl;
    if (tid == thread_count-1)
        this->done = true;
}

double Mandelbrotc::mandelbrot(const Vec2mp &n)
{
    double iter = 0.0;

    // BigFloat x = BigFloat(n.x.precision);
    // x.setValue(0.0);
    // BigFloat y = BigFloat(n.y.precision);
    // y.setValue(0.0);

    mpreal x = mpreal(0.0);
    mpreal y = mpreal(0.0);
    mpreal x2 = mpreal(0.0);
    mpreal y2 = mpreal(0.0);
    mpreal xtemp = mpreal(0.0);

    // while (x2 + y2 <= (1 << 16) && iter < max_iter) {
    //     y = 2.0 * x * y + n.y;
    //     x = x2 - y2 + n.x;
    //     x2 = x * x;
    //     y2 = y * y;

    //     iter++;
    // }

    while (x*x + y*y <= (1 << 16) && iter < max_iter) {
        xtemp = x*x - y*y + n.x;
        y = 2*x*y + n.y;
        x = xtemp;

        iter++;
    }

    if (iter == max_iter)
        return max_iter;

    double log_zn = (mpfr::log(x*x + y*y) / 2).toDouble();
    double nu = log(log_zn / log(2)) / log(2);

    iter = iter + 1 - nu;

    return iter;
    //return iter + 1 - log(log2(x2.toDouble()+y2.toDouble()));
}

void Mandelbrotc::updatePrec()
{
    //std::cout << "before updatePrec: " << this->f.x.get_prec() << std::endl;
    mpfr_prec_t bprec = this->f.x.get_prec();
    mpfr_prec_t aprec;

    std::cout << "before updatePrec def: " << mpfr::bits2digits(
              bprec) << std::endl;

    if (mpfr::abs(mpfr::abs(this->f.x) - mpfr::abs(this->t.x)) < diff ||
        (mpfr::abs(this->f.x) - mpfr::abs(this->t.y)) < diff) {

        aprec = mpfr::digits2bits(mpfr::bits2digits(bprec)+5);

        std::cout << "true" << std::endl;
    }
    else {
        return;
    }

    this->f.x.set_prec(aprec);
    this->f.y.set_prec(aprec);
    this->t.x.set_prec(aprec);
    this->t.y.set_prec(aprec);

    //std::cout << "after updatePrec: " << this->f.x.get_prec() << std::endl;
    std::cout << "after updatePrec def: " << mpfr::bits2digits(
            this->f.x.get_prec()) << std::endl;

    diff = mpfr::abs(mpfr::abs(this->f.x) - mpfr::abs(this->t.x));
    diff *= 0.01;
}