#include <stdint.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>

#include <quadmath.h>

#include <omp.h>

#define PREFIX "img2-"

#define VENTANA_ANCHO   (1920*2)
#define VENTANA_ALTO    (1080*2)
#define ANCHO   VENTANA_ANCHO
#define ALTO    VENTANA_ALTO
#define MAX_ITER 10000 //500 en video2

//__float128 xc = -1.740119020442992652289958266376L;
//__float128 yc = 0.028019209975637676263081582255L;
//__float128 xc = 0.264226442163892496529931097626331393L, yc = 0.002572261418806922467217290773078275;
__float128 xc = -1.749900374662927212679117139337847675Q, yc = 0.000000000002454985822810237255224056Q;

uint32_t computar_intensidad(__float128 x0, __float128 y0, __float128 zoom) {
    //x0 = x0 * zoom - 1.7400623825793398724575;  // Derecha: Achicar
    //y0 = y0 * zoom + 0.0281753397792111;  // Abajo: Achicar
    x0 = x0 * zoom + xc;  // Derecha: Achicar
    y0 = y0 * zoom + yc;  // Abajo: Achicar
    //x0 = x0 * zoom - 1.78;  // Derecha: Achicar
    //y0 = y0 * zoom + 0.0;  // Abajo: Achicar
    __float128 x = 0;
    __float128 y = 0;
    int n = 0;
    while(x*x + y*y <= 4 && n < MAX_ITER) {
        __float128 xtemp = x*x - y*y + x0;
        y = 2*x*y + y0;
        x = xtemp;
        n++;
    }

    if(n == MAX_ITER)
        return 0;

    //float h = (n + 1 - log(log2(sqrt(x*x+y*y)))) / MAX_ITER * 360;
    //float h = n * 360.0 / MAX_ITER;
    float h = (int)(n * 360.0 / 100) % 360;
    float xxx = (1 - fabs(fmodf(h / 60.0, 2) - 1));
    uint8_t xx = (xxx < 1 ? xxx : 1) * 255;
    uint8_t r = 0, g = 0, b = 0;

    if(h < 60)
        r = 255, g = xx;
    else if(h < 120)
        g = 255, r = xx;
    else if(h < 180)
        g = 255, b = xx;
    else if(h < 240)
        b = 255, g = xx;
    else if(h < 300)
        b = 255, r = xx;
    else
        r = 255, b = xx;

    return (r << 24) | (g << 16) | (b << 8);
}


int main() {
    __float128 zoom = 0.179;
    double factor = 0.99;

    size_t n = 0;
    for(double z = zoom; z >= 1e-21; z *= factor)
        n++;

    double *zooms = malloc(n * sizeof(double));
    zooms[0] = zoom;
    for(size_t i = 1; i < n; i++)
        zooms[i] = zooms[i - 1] * factor;

    zooms[0] = zooms[4160] * 0.5;
    printf("%.16e\n", zooms[4160]);
    n = 6;
    for(size_t i = 1; i < n; i++)
        zooms[i] = zooms[i - 1] * 0.5;

    #pragma omp parallel for schedule(dynamic)
    for(int i = 0; i < n; i++) {
        printf("%d %.16e\n", i, zooms[i]);

        char aux[100];
        sprintf(aux, "%s%05d.ppm", PREFIX, i);
        FILE *f = fopen(aux, "wt");
        fprintf(f, "P3\n%d %d\n255\n", ANCHO, ALTO);


        for(int vy = ALTO / 2; vy > - ALTO / 2; vy--)
            for(int vx = - ANCHO / 2; vx < ANCHO / 2; vx++) {
                uint32_t c = computar_intensidad(vx, vy, zooms[i]);
                fprintf(f, "%d %d %d\n", c >> 24, (c >> 16) & 0xFF, (c >> 8) & 0xFF);
            }
        fclose(f);

        sprintf(aux, "convert %s%05d.ppm %s%05d.png; rm %s%05d.ppm", PREFIX, i, PREFIX, i, PREFIX, i);
        system(aux);
    }

    free(zooms);

//    system("ffmpeg -y -framerate 24 -i " PREFIX "%05d.png -c:v libx264 -crf 0 output.mp4");

    return 0;
}