Фрактал Ляпунова

#include <stdio.h>
#include <math.h>

#define CLAMP(x) (((x) > 255) ? 255 : ((x) < 0) ? 0 : (x))

int seq_length = 12;
int seq[12] = {1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0};
double amin = 3.4;
double amax = 4.0;
double bmin = 2.5;
double bmax = 3.4;
int asize = 1000;
int bsize = 1500;
int nmax = 1000;                /* number of rounds */

/* for color generation; somewhat empirical, in order to match Wickerprint's original colors */
double lambda_min = -2.55;
double lambda_max = 0.3959;

int main(void)
{
  int ai, bi, n, m;
  double a, b, lambda, sum_log_deriv, prod_deriv, r, x, rgb_f[3];
  unsigned char rgb[3];
  FILE *fp;
  
  fp = fopen("lyapunov.ppm", "w");
  fprintf(fp, "P6\n%d %d\n255\n", bsize, asize); /* ppm header */
  
  for (ai = asize-1; ai >= 0; ai--) {
    fprintf(stderr, "\r%d/%d", asize-1-ai, asize);
    for (bi = 0; bi < bsize; bi++) {
      a = amin + (amax-amin)/asize*(ai+0.5);
      b = bmin + (bmax-bmin)/bsize*(bi+0.5);
      x = 0.5;
      /* one round without derivating, so that the value 0.5 is avoided */
      for (m = 0; m < seq_length; m++) {
        r = seq[m] ? b : a;
        x = r*x*(1-x);
      }
      sum_log_deriv = 0;
      for (n = 0; n < nmax; n++) {
        prod_deriv = 1;
        for (m = 0; m < seq_length; m++) {
          r = seq[m] ? b : a;
          /* avoid computing too many logarithms. One every round is acceptable. */
          prod_deriv *= r*(1-2*x); 
          x = r*x*(1-x);
        }
        sum_log_deriv += log(fabs(prod_deriv));
      }
      lambda = sum_log_deriv / (nmax*seq_length);
      
      if (lambda > 0) {
        rgb_f[2] = lambda/lambda_max;
        rgb_f[0] = rgb_f[1] = 0;
      } else {
        rgb_f[0] = 1 - pow(lambda/lambda_min, 2/3.0);
        rgb_f[1] = 1 - pow(lambda/lambda_min, 1/3.0);
        rgb_f[2] = 0;
      }
      rgb[0] = CLAMP(rgb_f[0]*255);
      rgb[1] = CLAMP(rgb_f[1]*255);
      rgb[2] = CLAMP(rgb_f[2]*255);
      fwrite(rgb, 3, 1, fp);
    }
  }
  fprintf(stderr, "\n");
  fclose(fp);
  return 0;
}