#include <complex>
#include <functional>
#include <vector>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <fstream>

using namespace std;

using Complex = complex<double>;

class Farbe {
  public:
    unsigned char r, g, b;
    Farbe() {}
    Farbe(unsigned char r_, unsigned char g_, unsigned char b_): r(r_), g(g_), b(b_) {}

    friend ostream& operator<<(ostream& stream, const Farbe& farbe) {
      return stream << farbe.r << farbe.g << farbe.b;
    }
};

class Konvergenz {
  private:
    static const vector<Farbe> farbtabelle;
    static constexpr unsigned long itmax = 200;
    static constexpr double eps = 1e-12;

    const function<Complex(Complex)>& f, fderiv;
    double a, b, c, d;
    int M, N;
    double D;

    vector<Complex> nullstellen{};
    vector<int> haeufigkeiten{{0}};

    vector<int> bild;
    ofstream bildstrm;

  public:
    Konvergenz (const function<Complex(Complex)>& f_, const function<Complex(Complex)>& fderiv_,
                double a_, double b_, double c_, double d_, int M_, int N_, double D_,
                string bildstrm_):
      f(f_), fderiv(fderiv_), a(a_), b(b_), c(c_), d(d_), M(M_), N(N_), D(D_), bildstrm(bildstrm_, ios::binary)
    {
      bildstrm << "P6" << " " << M << " " << N << " " << 255 << " ";
    }

    static bool newtonverfahren(Complex z0, Complex& z, function<Complex(Complex)> f, function<Complex(Complex)> fderiv) {
      unsigned long its = itmax;
      Complex z_prev;
      z = z0;
      do {
        z_prev = z;
        z = z - f(z) / fderiv(z);

        its--;
      } while (its > 0 && abs(z - z_prev) > eps * abs(z_prev));

      return its > 0;
    }

    Konvergenz& analysieren() {
      for (int n = N - 1; n >= 0; n--)
        for (int m = 0; m < M; m++) {
          Complex pos{a + m * (b - a) / (M - 1), c + n * (d - c) / (N - 1)};

          Complex nullstelle;
          bool konvergiert = newtonverfahren(pos, nullstelle, f, fderiv);

          vector<int>::size_type i = 0;
          if (konvergiert) {
            for (i = 0; i < nullstellen.size(); i++) {
              if (abs(nullstellen[i] - nullstelle) >= D)
                continue;

              haeufigkeiten[i + 1]++;
              break;
            }

            if (i >= nullstellen.size()) {
              nullstellen.push_back(nullstelle);
              haeufigkeiten.push_back(1);
            }

            i++;
          } else {
            haeufigkeiten[0]++;
          }

          bild.push_back(i);
        }

      for (int nullstelle: bild) {
        int farbe
          = nullstelle
          ? round(((double) (farbtabelle.size() - 2)) * ((double) nullstelle - 1) / ((double) (nullstellen.size() - 1))) + 1
          : 0;
        bildstrm << farbtabelle[farbe];
      }

      return *this;
    }

    friend ostream& operator<<(ostream& out, const Konvergenz& k) {
      out << "N/A(" << k.haeufigkeiten[0] << ")";
      for (vector<int>::size_type i = 0; i < k.nullstellen.size(); i++)
        out << ", " << k.nullstellen[i] << "(" << k.haeufigkeiten[i + 1] << ")";
      return out;
    }
};

const vector<Farbe> Konvergenz::farbtabelle{
  {255u, 255u, 255u},
  {255u, 0u, 0u}, {255u, 255u, 0u}, {0u, 255u, 0u}, {0u, 255u, 255u}, {0u, 0u, 255u}, {255u, 0u, 255u}
};

int main() {
  using namespace std::complex_literals;

  double a, b, c, d, D;
  int M, N;

  cout << "a, b, c, d, M, N, D: ";
  cin >> a >> b >> c >> d >> M >> N >> D;

  cout << setprecision(log10(1.0 / D));

  cout << Konvergenz{
    [] (Complex z) -> Complex { return z * z * z - 1.0; },
    [] (Complex z) -> Complex { return 3.0 * z * z; },
    a, b, c, d, M, N, D,
    "newton_p1.ppm"
  }.analysieren() << endl;

  cout << Konvergenz{
    [] (Complex z) -> Complex { return z * z * z + 3.0 * z * z - 1i; },
    [] (Complex z) -> Complex { return 3.0 * z * z + 6.0 * z; },
    a, b, c, d, M, N, D,
    "newton_p2.ppm"
  }.analysieren() << endl;

  cout << Konvergenz{
    [] (Complex z) -> Complex { return z * z * z * z * z - 2.0 * 1i * z * z * z * z - 13.0 * z * z * z + 14.0 * 1i * z * z + 24.0 * z - 1.0; },
    [] (Complex z) -> Complex { return 5.0 * z * z * z * z - 8.0 * 1i * z * z * z - 36.0 * z * z + 28.0 * 1i * z + 24.0; },
    a, b, c, d, M, N, D,
    "newton_p3.ppm"
  }.analysieren() << endl;

  cout << Konvergenz{
    [] (Complex z) -> Complex { return exp(z) * (z * z * z - 1.0); },
    [] (Complex z) -> Complex { return exp(z) * (z * z * z - 1.0) + exp(z) * 3.0 * z * z; },
    a, b, c, d, M, N, D,
    "newton_f4.ppm"
  }.analysieren() << endl;
}