#include <stdint.h> // uint32_t, uint64_t
#include <stdlib.h> // abs
#include <ctime> // time
#include <cmath> // INFINITY, powf
#include <algorithm> // min, max
#include <iostream> // cout
#include <iomanip> // setprecision
#include <fstream> // ofstream
using namespace std;

// Dla porządku wszystko umieszczamy wewnątrz osobnej przestrzeni nazw.
namespace smok {

// Punkt w dwóch wymiarach (x,y).
  class Punkt2D {
    public:
      Punkt2D(double x, double y): x(x), y(y) {}
// Kilka przydatnych operatorów.
// Do rozszerzania w miarę potrzeb.
// Mnożenie przez stałą (skalowanie)
      Punkt2D operator *(double skala) { return Punkt2D(x*skala, y*skala); }
// Dodawanie / odejmowanie (przesunięcie)
      Punkt2D operator -(const Punkt2D& p) { return Punkt2D(x - p.x, y - p.y); }
      Punkt2D operator +(const Punkt2D& p) { return Punkt2D(x + p.x, y + p.y); }

      double x;
      double y;
  };

// Prosty generator liczb pseudolosowych typu LCG
// (Linear Congruential Generator).
// Stałe na podstawie książki Numerical Recipes.
// Przydatny tam, gdzie chcemy mieć powtarzalność
// wyników. W tym programie - niespecjalnie.
  class RndLCG {
    protected:
      static const uint32_t a = 1664525;
      static const uint32_t c = 1013904223;
      uint32_t x;
    public:
      RndLCG(uint32_t seed): x(seed) {}
// Poniższa definicja operatora () czyni tę klasę "generatorem",
// tzn. wywołując obiekt klasy tak, jak byśmy wywoływali funkcję,
// dostajemy kolejne wygenerowane wartości.
// Możemy podać taką klasę jako argument dla std::generate
// (http://www.cplusplus.com/reference/algorithm/generate/)
// by łatwo wypełnić wygenerowanymi wartościami tablicę, czy
// wektor, czy listę, czy inny kontener.
      double operator()() {

// Poniższy wzór zawiera niejawną operację modulo 2³²,
// ponieważ operujemy na typie uint32 mającym dokładnie
// 32 bity długości.
        x = x * a + c;
        return x / 4294967296.0;
      }
  };

  class Smok {
    protected:
      Punkt2D p;
      RndLCG random;
    public:
// Generator pseudolosowy inicjujemy aktualnym czasem.
// Funkcja time() zwraca ilość sekund od 1 stycznia 1970.
// Jest to więc bardzo kiepska i łatwa do zgadnięcia
// wartość, ale dla tych zastosowań wystarczy.
//
// Inicjalizacja rand w tej klasie MUSI być zrobiona w
// ramach listy inicjalizacji, a nie wewnątrz w kodzie
// konstruktora, ponieważ RndLCG nie ma zdefiniowanego
// bezparametrowego konstruktora, i nie wiem czy dobrym
// pomysłem byłoby taki konstruktor zrobić (szczególnie
// tak tandetny jak użycie time(NULL)). Generowanie
// wartości seed to chyba dobre zadanie dla osobnej
// funkcji / klasy.
      Smok(double x = 1, double y = 1): random(time(NULL)), p(x, y) {}
// Generator wartości (analogicznie jak w RndLCG).
      Punkt2D operator()() {
        if (random() >= 0.5) {
          p.x = -0.4*p.x - 1;
          p.y = -0.4*p.y + 0.1;
        } else {
          double tmp = 0.76*p.x - 0.4*p.y;
          p.y = 0.4*p.x + 0.76*p.y;
          p.x = tmp;
        }
        return p;
      }

  };

// Prostokąt jest zawsze 2D :)
  class Prostokat {
    public:
// Definiowany przez 2 narożne punkty
      Punkt2D p1;
      Punkt2D p2;

      Prostokat(const Punkt2D& p1, const Punkt2D& p2): p1(p1), p2(p2) {}
// Szerokość i wysokość prostokąta zawsze dodatnia (abs)
      double szer() { return abs(p2.x - p1.x); }
      double wys() { return abs(p2.y - p1.y); }
  };

// Prostokąt ograniczający obiekty wewnątrz niego,
// tzw. Bounding Rectangle czy też Bounding Box.
// Dziedziczy ze zwykłego prostokąta, ale ma funkcję
// do dodawania punktów do niego.
  class ProstOgran: public Prostokat {
    private:
// long na wypadek, gdyby tych punktów było naprawdę dużo
      long ilosc;
      Punkt2D sm;
    public:
      ProstOgran(): Prostokat(Punkt2D(INFINITY, INFINITY), Punkt2D(-INFINITY, -INFINITY)),
        sm(Punkt2D(0, 0)), ilosc(0) {}
// Dodawanie punktów powiększa prostokąt.
      void dodajPunkt(const Punkt2D& p) {
        p1.x = min(p1.x, p.x);
        p1.y = min(p1.y, p.y);
        p2.x = max(p2.x, p.x);
        p2.y = max(p2.y, p.y);
// Średnie x i y (środek masy)
        ilosc++;
        sm.x += (p.x - sm.x) / ilosc;
        sm.y += (p.y - sm.y) / ilosc;
      }
      Punkt2D srodek_masy() {
        return sm;
      }
// Można rozszerzyć klasę o kolejne metody do dodawania
// innych obiektów.
//    void dodajProstokat(const Prostokat& pr) {
// [...]
// itp. itd.
  };

