#include <cmath>
#include <cstdlib>
#include <iostream>
#include <vector>

#include <mpfr.h>

using namespace std;

int main(int argc, char **argv)
{
  if (argc < 5)
  {
    cerr << "usage: " << argv[0] << " re im zoom iters" << endl;
    return 1;
  }
  int width = atoi(getenv("COLUMNS"));
  int height = atoi(getenv("LINES")) - 1;
  double zoom = atof(argv[3]);
  int iterations = atoi(argv[4]);
  int precision = log2(zoom * height) + 16;
  mpfr_t Cx, Cy;
  mpfr_init2(Cx, precision);
  mpfr_init2(Cy, precision);
  mpfr_set_str(Cx, argv[1], 10, MPFR_RNDN);
  mpfr_set_str(Cy, argv[2], 10, MPFR_RNDN);

  vector<double> Zxv, Zyv;
  // calculate reference
  {
    Zxv.reserve(iterations);
    Zyv.reserve(iterations);
    mpfr_t Zx, Zy, Zx2, Zy2, Zxy, Zxy2, Z2;
    mpfr_init2(Zx, precision);
    mpfr_init2(Zy, precision);
    mpfr_init2(Zx2, precision);
    mpfr_init2(Zy2, precision);
    mpfr_init2(Zxy, precision);
    mpfr_init2(Zxy2, precision);
    mpfr_init2(Z2, precision);
    // z = 0
    mpfr_set_ui(Zx, 0, MPFR_RNDN);
    mpfr_set_ui(Zy, 0, MPFR_RNDN);
    for (int i = 0; i < iterations; ++i)
    {
      Zxv.push_back(mpfr_get_d(Zx, MPFR_RNDN));
      Zyv.push_back(mpfr_get_d(Zy, MPFR_RNDN));
      mpfr_add(Zxy, Zx, Zy, MPFR_RNDN);
      mpfr_sqr(Zx2, Zx, MPFR_RNDN);
      mpfr_sqr(Zy2, Zy, MPFR_RNDN);
      mpfr_sqr(Zxy2, Zxy, MPFR_RNDN);
      mpfr_add(Z2, Zx2, Zy2, MPFR_RNDN);
      if (mpfr_get_d(Z2, MPFR_RNDZ) > 4) // note rounding mode
      {
        break;
      }
      mpfr_sub(Zx, Zx2, Zy2, MPFR_RNDN);
      mpfr_sub(Zy, Zxy2, Z2, MPFR_RNDN);
      mpfr_add(Zx, Zx, Cx, MPFR_RNDN);
      mpfr_add(Zy, Zy, Cy, MPFR_RNDN);
    }
    mpfr_clear(Z2);
    mpfr_clear(Zxy2);
    mpfr_clear(Zy2);
    mpfr_clear(Zx2);
    mpfr_clear(Zy);
    mpfr_clear(Zx);
  }

  int image[height][width];
  // calculate image
  {
    double *Zxp = &Zxv[0];
    double *Zyp = &Zyv[0];
    int M = Zxv.size() - 1;
    double pixel_spacing_y = 4 / (zoom * height);
    double pixel_spacing_x = pixel_spacing_y / 2;
    #pragma omp parallel for schedule(dynamic, 1)
    for (int y = 0; y < height; ++y)
    {
      for (int x = 0; x < width; ++x)
      {
        double jx = 0.5; // pseudo-random jitter
        double jy = 0.5; // pseudo-random jitter
        double cx = (x + jx - 0.5 * width) * pixel_spacing_x;
        double cy = (y + jy - 0.5 * height) * pixel_spacing_y;
        double zx = 0;
        double zy = 0;
        double zx2 = zx * zx;
        double zy2 = zy * zy;
        image[y][x] = 0;
        int m = 0;
        double Zx = Zxp[m];
        double Zy = Zyp[m];
        for (int n = 0; n < iterations; ++n)
        {
          double zxn = 2 * (zx * Zx - zy * Zy) + zx2 - zy2 + cx;
          double zyn = 2 * (zx * Zy + zy * Zx + zx * zy) + cy;
          zx = zxn;
          zy = zyn;
          m++;
          Zx = Zxp[m];
          Zy = Zyp[m];
          double Zzx = Zx + zx;
          double Zzy = Zy + zy;
          double Zzx2 = Zzx * Zzx;
          double Zzy2 = Zzy * Zzy;
          double Zz2 = Zzx2 + Zzy2;
          if (Zz2 > 4)
          {
            image[y][x] = n + 1;
            break;
          }
          zx2 = zx * zx;
          zy2 = zy * zy;
          double z2 = zx2 + zy2;
          if (Zz2 < z2 || m == M)
          {
            m = 0;
            zx = Zzx;
            zy = Zzy;
            zx2 = Zzx2;
            zy2 = Zzy2;
            Zx = Zxp[m];
            Zy = Zyp[m];
          }
        }
      }
    }
  }

  // output image
  {
    for (int y = 0; y < height; ++y)
    {
      for (int x = 0; x < width; ++x)
      {
        cout << "\33[48;5;" << image[y][x] << "m" << " ";
      }
      cout << endl;
    }
  }

  return 0;
}