home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume23 / lyapunov / part01 / lyapunov.c next >
Encoding:
C/C++ Source or Header  |  1991-09-27  |  10.9 KB  |  548 lines
  1. /*
  2.  * by Garrett Wollman based on a function posted to Usenet by Ed Kubaitis
  3.  */
  4.  
  5. /* This program is Copyright (C) 1991, Garrett A. Wollman.  This
  6.  * program may be modified and distributed for any purpose and without
  7.  * fee, provided that this notice remains on all such copies
  8.  * unaltered.  Binary distributions not including this source file are
  9.  * prohibited.  Modified versions must be marked with the name of the
  10.  * modifier and the date of modification.
  11.  *
  12.  * Under no circumstances shall Garrett A. Wollman or the University
  13.  * of Vermont and State Agricultural College be held liable for any
  14.  * damages resulting from the use or misuse of this program, whether
  15.  * the author is aware of such possibility or not.  This program is
  16.  * warranted solely to occupy disk space.
  17.  *
  18.  * I'm sorry for all this legal garbage, but it is sadly necessary
  19.  * these days.
  20.  */
  21.  
  22. #include <math.h>
  23. #ifndef SUN_BROKEN_STDLIB
  24. #include <stdlib.h>
  25. #endif
  26. #include <stdio.h>
  27. #ifdef sgi
  28. #ifndef NOMULTI
  29. #include <sys/types.h>
  30. #include <sys/prctl.h>
  31. #include <sys/wait.h>
  32.  
  33. #define MULTIPROC
  34. #define DO_SPROC do_sproc
  35. #define MAXPROCS 8
  36. #endif
  37. #endif
  38.  
  39. /* other multiprocessing hosts please add your info here */
  40.  
  41.  
  42. /*
  43.  * Global variables
  44.  */
  45. #ifdef MULTIPROC
  46. int nprocs = 1;
  47. #endif
  48.  
  49. char *whoami;
  50.  
  51. extern int getopt(int,char **,const char *);
  52. extern char *optarg;
  53. extern int optind;
  54. extern int opterr;
  55. extern void perror(const char *);
  56.  
  57. #ifdef SUN_BROKEN_STDLIB
  58. extern volatile void exit(int);
  59. extern void *malloc(int);
  60. extern int atoi(const char *);
  61. #endif
  62.  
  63. extern int strlen(const char *);
  64.  
  65. #ifdef NO_STDIO_PROTOS
  66. extern int fwrite(const char *,int,int,FILE *);
  67. extern int fputc(int,FILE *);
  68. extern int fputs(const char *,FILE *);
  69. extern int pclose(FILE *);
  70. extern void fclose(FILE *);
  71. extern int fprintf(FILE *,const char *,...);
  72. extern int sscanf(char *,const char *,...);
  73. #endif
  74.  
  75. #ifdef MULTIPROC
  76. #define OPTIONS "r:c:v:m:d:s:g:5pt#:o:l:"
  77. #else
  78. #define OPTIONS "r:c:v:m:d:s:g:5pto:l:"
  79. #endif
  80.  
  81. int rows = 512,cols = 512;
  82. double amin,amax;
  83. double bmin,bmax;
  84. double aincr,curA;
  85. double bincr;
  86.  
  87. #define nColors 256        /* don't even try to change this */
  88. struct {
  89.   unsigned char red;
  90.   unsigned char green;
  91.   unsigned char blue;
  92. } colormap[nColors];
  93.  
  94. int Settle = 600;
  95. int Dwell = 2000;
  96. int *Rvec;
  97. double minLya = -5;
  98. int showpos = 0;
  99. double initGuess = 0.5;
  100.  
  101. /*
  102.  * Colormap functions
  103.  */
  104.  
  105. /*
  106.  * set up our initial shades of grey
  107.  */
  108. void init_colormap(void) {
  109.   int i;
  110.   for(i=0; i < nColors; i++)
  111.     colormap[i].red = colormap[i].green = colormap[i].blue = i;
  112. }
  113.  
  114. void read_colormap(const char *name) {
  115.   FILE *cmap;
  116.   char buf[255];        /* magic number */
  117.   int a,b,c,i;
  118.   
  119.   cmap = fopen(name,"r");
  120.   if(!cmap) {
  121.     perror(name);
  122.     exit(-1);
  123.   }
  124.  
  125.   i = 0;
  126.   while(!feof(cmap) && i < nColors) {
  127.     fgets(buf,sizeof buf,cmap);
  128.     /*
  129.      * unparseable lines are comments
  130.      * as is anything after the <r g b> value on a single line
  131.      *
  132.      * per Fractint 15.1 standard
  133.      */
  134.     if(sscanf(buf,"%d %d %d",&a,&b,&c) != 3)
  135.       continue;
  136.  
  137.     colormap[i].red = a;
  138.     colormap[i].green = b;
  139.     colormap[i].blue = c;
  140.     i++;
  141.   }
  142.   fclose(cmap);
  143. }
  144.  
  145.  
  146. /*
  147.  * Lyapunov exponent and coloring calculation
  148.  *
  149.  * original function by Ed Kubaitis, used by permission
  150.  */
  151.  
  152. /*
  153.  * Speedup...
  154.  *
  155.  * Erick B. Wong noticed that $\log_2 d(x_1)+\dots +\log_2 d(x_n)$
  156.  * is, by a principle from high-school algebra, the same as
  157.  * $\log_2 (d(x_1) + \dots + d(x_n))$.  We take advantage of this
  158.  * by unrolling the dwell loop by four (hence the rounding in main())
  159.  * and accumulating by multiplication before taking the log.  Since
  160.  * log() is extremely expensive, this should lead to a considerable
  161.  * speedup, while still allowing complete flexibility in selecting
  162.  * dwell values.
  163.  *
  164.  * Presumably, one might unroll this loop even further, with a 
  165.  * concomitant increase in textual complexity.
  166.  */
  167.  
  168. int lya(double a,double b) {
  169.   double t = 0;
  170.   double r, lya;
  171.   int d;
  172.   double acc;
  173.   double x;
  174.  
  175.   x = initGuess;        /* initialize x */
  176.  
  177.   for(d=0; d < Settle; d++) {
  178.     r = (Rvec[d]) ? b : a;
  179.     x = r*x*(1-x);
  180.   }
  181.  
  182. #if 1
  183.   for(d=0; d <= Dwell; d+= 4) {
  184.     r = (Rvec[d]) ? b : a;
  185.     x = r * x * (1-x);
  186.     acc = fabs(r - 2*r*x);
  187.  
  188.     r = (Rvec[d+1]) ? b : a;
  189.     x = r * x * (1-x);
  190.     acc *= fabs(r - 2*r*x);
  191.  
  192.     r = (Rvec[d+2]) ? b : a;
  193.     x = r * x * (1-x);
  194.     acc *= fabs(r - 2*r*x);
  195.  
  196.     r = (Rvec[d+3]) ? b : a;
  197.     x = r * x * (1-x);
  198.     t += log(acc * fabs(r - 2*r*x));
  199.  
  200.     if(fabs(x-0.5) < 1e-20) {
  201.       d += 3;
  202.       break;
  203.     }
  204.   }
  205. #else
  206.   for(d = 0; d <= Dwell; d += 4) {
  207.     r = (Rvec[d]) ? b : a;
  208.     x = r * x * (1-x);
  209.     t += log(fabs(r-2*r*x));
  210.  
  211.     r = (Rvec[d+1]) ? b : a;
  212.     x = r * x * (1-x);
  213.     t += log(fabs(r-2*r*x));
  214.  
  215.     r = (Rvec[d+2]) ? b : a;
  216.     x = r * x * (1-x);
  217.     t += log(fabs(r-2*r*x));
  218.  
  219.     r = (Rvec[d+3]) ? b : a;
  220.     x = r * x * (1-x);
  221.     t += log(fabs(r-2*r*x));
  222.  
  223.     if(fabs(x-0.5) < 1e-20) {
  224.       d += 3;
  225.       break;
  226.     }
  227.   }
  228. #endif
  229.  
  230.   lya = (t * 1.442695041)/d;
  231.  
  232.   if(showpos)            
  233.     lya *= -1;
  234.   if(lya < minLya)        /* sanity check! */
  235.     lya = minLya;
  236.  
  237.   if(lya < 0) {
  238.     return (int)(lya / minLya * (double)nColors);
  239.   } else {
  240.     return 0;
  241.   }
  242. }
  243.  
  244.  
  245. /*
  246.  * Usage message
  247.  */
  248. #ifdef __GNUC__
  249. volatile
  250. #endif
  251. void usage(void) {
  252.   fprintf(stderr,
  253.       "%s: usage:\n\n"
  254.       "%s [-r rows] [-c cols] [-v program] [-p] [-t]\n\t"
  255. #ifdef MULTIPROC
  256.       "[-# processors] "
  257. #endif
  258.       "[-o file] [-l colormap]\n\t"
  259.       "[-m minLya] [-d Dwell] [-s Settle] [-g initGuess]\n\t"
  260.       "[-5] amin amax bmin bmax bits\n",whoami,whoami);
  261.   exit(-1);
  262. }
  263.  
  264.  
  265. /*
  266.  * SGI-specific multiprocessing code
  267.  */
  268. #ifdef sgi
  269. #ifndef NOMULTI
  270. struct info {
  271.   int *row;
  272.   double curA;
  273. };
  274.  
  275. void real_doit(double curA,int *row) {
  276.   int onCol;
  277.   double curB;
  278.  
  279.   for(onCol = cols - 1,curB = bmax; onCol; onCol--,curB -= bincr) {
  280.     row[onCol] = lya(curA,curB);
  281.   }
  282. }
  283.  
  284. void doit(void *info) {
  285.   struct info *myinfo = info;
  286.  
  287.   real_doit(myinfo->curA,myinfo->row);
  288.   exit(0);
  289. }
  290.  
  291. void do_sproc(int *manyRows[]) {
  292.   int pids[MAXPROCS];
  293.   struct info infos[MAXPROCS];
  294.   int i;
  295.   double aval = curA;
  296.  
  297.   for(i=0; i < nprocs - 1; i++,aval -= aincr) {
  298.     infos[i].row = manyRows[i];
  299.     infos[i].curA = aval;
  300.     pids[i] = sproc(doit,PR_SADDR,(void *)&infos[i]);
  301.   }
  302.   real_doit(aval,manyRows[i]);
  303.  
  304.   /*
  305.    * We know we started (nprocs - 1) kids, so we wait for (nprocs - 1) kids
  306.    * to die.
  307.    */
  308.   for(i=0; i < nprocs - 1; i++) {
  309.     (void)wait((int *)0);
  310.   }
  311. }
  312. #endif
  313. #endif
  314.  
  315.  
  316. int main(int argc,char **argv) {
  317.   int c,i,j;
  318.   char *viewprogram = NULL;
  319.   char *outname = NULL;
  320.   FILE *outfile;
  321.  
  322.   char *bits;
  323.  
  324.   double curB;
  325.   int onRow,onCol;
  326. #ifdef MULTIPROC
  327.   int *manyRows[MAXPROCS];
  328. #endif
  329.   int *oneRow;
  330.  
  331.   int usepgm = 0;
  332.   int outtext = 0;
  333.  
  334.   init_colormap();
  335.   opterr = 0;
  336.  
  337.   whoami = *argv;
  338.   while((c = getopt(argc,argv,OPTIONS)) != EOF) {
  339.     switch(c) {
  340.     case 'r':
  341.       rows = atoi(optarg);
  342.       if(!rows)
  343.     usage();
  344.       break;
  345.  
  346.     case 'c':
  347.       cols = atoi(optarg);
  348.       if(!cols)
  349.     usage();
  350.       break;
  351.  
  352. #ifdef MULTIPROC
  353.     case '#':
  354.       nprocs = atoi(optarg);
  355.       if(!nprocs)
  356.     usage();
  357.       break;
  358. #endif
  359.  
  360.     case 'o':
  361.       outname = optarg;
  362.       viewprogram = NULL;
  363.       break;
  364.  
  365.     case 'v':
  366.       viewprogram = optarg;
  367.       outname = NULL;
  368.       break;
  369.  
  370.     case 'l':
  371.       read_colormap(optarg);
  372.       break;
  373.  
  374.     case 'm':
  375.       if(!sscanf(optarg,"%lf",&minLya))
  376.     usage();
  377.       break;
  378.  
  379.     case 'd':
  380.       Dwell = atoi(optarg);
  381.       if(!Dwell)
  382.     usage();
  383.       /*
  384.        * By forcing the Dwell to be a multiple of four, we can use the
  385.        * speedup noted in lya().  We *could* just do this silently,
  386.        * by if it might make a difference to someone...
  387.        */
  388.       if(Dwell % 4) {
  389.     fprintf(stderr,"Dwell value of %d being rounded off to %d.\n",
  390.         Dwell,
  391.         Dwell += 4 - (Dwell % 4));
  392.       }
  393.       break;
  394.  
  395.     case 's':
  396.       Settle = atoi(optarg);
  397.       if(!Settle)        /* for settle close to 0, use 1 */
  398.     usage();
  399.       break;
  400.  
  401.     case 'g':
  402.       if(!sscanf(optarg,"%lf",&initGuess))
  403.     usage();
  404.       break;
  405.  
  406.     case '5':
  407.       usepgm = 1;
  408.       break;
  409.  
  410.     case 'p':
  411.       showpos = 1;
  412.       break;
  413.  
  414.     case 't':
  415.       outtext = 1;
  416.       break;
  417.  
  418.     default:            /* includes '?' (and '#' when !MULTIPROC) */
  419.       usage();
  420.     }
  421.   }
  422.  
  423.   if(!argv[optind])
  424.     usage();
  425.  
  426.   /*
  427.    * Note to language nit-pickers... The code below is legal because
  428.    * the || operator is a sequence point.
  429.    */
  430.   if(!sscanf(argv[optind++],"%lf",&amin) || !argv[optind] ||
  431.      !sscanf(argv[optind++],"%lf",&amax) || !argv[optind] ||
  432.      !sscanf(argv[optind++],"%lf",&bmin) || !argv[optind] ||
  433.      !sscanf(argv[optind++],"%lf",&bmax) || !argv[optind])
  434.     usage();
  435.  
  436.   bits = argv[optind];
  437.  
  438.   if(Dwell < Settle)
  439.     Dwell = Settle;        /* a bit of sanity */
  440.  
  441.   j = strlen(bits);
  442.  
  443.   Rvec = malloc((Dwell+1) * sizeof *Rvec);
  444.   if(!Rvec) {
  445.     fputs("Cannot allocate space for Markus vector!\n",stderr);
  446.     exit(-1);
  447.   }
  448.  
  449.   for(i=0; i <= Dwell; i++) {
  450.     if(bits[i % j] == 'a')
  451.       Rvec[i] = 0;
  452.     else if(bits[i % j] == 'b')
  453.       Rvec[i] = 1;
  454.     else {
  455.       fprintf(stderr,"Invalid character \\%o ('%c') in bit string\n",
  456.           (int)(unsigned char)bits[i % j], bits[i % j]);
  457.       exit(-1);
  458.     }
  459.   }
  460.  
  461. #ifdef MULTIPROC
  462.   for(i=0; i < nprocs; i++) {
  463.     manyRows[i] = malloc(cols * sizeof *manyRows[i]);
  464.     if(!manyRows[i]) {
  465.       fputs("Cannot allocate space for row buffers!\n",stderr);
  466.       exit(-1);
  467.     }
  468.   }
  469. #else 
  470.   oneRow = malloc(cols * sizeof *oneRow);
  471.   if(!oneRow) {
  472.     fputs("Cannot allocate space for single-row buffer!\n",stderr);
  473.     exit(-1);
  474.   }
  475. #endif
  476.  
  477.   aincr = (amax - amin) / rows;
  478.   bincr = (bmax - bmin) / cols;
  479.  
  480.   if(viewprogram) {
  481.     outfile = popen(viewprogram,"w");
  482.     if(!outfile) {
  483.       fprintf(stderr,"Couldn't popen(\"%s\",\"w\")!\n",viewprogram);
  484.       exit(-1);
  485.     }
  486.   } else if(outname) {
  487.     outfile = fopen(outname,"w");
  488.     if(!outfile) {
  489.       perror(outname);
  490.       exit(-1);
  491.     }
  492.   } else {
  493.     outfile = stdout;
  494.   }
  495.  
  496.   /*
  497.    * Print P[PG]M header...
  498.    */
  499.   fprintf(outfile,"P%d %d %d\n",(usepgm?2:3) + (outtext?0:3),cols,rows);
  500.   fprintf(outfile,"%d\n",nColors-1);
  501.  
  502. #ifndef MULTIPROC
  503.   for(onRow = 0, curA = amax; onRow < rows; curA -= aincr, onRow++) {
  504.     for(onCol = cols - 1, curB = bmax; onCol; curB -= bincr, onCol--) {
  505.       oneRow[onCol] = lya(curA,curB);
  506.     }
  507.     
  508. #else
  509.   for(onRow = 0, curA = amax; onRow < rows;
  510.       curA -= aincr*nprocs, onRow += nprocs) {
  511.     DO_SPROC(manyRows);
  512.  
  513.     for(j = 0; j < nprocs; j++) {
  514.       oneRow = manyRows[j];
  515. #endif
  516.     for(i=0; i < cols; i++) {
  517.       unsigned char c;
  518.       c = (unsigned char)oneRow[i];
  519.       if(outtext) {
  520.     if(usepgm)
  521.       fprintf(outfile,"%d ",c);
  522.     else
  523.       fprintf(outfile,"%d %d %d ",colormap[c].red,colormap[c].green,
  524.           colormap[c].blue);
  525.       } else if(usepgm) {
  526.     fwrite((char *)&c,sizeof c,1,outfile);
  527.       } else {
  528.     fwrite((char *)&colormap[c],sizeof colormap[0],1,outfile);
  529.       }
  530.     }
  531. #ifdef MULTIPROC
  532.     }
  533. #endif
  534.     if(!(onRow % 16)) {
  535.       if(outtext)
  536.     fputc('\n',outfile);
  537.       fprintf(stderr,"Done row %d\n",onRow);
  538.     }
  539.   }
  540.  
  541.   if(viewprogram)
  542.     pclose(outfile);
  543.   else if(outname)
  544.     fclose(outfile);
  545.  
  546.   return 0;
  547. }
  548.