/*
 * by Garrett Wollman based on a function posted to Usenet by Ed Kubaitis
 */

/* This program is Copyright (C) 1991, Garrett A. Wollman.  This
 * program may be modified and distributed for any purpose and without
 * fee, provided that this notice remains on all such copies
 * unaltered.  Binary distributions not including this source file are
 * prohibited.  Modified versions must be marked with the name of the
 * modifier and the date of modification.
 *
 * Under no circumstances shall Garrett A. Wollman or the University
 * of Vermont and State Agricultural College be held liable for any
 * damages resulting from the use or misuse of this program, whether
 * the author is aware of such possibility or not.  This program is
 * warranted solely to occupy disk space.
 *
 * I'm sorry for all this legal garbage, but it is sadly necessary
 * these days.
 */

#include <math.h>
#ifndef SUN_BROKEN_STDLIB
#include <stdlib.h>
#endif
#include <stdio.h>
#ifdef sgi
#ifndef NOMULTI
#include <sys/types.h>
#include <sys/prctl.h>
#include <sys/wait.h>

#define MULTIPROC
#define DO_SPROC do_sproc
#define MAXPROCS 8
#endif
#endif

/* other multiprocessing hosts please add your info here */


/*
 * Global variables
 */
#ifdef MULTIPROC
int nprocs = 1;
#endif

char *whoami;

extern int getopt(int,char **,const char *);
extern char *optarg;
extern int optind;
extern int opterr;
extern void perror(const char *);

#ifdef SUN_BROKEN_STDLIB
extern volatile void exit(int);
extern void *malloc(int);
extern int atoi(const char *);
#endif

extern int strlen(const char *);

#ifdef NO_STDIO_PROTOS
extern int fwrite(const char *,int,int,FILE *);
extern int fputc(int,FILE *);
extern int fputs(const char *,FILE *);
extern int pclose(FILE *);
extern void fclose(FILE *);
extern int fprintf(FILE *,const char *,...);
extern int sscanf(char *,const char *,...);
#endif

#ifdef MULTIPROC
#define OPTIONS "r:c:v:m:d:s:g:5pt#:o:l:"
#else
#define OPTIONS "r:c:v:m:d:s:g:5pto:l:"
#endif

int rows = 512,cols = 512;
double amin,amax;
double bmin,bmax;
double aincr,curA;
double bincr;

#define nColors 256            /* don't even try to change this */
struct {
  unsigned char red;
  unsigned char green;
  unsigned char blue;
} colormap[nColors];

int Settle = 600;
int Dwell = 2000;
int *Rvec;
double minLya = -5;
int showpos = 0;
double initGuess = 0.5;

/*
 * Colormap functions
 */

/*
 * set up our initial shades of grey
 */
void init_colormap(void) {
  int i;
  for(i=0; i < nColors; i++)
    colormap[i].red = colormap[i].green = colormap[i].blue = i;
}

void read_colormap(const char *name) {
  FILE *cmap;
  char buf[255];               /* magic number */
  int a,b,c,i;
  
  cmap = fopen(name,"r");
  if(!cmap) {
    perror(name);
    exit(-1);
  }

  i = 0;
  while(!feof(cmap) && i < nColors) {
    fgets(buf,sizeof buf,cmap);
    /*
     * unparseable lines are comments
     * as is anything after the <r g b> value on a single line
     *
     * per Fractint 15.1 standard
     */
    if(sscanf(buf,"%d %d %d",&a,&b,&c) != 3)
      continue;

    colormap[i].red = a;
    colormap[i].green = b;
    colormap[i].blue = c;
    i++;
  }
  fclose(cmap);
}


/*
 * Lyapunov exponent and coloring calculation
 *
 * original function by Ed Kubaitis, used by permission
 */

/*
 * Speedup...
 *
 * Erick B. Wong noticed that $\log_2 d(x_1)+\dots +\log_2 d(x_n)$
 * is, by a principle from high-school algebra, the same as
 * $\log_2 (d(x_1) + \dots + d(x_n))$.  We take advantage of this
 * by unrolling the dwell loop by four (hence the rounding in main())
 * and accumulating by multiplication before taking the log.  Since
 * log() is extremely expensive, this should lead to a considerable
 * speedup, while still allowing complete flexibility in selecting
 * dwell values.
 *
 * Presumably, one might unroll this loop even further, with a 
 * concomitant increase in textual complexity.
 */

