/* kohaku1.c - Don Yang (uguu.org)
   Julia set renderer

   10/23/05
                                                                           */
/*@ -realcompare -usedef @*/

#include<stdio.h>
#include<stdlib.h>
#include<time.h>

#define ITERATION_COUNT    26
#define X_STEPS            79
#define Y_STEPS            24
#define 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(double z_r, double z_i, double c_r, double c_i)
{
   double r;
   int i;

   for(i = 0; i < ITERATION_COUNT && z_r * z_r + z_i * z_i <= 4.0; i++)
   {
      r = z_r;
      z_r = z_r * z_r - z_i * z_i + c_r;
      z_i = z_i * r * 2.0 + c_i;
   }
   return i;
}

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

/* Check if a particular set of coordinates is interesting by checking
   for high variety of values in sampled points.                           */
static int interesting(double x0, double x1, double y0, double y1,
                       double c_r, double c_i)
{
   double v[5], a, d, *p;
   int i;

   p = v;
   a  = (*p++ = (double)render_point(x0, y0, c_r, c_i));
   a += (*p++ = (double)render_point(x0, y1, c_r, c_i));
   a += (*p++ = (double)render_point(x1, y0, c_r, c_i));
   a += (*p++ = (double)render_point(x1, y1, c_r, c_i));
   a += (*p++ = (double)render_point((x0 + x1) * 0.5, (y0 + y1) * 0.5,
                                     c_r, c_i));
   a /= 5.0;

   d = 0.0;
   for(i = 0; i < 5; i++)
   {
      --p;
      d += (*p - a) * (*p - a);
   }
   return (int)((d / 5.0) > (ITERATION_COUNT * 4.6));
}

/* Render image */
static void render_image(double x0, double x1, int xsteps,
                         double y0, double y1, int ysteps,
                         double c_r, double c_i)
{
   int x, y;
   for(y = 0; y < ysteps; y++)
   {
      for(x = 0; x < xsteps; x++)
      {
         (void)putchar(iter_char(render_point(
            ((double)x * (x1 - x0)) / (xsteps - 1) + x0,
            ((double)y * (y0 - y1)) / (ysteps - 1) + y1, c_r, c_i)));
      }
      (void)putchar('\n');
   }
}

/* Random numbers */
static double random_range(double x0, double x1)
{
#ifdef __GNUC__
   return ((double)(random() & 0x3fffff) / 0x3fffff) * (x1 - x0) + x0;
#else
   return ((double)(rand() & 0x7fff) / 0x7fff) * (x1 - x0) + x0;
#endif
}

/* Generate coordinates randomly until an interesting one is found,
   then render image.                                                   */
static void render_random_image(void)
{
   double x0, y0, s, c_r, c_i;
   do
   {
      x0 = random_range(-1.0, 1.0);
      y0 = random_range(-1.0, 1.0);
      s = random_range(0.1, 1.0);
      c_r = random_range(-1.0, 1.0);
      c_i = random_range(-1.0, 1.0);
   } while( interesting(x0, x0 + s * ASPECT_RATIO, y0, y0 + s, c_r, c_i)==0 );
   render_image(x0, x0 + s * ASPECT_RATIO, X_STEPS,
                y0, y0 + s, Y_STEPS, c_r, c_i);
}


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