#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <cstdint>
#include <gsl/gsl_poly.h>

unsigned long long cur_root;
unsigned long long capacity;

void polynomial_recursive(double **roots, double min_x, double max_x, int min_degree, int max_degree, int cur_degree, double *coeff)
{
    for (double c = cur_degree==max_degree?1:min_x; c <= max_x; c++)
    {
            double *new_coeff;

            new_coeff = (double *) malloc (sizeof (double) * (cur_degree + 1));
            if (new_coeff == NULL)
            {
                    perror ("malloc failed");
            }
            memcpy (new_coeff, coeff, cur_degree * sizeof (double));
            new_coeff[cur_degree] = c;

            if ( cur_degree == max_degree )
            {
                    double z[cur_degree*2];

                    gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (cur_degree + 1);
                    gsl_poly_complex_solve (new_coeff, cur_degree + 1, w, z);
                    gsl_poly_complex_workspace_free (w);

                    for (int i = 0; i < cur_degree; i++)
                    {
                            if (cur_root + 2 >= capacity)
                            {
                                    capacity *= 2;
                                    double *tmp = (double *) realloc (*roots, capacity * sizeof (double));
                                    if (tmp == NULL)
                                    {
                                            perror ("realloc failed");
                                    }
                                    else
                                    {
                                            *roots = tmp;
                                    }
                            }
                            (*roots) [cur_root++] = z[2*i];
                            (*roots) [cur_root++] = z[2*i+1];
                    }
            } else {
                    #pragma omp task shared(new_coeff)
                    polynomial_recursive (roots, min_x, max_x, min_degree, max_degree, cur_degree + 1, new_coeff);
                    #pragma omp taskwait
            }
            free (new_coeff);
    }
}

unsigned long long getDensityForPoint(double x, double y, double *points, int size, double radius)
{
    unsigned long long count = 0;
    for (int i = 0; i < size; i += 2)
      {
        double x2 = points[i];
        double y2 = points[i+1];
        if (x != x2 || y != y2)
          {
            double dist = fabs(x2 - x) + fabs(y2 - y);;
            if (dist <= radius)
              {
                count++;
              }
          }
      }
    return count;
}

double
calculateDensityScalarFromDensity(unsigned long min, unsigned long max, unsigned long density)
{
    return 0.0 + ((1.5 - 0.0) / (double)(max - min)) * (density - min);
}

int
main(void)
{
    double *roots;

    capacity = 100;
    cur_root = 0;
    roots = (double *) malloc (capacity * sizeof (double));
    if (roots == NULL)
      {
        perror ("malloc failed");
      }

    #pragma omp parallel
    {
        #pragma omp single
        {
            polynomial_recursive (&roots, -4, 4, 2, 5, 0, NULL);
        }
    }

    printf("%llu\n", cur_root);

    int w=1920,h=1080;

    uint64_t *buckets = (uint64_t*)calloc(w*h,sizeof(buckets));
    double min_x,max_x,min_y,max_y;
    min_x = max_x = roots[0];
    min_y = max_y = roots[1];

    for ( int i = 0; i < cur_root; i+=2 )
    {
            min_x = roots[i]<min_x?roots[i]:min_x;
            max_x = roots[i]>max_x?roots[i]:max_x;
            min_y = roots[i+1]<min_y?roots[i+1]:min_y;
            max_y = roots[i+1]>max_y?roots[i+1]:max_y;
    }

    printf("range: (%f,%f)-(%f,%f)\n",min_x,min_y,max_x,max_y);

    double pad_x = (max_x-min_x)/100;
    double pad_y = (max_y-min_y)/100;
    min_x -= pad_x;
    max_x += pad_x;
    min_y -= pad_y;
    max_y += pad_y;

    for ( int i = 0; i < cur_root; i+=2 )
    {
            int x = (roots[i]-min_x)/(max_x-min_x)*w;
            int y = (roots[i+1]-min_y)/(max_y-min_y)*h;
            buckets[x+w*y]++;
    }

    uint64_t v_min,v_max;
    v_max = v_min = buckets[0];
    for ( int i = 0; i < w*h; i++ )
    {
            v_min = buckets[i]<v_min?buckets[i]:v_min;
            v_max = buckets[i]>v_max?buckets[i]:v_max;
    }

    printf("max %llu, min %llu\n", v_max, v_min);
    FILE *gnuplotPipe = popen ("gnuplot -persistent", "w");
    FILE *temp = fopen("data.temp", "w");

    for ( int i = 0; i < w*h; i++ )
    {
            uint64_t val = buckets[i];
            if ( !val ) continue;
            int x = i%w;
            int y = i/w;
            double fx = min_x + (max_x-min_x) * x / w;
            double fy = min_y + (max_y-min_y) * y / h;
            fprintf (temp, "%lf %lf %lf\n", fx, fy, calculateDensityScalarFromDensity(v_min, v_max, val));
    }

    free(buckets);
    free(roots);

    fprintf (gnuplotPipe, "set zrange [0.0:1.5]\nset cbrange [0.0:1.5]\nset term pngcairo\nset terminal png size %d,%d\nset output \"plot.png\"\nset palette defined (0 \"black\", 0.5 \"red\", 1 \"yellow\", 1.5 \"white\")\nplot 'data.temp' with points pointtype 5 pointsize 0.001 palette\n",w,h);
}