/* kohaku0.cc - Don Yang (uguu.org)
   Julia set renderer

   Done functionally and not efficiently.

   10/13/05
                                                                           */

#include<stdio.h>
#include<time.h>
#include<complex>
#include<vector>
using namespace std;

// Constants
static const int iteration_count = 26;
static const int x_steps = 79;
static const int y_steps = 24;
static const double aspect_ratio = 1.5;


// Render a point by exponentiating it repeatedly.  Returns number of
// iterations it takes for the number to diverge, clamped by
// iteration_count.
static int render_point(const complex<double> &z,
                        const complex<double> &c, int i = 0)
{
   return (i > iteration_count || abs(z) > 2.0)
      ? i : render_point(z * z + c, c, i + 1);
}


// Convert iteration value to character.
static char iter_char(int i)
{
   return (i >= 26 || i < 0)
      ? ' ' : static_cast<char>(static_cast<int>('A') + i);
}


// Check if a particular set of coordinates is interesting by checking
// for high variety of values in sampled points.
static double average(const vector<double> &series)
{
   double sum = 0.0;
   for(vector<double>::const_iterator i = series.begin();
       i != series.end(); ++i)
      sum += *i;
   return sum / series.size();
}
static double stddev(const vector<double> &series)
{
   double a = average(series), sum_sq = 0.0;
   for(vector<double>::const_iterator i = series.begin();
       i != series.end(); ++i)
      sum_sq += (*i - a) * (*i - a);
   return sum_sq / series.size();
}
static bool interesting(double x0, double x1, double y0, double y1,
                        const complex<double> &c)
{
   vector<double> series;
   series.push_back(static_cast<double>(
      render_point(complex<double>(x0, y0), c)));
   series.push_back(static_cast<double>(
      render_point(complex<double>(x0, y1), c)));
   series.push_back(static_cast<double>(
      render_point(complex<double>(x1, y0), c)));
   series.push_back(static_cast<double>(
      render_point(complex<double>(x1, y1), c)));
   series.push_back(static_cast<double>(
      render_point(complex<double>((x0 + x1) * 0.5, (y0 + y1) * 0.5), c)));
   return stddev(series) > iteration_count * 4.6;
}


// Generate a list of evenly spaced values
static void generate_series(double x0, double x1, int steps,
                            vector<double> *series)
{
   for(int i = 0; i < steps; i++)
      series->push_back(x0 + i * (x1 - x0) / static_cast<double>(steps - 1));
}


// Render image
static void render_scanline(double x0, double x1, int xsteps,
                            double y, const complex<double> &c)
{
   vector<double> series;
   generate_series(x0, x1, xsteps, &series);
   for(vector<double>::const_iterator x = series.begin();
       x != series.end(); ++x)
      (void)putchar(iter_char(render_point(complex<double>(*x, y), c)));
   (void)putchar('\n');
}
static void render_image(double x0, double x1, int xsteps,
                         double y0, double y1, int ysteps,
                         const complex<double> &c)
{
   vector<double> series;
   generate_series(y1, y0, ysteps, &series);
   for(vector<double>::const_iterator y = series.begin();
       y != series.end(); ++y)
      render_scanline(x0, x1, xsteps, *y, c);
}


// Random numbers
static double random_range(double x0, double x1)
{
#ifdef __GNUC__
   return (static_cast<double>(random() & 0x3fffff) / 0x3fffff) * (x1-x0) + x0;
#else
   return (static_cast<double>(rand() & 0x7fff) / 0x7fff) * (x1 - x0) + x0;
#endif
}
static void random_complex(complex<double> *c)
{
   *c = complex<double>(random_range(-1.0, 1.0), random_range(-1.0, 1.0));
}


// Generate coordiates repeatedly until an interesting one is found,
// then render image.
static void render_random_image(void)
{
   double x0, y0, s;
   complex<double> c;
   do
   {
      x0 = random_range(-1.0, 1.0);
      y0 = random_range(-1.0, 1.0);
      s = random_range(0.1, 1.0);
      random_complex(&c);
   } while( !interesting(x0, x0 + s * aspect_ratio, y0, y0 + s, c) );
   render_image(x0, x0 + s * aspect_ratio, x_steps, y0, y0 + s, y_steps, c);
}

int main(void)
{
#ifdef __GNUC__
   srandom(time(NULL));
#else
   srand(time(NULL));
#endif
   render_random_image();
   return 0;
}