  class Obraz {
    protected:
      int s;
      int w;
// Nasz obraz składa się z 64-bitowych wartości,
// co pozwala z każdej wartości zrobić kwadrat
// o wielkości 8 na 8 pikseli.
      uint64_t* bitmap;
    public:
      Obraz(int s, int w): s(s), w(w),
// Alokacja pamięci na obraz.
// Forma new ...() inicjalizuje pamięć w przeciwieństwie
// do formy bez nawiasów.
        bitmap( new uint64_t[s*w]() ) {}
// Ponieważ w konstruktorze przydzielamy pamięć
// na obraz, w destruktorze musimy ją zwolnić.
// Używamy delete[] bo przydzielaliśmy przez new ...[].
      ~Obraz() { delete[] bitmap; }

      void dodajPunkt(const Punkt2D& p) {
// Jeśli punkt wychodzi poza obraz to nie rysujemy
        if (p.x < 0 || p.x >= s) return;
        if (p.y < 0 || p.y >= w) return;

// Wartości całkowite współrzędnych x i y określają
// indeks w tablicy.
        int indeks = (int)p.y * s + (int)p.x;
// Wartości ułamkowe określają numer bitu od 0 do 63
        int x_s = 7 & int(p.x * 8);
        int y_s = 7 & int(p.y * 8);
        int bit = y_s * 8 + x_s;
// "Zapalamy" wybrany bit.
// Rzutowanie na uint64_t jest konieczne!
// Inaczej kompilator weźmie zwykłego inta,
// w którym nasza wartość się nie zmieści.
        bitmap[indeks] |= (uint64_t)1 << bit;
      }

      int szer() { return s; }
      int wys() { return w; }

// Obliczenie ilości zapalonych bitów.
      int ilosc_bitow(uint64_t x) {
        int ile;
        for (ile = 0; x; ile++) x &= x-1;
        return ile;
      }

// Zwraca jasność punktu od 0 do 1 (włącznie)
// pod współrzędnymi (x,y).
      float jasnosc(int x, int y) {
        int indeks = y * s + x;
        return ilosc_bitow(bitmap[indeks]) / 64.0f; // f na końcu oznacza float
      }
  };

// Konwersja z jasności liniowej od 0 do 1 na
// sRGB, ponieważ tego oczekują monitory.
// Bez tego obraz będzie zbyt ciemny.
  float to_sRGB(float j) {
    if (j <= 0.0031308f) {
      return 12.92f * j;
    } else {
      return 1.055f * powf(j, 0.4166667f) - 0.055f;
    }
  }

  int main() {
    long iter_oblicz = 5000;
    long iter_obraz = 50000000;
// Maksymalna szerokość i wysokość obrazka.
// Jeden z wymiarów może wyjść mniejszy.
    int obr_szer = 800;
    int obr_wys = 600;
// Dodatkowe pomniejszenie jako margines błędu w przypadku,
// gdy do obliczeń bierzemy mniej iteracji.
    double skala_dod = 0.90;

    Smok smok;

// Pierwsze 100 wyników odrzucamy
    for(int i = 0; i < 100; i++) smok();

// Robimy kopię smoka do oszacowania wymiarów
    Smok smok_obl = smok;

    ProstOgran pr_o;
    for(long i = 0; i < iter_oblicz; i++) pr_o.dodajPunkt(smok_obl());

// Wypisanie obliczonych wartości.
// setprecision(2) nam ładnie zaokrągla do 2 miejsc po przecinku
    cout << setprecision(2);
    cout << "Min x: " << pr_o.p1.x << "; Min y: " << pr_o.p1.y << endl;
    cout << "Max x: " << pr_o.p2.x << "; Max y: " << pr_o.p2.y << endl;
    cout << "Szerokość: " << pr_o.szer() << "; Wysokość: " << pr_o.wys() << endl;
    cout << "Środek masy: " << pr_o.srodek_masy().x << ", " << pr_o.srodek_masy().y << endl;

// Obliczenie skali w poziomie i w pionie
    double skala_x = obr_szer / pr_o.szer();
    double skala_y = obr_wys / pr_o.wys();
// Ogólna skala to minimum z tych dwóch wartości
    double skala = min(skala_x, skala_y);
// Ostateczne wymiary obrazu na podstawie skali
    Obraz obr(pr_o.szer() * skala, pr_o.wys() * skala);
// Wyśrodkowanie po dodatkowym pomniejszeniu
    Punkt2D wysrodkowanie = Punkt2D(obr.szer(), obr.wys()) * (1 - skala_dod) * 0.5;
// Rysowanie smoka.
    for(long i = 0; i < iter_obraz; i++) obr.dodajPunkt((smok() - pr_o.p1) * skala * skala_dod + wysrodkowanie);
// Utworzenie pliku z obrazkiem.
    ofstream plik;
    plik.open("smok.pgm");
// P2 oznacza obrazek w odcieniach szarości.
    plik << "P2" << endl;
    plik << obr.szer() << " " << obr.wys() << endl;
// 256 poziomów szarości (od 0 do 255 włącznie).
    plik << 255 << endl;
// Konieczne dodanie "fixed" by nam nie zapisywał liczb
// w postaci wykładniczej!
    plik << fixed << setprecision(0);
    for(int y = 0; y < obr.wys(); y++) {
      for(int x = 0; x < obr.szer(); x++) {
// Konwertujemy do sRGB i zakresu 0-255
        plik << to_sRGB(obr.jasnosc(x, y)) * 255 << " ";
      }
      plik << endl;
    }
    plik.close();
    return 0;
  }

// koniec namespace
}

int main() {
  return smok::main();
}