//-----------------------------------------------------
// Produces Mandelbrot plot in the range [-2:1]x[-1:1]
// It uses formula z = z*z + z0 iteratively until
// asb(z) gets bigger than 2 
// (deciding that z0 is not in mandelbrot set)
//-----------------------------------------------------
#include <iostream>
#include <complex>
using namespace std;

double Mandelb(complex<double>& z0, int max_iterations=1000) // default value
{
  complex<double> z=0; // initial value of z
  for (int i=0; i<max_iterations; i++){
    if (norm(z)>4.) return i;
    z = z*z + z0;     // if |z|>2 the point is not part of mandelbrot set
  }
  return max_iterations*1e3;
}

double Mandelb_optimized(double x0, double y0, int max_iterations=1000)
{
  double x = 0, y = 0; // current coordiates are set to starting coordinate
  double x2 = x*x, y2 = y*y;
  for (int i=0; i<max_iterations; i++){
    if (x2+y2>4) return i; // if |z|>2 the point is not part of mandelbrot set
    y = 2*x*y + y0;        // z = z*z+z0
    x = x2 - y2 + x0;
    x2 = x*x;              // x^2 is precomputed for speed
    y2 = y*y;              // y^2
  }
  return max_iterations*1e3;
}


int main()   // the main program should always return int: 0 - if no error occured,
             //                                      non-zero  if an error occured.
{
  int Nx=400;
  int Ny=400;
  for (int i=0; i<Nx; i++){
    for (int j=0; j<Ny; j++){
      double x = -2. + 3*i/(Nx-1.);
      double y = -1. + 2*j/(Ny-1.);
      complex<double> z0(x, y);
      //cout<<x<<" "<<y<<" "<<1/Mandelb(z0)<<endl;
      cout<<x<<" "<<y<<" "<<1/Mandelb_optimized(x,y)<<endl;
    }
  }
  return 0; // no error
}