int lya(double a,double b) {
  double t = 0;
  double r, lya;
  int d;
  double acc;
  double x;

  x = initGuess;               /* initialize x */

  for(d=0; d < Settle; d++) {
    r = (Rvec[d]) ? b : a;
    x = r*x*(1-x);
  }

#if 1
  for(d=0; d <= Dwell; d+= 4) {
    r = (Rvec[d]) ? b : a;
    x = r * x * (1-x);
    acc = fabs(r - 2*r*x);

    r = (Rvec[d+1]) ? b : a;
    x = r * x * (1-x);
    acc *= fabs(r - 2*r*x);

    r = (Rvec[d+2]) ? b : a;
    x = r * x * (1-x);
    acc *= fabs(r - 2*r*x);

    r = (Rvec[d+3]) ? b : a;
    x = r * x * (1-x);
    t += log(acc * fabs(r - 2*r*x));

    if(fabs(x-0.5) < 1e-20) {
      d += 3;
      break;
    }
  }
#else
  for(d = 0; d <= Dwell; d += 4) {
    r = (Rvec[d]) ? b : a;
    x = r * x * (1-x);
    t += log(fabs(r-2*r*x));

    r = (Rvec[d+1]) ? b : a;
    x = r * x * (1-x);
    t += log(fabs(r-2*r*x));

    r = (Rvec[d+2]) ? b : a;
    x = r * x * (1-x);
    t += log(fabs(r-2*r*x));

    r = (Rvec[d+3]) ? b : a;
    x = r * x * (1-x);
    t += log(fabs(r-2*r*x));

    if(fabs(x-0.5) < 1e-20) {
      d += 3;
      break;
    }
  }
#endif

  lya = (t * 1.442695041)/d;

  if(showpos)                  
    lya *= -1;
  if(lya < minLya)             /* sanity check! */
    lya = minLya;

  if(lya < 0) {
    return (int)(lya / minLya * (double)nColors);
  } else {
    return 0;
  }
}


/*
 * Usage message
 */
#ifdef __GNUC__
volatile
#endif
void usage(void) {
  fprintf(stderr,
         "%s: usage:\n\n"
         "%s [-r rows] [-c cols] [-v program] [-p] [-t]\n\t"
#ifdef MULTIPROC
         "[-# processors] "
#endif
         "[-o file] [-l colormap]\n\t"
         "[-m minLya] [-d Dwell] [-s Settle] [-g initGuess]\n\t"
         "[-5] amin amax bmin bmax bits\n",whoami,whoami);
  exit(-1);
}


/*
 * SGI-specific multiprocessing code
 */
#ifdef sgi
#ifndef NOMULTI
struct info {
  int *row;
  double curA;
};

void real_doit(double curA,int *row) {
  int onCol;
  double curB;

  for(onCol = cols - 1,curB = bmax; onCol; onCol--,curB -= bincr) {
    row[onCol] = lya(curA,curB);
  }
}

void doit(void *info) {
  struct info *myinfo = info;

  real_doit(myinfo->curA,myinfo->row);
  exit(0);
}

void do_sproc(int *manyRows[]) {
  int pids[MAXPROCS];
  struct info infos[MAXPROCS];
  int i;
  double aval = curA;

  for(i=0; i < nprocs - 1; i++,aval -= aincr) {
    infos[i].row = manyRows[i];
    infos[i].curA = aval;
    pids[i] = sproc(doit,PR_SADDR,(void *)&infos[i]);
  }
  real_doit(aval,manyRows[i]);

  /*
   * We know we started (nprocs - 1) kids, so we wait for (nprocs - 1) kids
   * to die.
   */
  for(i=0; i < nprocs - 1; i++) {
    (void)wait((int *)0);
  }
}
#endif
#endif


int main(int argc,char **argv) {
  int c,i,j;
  char *viewprogram = NULL;
  char *outname = NULL;
  FILE *outfile;

  char *bits;

  double curB;
  int onRow,onCol;
#ifdef MULTIPROC
  int *manyRows[MAXPROCS];
#endif
  int *oneRow;

  int usepgm = 0;
  int outtext = 0;

  init_colormap();
  opterr = 0;

  whoami = *argv;
  while((c = getopt(argc,argv,OPTIONS)) != EOF) {
    switch(c) {
    case 'r':
      rows = atoi(optarg);
      if(!rows)
       usage();
      break;

    case 'c':
      cols = atoi(optarg);
      if(!cols)
       usage();
      break;

#ifdef MULTIPROC
    case '#':
      nprocs = atoi(optarg);
      if(!nprocs)
       usage();
      break;
#endif

    case 'o':
      outname = optarg;
      viewprogram = NULL;
      break;

    case 'v':
      viewprogram = optarg;
      outname = NULL;
      break;

    case 'l':
      read_colormap(optarg);
      break;

    case 'm':
      if(!sscanf(optarg,"%lf",&minLya))
       usage();
      break;

    case 'd':
      Dwell = atoi(optarg);
      if(!Dwell)
       usage();
      /*
       * By forcing the Dwell to be a multiple of four, we can use the
       * speedup noted in lya().  We *could* just do this silently,
       * by if it might make a difference to someone...
       */
      if(Dwell % 4) {
       fprintf(stderr,"Dwell value of %d being rounded off to %d.\n",
               Dwell,
               Dwell += 4 - (Dwell % 4));
      }
      break;

    case 's':
      Settle = atoi(optarg);
      if(!Settle)              /* for settle close to 0, use 1 */
       usage();
      break;

    case 'g':
      if(!sscanf(optarg,"%lf",&initGuess))
       usage();
      break;

    case '5':
      usepgm = 1;
      break;

    case 'p':
      showpos = 1;
      break;

    case 't':
      outtext = 1;
      break;

    default:                   /* includes '?' (and '#' when !MULTIPROC) */
      usage();
    }
  }

  if(!argv[optind])
    usage();

  /*
   * Note to language nit-pickers... The code below is legal because
   * the || operator is a sequence point.
   */
  if(!sscanf(argv[optind++],"%lf",&amin) || !argv[optind] ||
     !sscanf(argv[optind++],"%lf",&amax) || !argv[optind] ||
     !sscanf(argv[optind++],"%lf",&bmin) || !argv[optind] ||
     !sscanf(argv[optind++],"%lf",&bmax) || !argv[optind])
    usage();

  bits = argv[optind];

  if(Dwell < Settle)
    Dwell = Settle;            /* a bit of sanity */

  j = strlen(bits);

  Rvec = malloc((Dwell+1) * sizeof *Rvec);
  if(!Rvec) {
    fputs("Cannot allocate space for Markus vector!\n",stderr);
    exit(-1);
  }

  for(i=0; i <= Dwell; i++) {
    if(bits[i % j] == 'a')
      Rvec[i] = 0;
    else if(bits[i % j] == 'b')
      Rvec[i] = 1;
    else {
      fprintf(stderr,"Invalid character \\%o ('%c') in bit string\n",
             (int)(unsigned char)bits[i % j], bits[i % j]);
      exit(-1);
    }
  }

#ifdef MULTIPROC
  for(i=0; i < nprocs; i++) {
    manyRows[i] = malloc(cols * sizeof *manyRows[i]);
    if(!manyRows[i]) {
      fputs("Cannot allocate space for row buffers!\n",stderr);
      exit(-1);
    }
  }
#else 
  oneRow = malloc(cols * sizeof *oneRow);
  if(!oneRow) {
    fputs("Cannot allocate space for single-row buffer!\n",stderr);
    exit(-1);
  }
#endif

  aincr = (amax - amin) / rows;
  bincr = (bmax - bmin) / cols;

  if(viewprogram) {
    outfile = popen(viewprogram,"w");
    if(!outfile) {
      fprintf(stderr,"Couldn't popen(\"%s\",\"w\")!\n",viewprogram);
      exit(-1);
    }
  } else if(outname) {
    outfile = fopen(outname,"w");
    if(!outfile) {
      perror(outname);
      exit(-1);
    }
  } else {
    outfile = stdout;
  }

  /*
   * Print P[PG]M header...
   */
  fprintf(outfile,"P%d %d %d\n",(usepgm?2:3) + (outtext?0:3),cols,rows);
  fprintf(outfile,"%d\n",nColors-1);

#ifndef MULTIPROC
  for(onRow = 0, curA = amax; onRow < rows; curA -= aincr, onRow++) {
    for(onCol = cols - 1, curB = bmax; onCol; curB -= bincr, onCol--) {
      oneRow[onCol] = lya(curA,curB);
    }
    
#else
  for(onRow = 0, curA = amax; onRow < rows;
      curA -= aincr*nprocs, onRow += nprocs) {
    DO_SPROC(manyRows);

    for(j = 0; j < nprocs; j++) {
      oneRow = manyRows[j];
#endif
    for(i=0; i < cols; i++) {
      unsigned char c;
      c = (unsigned char)oneRow[i];
      if(outtext) {
       if(usepgm)
         fprintf(outfile,"%d ",c);
       else
         fprintf(outfile,"%d %d %d ",colormap[c].red,colormap[c].green,
                 colormap[c].blue);
      } else if(usepgm) {
       fwrite((char *)&c,sizeof c,1,outfile);
      } else {
       fwrite((char *)&colormap[c],sizeof colormap[0],1,outfile);
      }
    }
#ifdef MULTIPROC
    }
#endif
    if(!(onRow % 16)) {
      if(outtext)
       fputc('\n',outfile);
      fprintf(stderr,"Done row %d\n",onRow);
    }
  }

  if(viewprogram)
    pclose(outfile);
  else if(outname)
    fclose(outfile);

  return 0;
}