home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 275_02 / prob31.c < prev    next >
Text File  |  1979-12-31  |  20KB  |  706 lines
  1.  
  2. /* prob31.c                        */
  3. /* program for LCA31 option 't'                */
  4. /* calculate probabilities related to evolution        */
  5. /* Harold V. McIntosh, 10 August 1987            */
  6.  
  7. /* references:                            */
  8. /*                                */
  9. /* W. John Wilbur, David J. Lipman and Shihab A. Shamma        */
  10. /* On the prediction of local patterns in cellular automata    */
  11. /* Physica 19D 397-410 (1986)                    */ 
  12. /*                                */
  13. /* Howard A. Gutowitz, Jonathan D. Victor and Bruce W. Knight    */
  14. /* Local structure theory for cellular automata            */
  15. /* Physica 28D 18-48 (1987)                    */
  16.  
  17. /*    Copyright (C) 1987    */
  18. /*    Copyright (C) 1988    */
  19. /*    Harold V. McIntosh    */
  20. /*    Gerardo Cisneros S.    */
  21.  
  22. # define BROW       13   /* row for bar charts        */
  23. # define EROW       22   /* row for evolution synopsis    */
  24. # define TROW      99   /* column for t-triangle    */
  25. # define TCOL     205   /* column for t-triangle    */
  26.  
  27. /* edit the probability screen    */
  28. edtri() {char c;
  29.   videomode(COLGRAF);
  30.   videopalette(YELREGR);
  31.  
  32.   while (0<1) {
  33.   videocursor(0,0,0);
  34.   scrrul();
  35.   videocursor(0,0,36);
  36.   videoputc('?',2);
  37.   c=kbdin();
  38.   if (c == '\015') break;
  39.   videocursor(0,0,38);
  40.   videoputc(c,2);
  41.   videocursor(0,0,36);
  42.   videoputc(' ',0);
  43.   switch (c) {
  44.     case '+': videopalette(WHCYMAG); break;
  45.     case '-': videopalette(YELREGR); break;
  46.     case 'a': nblclr(); asfreq(TROW,TCOL,3); break;
  47.     case 'e': pevolve(); break;
  48.     case 'g': lifreq(50,TROW,TCOL,2); break;
  49.     case 'G': lifreq(200,TROW,TCOL,1); break;
  50.     case 'm': nblclr(); moncar(); break;
  51.     case 's': spdiff(50,6); break;
  52.     case 't': pdiff(100,6); break;
  53.     case '1': nblclr(); oneblfreq(8*BROW,300,48); break;
  54.     case '2': nblclr(); twoblfreq(8*BROW,300,48); break;
  55.     case '3': nblclr(); thrblfreq(8*BROW,300,48); break;
  56.     case '4': nblclr(); foublfreq(8*BROW,300,48); break;
  57.     case '5': nblclr(); fivblfreq(8*BROW,300,48); break;
  58.     case 'z': nblclr(); break;
  59.     case '/': videomode(COLGRAF); videopalette(YELREGR); break;
  60.     case '?': nblclr(); trmenu(); break;
  61.     default: break;
  62.     }; /* end switch */
  63.   };   /* end while  */
  64.   videopalette(WHCYMAG);
  65.   videomode(T80X25);
  66. }      /* end edtri  */
  67.  
  68. /* show menu */
  69. trmenu() {
  70.   videocursor(0,BROW,0);
  71.   printf("a     - a priori estimates\n");
  72.   printf("m,g,G - sample evolution\n");
  73.   printf("s     - graph probabilities in stereo\n");
  74.   printf("t     - graph 2 block probabilities\n");
  75.   printf("12345 - n-block bar charts\n");
  76.   printf("+-    - change color pallette\n");
  77.   printf("e     - 16 lines evolution\n");
  78.   printf("z/?   - clear panel, screen; show menu\n");
  79.   printf("<cr>  - exit\n");
  80. }
  81.  
  82. /* show sixteen lines of evolution on screen */
  83. pevolve() {int i, j;
  84.   videoscroll(EROW,0,EROW+1,40,0,0);
  85.   asctobin();
  86.   for (j=8*EROW; j<8*EROW+16; j++) {
  87.     for (i=0; i<AL; i++) videodot(j,i,arr1[i]);
  88.     onegen(AL);
  89.     };
  90. }
  91.  
  92. /* Clear a space for the n-block frequencies */
  93. nblclr() {videoscroll(BROW,0,BROW+8,40,0,0);}
  94.  
  95. /* make a frame for a graph          */
  96. /* (x,y) = lower left corner; e.g. (0,0) */
  97. /* n     = dimension of frame            */
  98. gfram(x,y,n) int x, y, n; {int i;
  99.  
  100. for (i=0; i<=n; i++) videodot(199-y-i,x,0);
  101. for (i=0; i<=n; i++) videodot(199-y-i,x+n,0);
  102. for (i=0; i<=n; i++) videodot(199-n-y,x+i,0);
  103. for (i=0; i<=n; i++) videodot(199-y,x+i,0);
  104.  
  105. for (i=0; i<=n; i+=10) videodot(199-y-i,x,3);
  106. for (i=0; i<=n; i+=10) videodot(199-y-i,x+n,3);
  107. for (i=0; i<=n; i+=10) videodot(199-n-y,x+i,3);
  108. for (i=0; i<=n; i+=10) videodot(199-y,x+i,3);
  109. }
  110.  
  111. /* display state frequencies in arr1 as a dot in triangle at (u,v) */
  112. /* n = generations of evolution                       */
  113. /* (u,v) = corner of triangle                       */
  114. /* l = color of dot                           */
  115. lifreq(n,u,v,l) int n, u, v, l; {
  116. int i, ii;
  117. int stat[KK];
  118. double staf[KK], s;
  119.  
  120. s=1.0/((double)(AL));
  121. asctobin();
  122. for (ii=0; ii<n; ii++) {
  123.   onegen(AL);
  124.   for (i=0; i<KK; i++) stat[i]=0;
  125.   for (i=0; i<AL; i++) (stat[arr1[i]])++;
  126.   for (i=0; i<KK; i++) staf[i]=s*((double)stat[i]);
  127.   videotrdot(u,v,staf[0],staf[1],staf[2],l);
  128.   };
  129. }
  130.  
  131. /* plot a single point on a triangular grid */
  132. videotrdot(u,v,x,y,z,l) double x, y, z; int u, v, l; {double s;
  133. s=199.0/(x+y+z); x=x*s; y=y*s; z=z*s;
  134. videodot(u-(int)(0.433*z),v+(int)(0.250*(y+y+z)),l);
  135. }
  136.  
  137. /* plot a stereo point on a triangular grid */
  138. /* (u,v) screen location lower left corner  */
  139. /* (x,y,z) triangular coordinates        */
  140. /* t       height                */
  141. /* l       color                */
  142. videostdot(u,v,x,y,z,t,l) double x, y, z, t; int u, v, l; {double s;
  143. s=199.0/(x+y+z); x=x*s; y=y*s; z=z*s;
  144. videodot(u-(int)(0.433*z),v+(int)(0.250*(y+y+z+t)),l);
  145. }
  146.  
  147. /* set up a triangular gridwork */
  148. /* (u,v) = corner of triangle    */
  149. /* n = number of rulings    */
  150. /* l = color of lines        */
  151. videotrgrid(u,v,n,l) int u, v, n, l; {
  152. int i, j;
  153. double a, b, c;
  154. if(n==0) return;
  155. for (i=0; i<=n;  i++) {
  156. for (j=0; i+j<n; j++) {
  157.   a=((double)(i)/(double)(n));
  158.   b=((double)(j)/(double)(n));
  159.   c=1.0-a-b;
  160.   videotrdot(u,v,a,b,c,l);
  161.   };};
  162. }
  163.  
  164. /* "Monte Carlo" determination of probabilities */
  165. moncar() {
  166. int i, b[KK];
  167. double bf[KK];
  168.   asctobin();
  169.   onegen(AL);
  170.   for (i=0; i<KK; i++) b[i]=0;
  171.   for (i=0; i<AL; i++) b[arr1[i]]++;
  172.   for (i=0; i<KK; i++) bf[i]=((double)(b[i]))/((double)(AL));
  173.   videocursor(0,BROW+8,0);
  174.   printf("(montecarlo) "); 
  175.   for (i=0; i<KK; i++) printf("%2d:%5.3f",i,bf[i]);
  176. }
  177.  
  178. /* Generate coefficients of the Bernstein Polynomial */
  179. berncoef() {int i, i0, i1, i2;
  180. for ( i=0;  i<KK;  i++) {
  181. for (i0=0; i0<KK; i0++) {
  182. for (i1=0; i1<KK; i1++) {
  183. for (i2=0; i2<KK; i2++) {
  184.   bp[i][i0][i1][i2]=0.0;
  185.   };};};};
  186. for (i0=0; i0<KK; i0++) {
  187. for (i1=0; i1<KK; i1++) {
  188. for (i2=0; i2<KK; i2++) {
  189.   bp[ascrule[i0][i1][i2]-'0'][i0][i1][i2]+=1.0;
  190.   };};};
  191. }
  192.  
  193. /* evaluate the nth generation Bernstein polynomial at point p */
  194. double bern(i,x,y,z) int i; double x, y, z; {double s;
  195. s=x*x*x*bp[i][0][0][0]+y*y*y*bp[i][1][1][1]+z*z*z*bp[i][2][2][2];
  196. s+=x*x*y*(bp[i][0][0][1]+bp[i][0][1][0]+bp[i][1][0][0]);
  197. s+=x*x*z*(bp[i][0][0][2]+bp[i][0][2][0]+bp[i][2][0][0]);
  198. s+=x*y*y*(bp[i][0][1][1]+bp[i][1][0][1]+bp[i][1][1][0]);
  199. s+=x*z*z*(bp[i][0][2][2]+bp[i][2][0][2]+bp[i][2][2][0]);
  200. s+=y*y*z*(bp[i][1][1][2]+bp[i][1][2][1]+bp[i][2][1][1]);
  201. s+=y*z*z*(bp[i][1][2][2]+bp[i][2][1][2]+bp[i][2][2][1]);
  202. s+=x*y*z*(bp[i][0][1][2]+bp[i][0][2][1]+bp[i][1][0][2]+
  203.       bp[i][1][2][0]+bp[i][2][0][1]+bp[i][2][1][0]);
  204. return s;
  205. }
  206.  
  207. /* graph the probability differential for state l over a triangle */
  208. pdiff(n,l) int n, l; {int i, j; double a, b, c, en, t, sqsu();
  209. if (n==0) return;
  210. en=1.0/((double)(n));
  211. berncoef();
  212. for (i=0; i<=n; i++) {
  213. for (j=0; i+j<=n; j++) {
  214.   a=((double)(i))*en;
  215.   b=((double)(j))*en;
  216.   c=1.0-a-b;
  217.   t=sqsu(a,b,c);
  218.   videotrdot(99,205,a,b,c,sqco(t,i));
  219.   if (kbdst()) {kbdin(); return;};
  220.   };};
  221. }
  222.  
  223. /* stereo graph of probability differential for state l over a triangle */
  224. spdiff(n,l) int n, l; {int i, j; double a, b, c, en, t, sqsu();
  225. if (n==0) return;
  226. en=1.0/((double)(n));
  227. berncoef();
  228. for (i=0; i<=n; i++) {
  229. for (j=0; i+j<=n; j++) {
  230.   a=((double)(i))*en;
  231.   b=((double)(j))*en;
  232.   c=1.0-a-b;
  233.   t=sqsu(a,b,c);
  234.   videostdot(99,5,a,b,c,-12.0*t,sqco(t,i));
  235.   videostdot(99,105,a,b,c,12.0*t,sqco(t,i));
  236.   if (kbdst()) {kbdin(); return;};
  237.   };};
  238. }
  239.  
  240. /* calculate the squaresum of the three polynomials at a point */
  241. double sqsu(x,y,z) double x, y, z; {double u, v, w, bern();
  242. u=bern(0,x,y,z)-x;
  243. v=bern(1,x,y,z)-y;
  244. w=bern(2,x,y,z)-z;
  245. return (u*u+v*v+w*w);
  246. }
  247.  
  248. /* generate squaresum compatible contour values */
  249. int sqco(t,i) int i; double t; {int l;
  250. if (t<0.0075) l=0; else {
  251. if (t<0.0750) {if (i%2) l=0; else l=3;} else {
  252. if (t<0.1250) l=3; else {
  253. if (t<0.2500) {if (i%2) l=3; else l=2;} else {
  254. if (t<0.5000) l=2; else {
  255. if (t<0.7500) {if (i%2) l=2; else l=1;} else {
  256. if (t<1.0000) l=1; else {
  257.           if (i%2) l=1; else l=0;
  258.    }}}}}}}
  259. return l;
  260. }
  261.  
  262. /* graph the frequencies of ascrule in color l */
  263. /* put dot at coordinates (u,v) */
  264. asfreq(u,v,l) int u, v, l; {
  265. int    i, j, i0, i1, i2;
  266. int    stat[KK], stal, stac, star;        /* statistic counts */
  267. double staf[KK], pp;
  268.  
  269. pp=1.0/((double)(KK*KK*KK));
  270. stal=0; stac=0; star=0;
  271. for (i=0; i<KK; i++) stat[i]=0;
  272. for (i0=0; i0<KK; i0++) {
  273. for (i1=0; i1<KK; i1++) {
  274. for (i2=0; i2<KK; i2++) {
  275.   j=ascrule[i0][i1][i2]-'0';
  276.   stat[j]++;
  277.   if(j==i0) star++;
  278.   if(j==i1) stac++;
  279.   if(j==i2) stal++;
  280.   };};};
  281. videocursor(0,BROW,0);
  282. for (i=0; i<KK; i++) printf("%1d    - %4.2f\n",i,((double)stat[i])*pp);
  283. printf("\n");
  284. printf("left  - %4.2f\n",((double)stal)*pp);
  285. printf("still - %4.2f\n",((double)stac)*pp);
  286. printf("right - %4.2f\n",((double)star)*pp);
  287. for (i=0; i<KK; i++) staf[i]=((double)stat[i])*pp;
  288. videotrdot(u,v,staf[0],staf[1],staf[2],l);
  289. }
  290.  
  291. /* evaluate the one-block probabilities after one generation */
  292. onebl(x,a) double x[KK], a[KK]; {int i0, i1, i2, j0;
  293.  
  294. for (j0=0; j0<KK; j0++) x[j0]=0.0;
  295.  
  296. for (i0=0; i0<KK; i0++) {
  297. for (i1=0; i1<KK; i1++) {
  298. for (i2=0; i2<KK; i2++) {
  299.   j0=ascrule[i0][i1][i2]-'0';
  300.   x[j0]+=a[i0]*a[i1]*a[i2];
  301.   };};};
  302. }
  303.  
  304. /* iterate the 1-block parameters to find self-consistent values */
  305. /* graph the iterative steps in bar-chart form             */
  306. /* ll - initial line    */
  307. /* mm - length of line    */
  308. /* nn - number of lines    */
  309. oneblfreq(ll,mm,nn) int ll, mm, nn; {
  310. int    ii, i, l, m, n;
  311. double op[KK], np[KK];
  312. double d, e, f, s;
  313.  
  314. m=0;
  315. f=(double)mm;
  316. s=1.0/((double)(KK));
  317. n=(int)(f*s);
  318. videodot(ll,m++,3);
  319. for (i=0; i<KK; i++) {op[i]=s; for (l=0; l<n; l++) videodot(ll,m++,i);};
  320.  
  321. for (ii=1; ii<=nn; ii++) {
  322.   e=0.0; m=0;
  323.   onebl(np,op);
  324.   videodot(ll+ii,m++,3);
  325.   for (i=0; i<KK; i++) {
  326.     n=(int)(f*np[i]); if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,i);
  327.     d=op[i]-np[i];
  328.     e+=d*d;
  329.     op[i]=np[i];
  330.     };
  331.   videodot(ll+ii,m++,3);
  332.   if (op[0]<=0.001) break;
  333.   if (op[1]<=0.001) break;
  334.   if (e<=0.0000001) break;
  335.   if (kbdst()) {kbdin(); break;};
  336.   };
  337. videocursor(0,BROW+8,0);
  338. printf("(1-block) "); 
  339. for (i=0; i<KK; i++) printf("%2d:%5.3f ",i,op[i]);
  340. videotrdot(TROW,TCOL,op[0],op[1],op[2],l);
  341. }
  342.  
  343. /* evaluate the two-block probabilities after one generation */
  344. twobl(x,a) double x[KK][KK], a[KK][KK]; {
  345. int    i0, i1, i2, i3;
  346. int    j0, j1;
  347. double w, b[KK];
  348.  
  349. for (j0=0; j0<KK; j0++) {for (j1=0; j1<KK; j1++) x[j0][j1]=0.0;};
  350.  
  351. for (j0=0; j0<KK; j0++) {b[j0]=0.0; for (j1=0; j1<KK; j1++) b[j0]+=a[j0][j1];};
  352.  
  353. for (i0=0; i0<KK; i0++) {
  354. for (i1=0; i1<KK; i1++) {
  355. for (i2=0; i2<KK; i2++) {
  356. for (i3=0; i3<KK; i3++) {
  357.   j0=ascrule[i0][i1][i2]-'0';
  358.   j1=ascrule[i1][i2][i3]-'0';
  359.   w=a[i0][i1]*a[i1][i2]*a[i2][i3];
  360.   if (w!=0.0) w/=b[i1]*b[i2];
  361.   x[j0][j1]+=w;
  362.   };};};};
  363. }
  364.  
  365. /* iterate the 2-block parameters to find self-consistent values */
  366. /* graph the iterative steps in bar-chart form             */
  367. /* ll - initial line    */
  368. /* mm - length of line    */
  369. /* nn - number of lines    */
  370. twoblfreq(ll,mm,nn) int ll, mm, nn; {
  371. int    ii, i, j, l, m, n;
  372. double op[KK][KK], np[KK][KK];
  373. double b[KK], d, e, f, s;
  374.  
  375. m=0;
  376. f=(double)mm;
  377. s=1.0/((double)(KK*KK));
  378. n=(int)(f*s);
  379. videodot(ll,m++,3);
  380. for (i=0; i<KK; i++) {
  381. for (j=0; j<KK; j++) {
  382.   op[i][j]=s;
  383.   for (l=0; l<n; l++) videodot(ll,m++,j);
  384.   };};
  385. videodot(ll,m++,3);
  386.  
  387. for (ii=1; ii<=nn; ii++) {
  388.   e=0.0; m=0;
  389.   twobl(np,op);
  390.   videodot(ll+ii,m++,3);
  391.   for (i=0; i<KK; i++) {
  392.   for (j=0; j<KK; j++) {
  393.     n=(int)(f*np[i][j]);
  394.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,j);
  395.     d=op[i][j]-np[i][j];
  396.     e+=d<0.0?-d:d;
  397.     op[i][j]=np[i][j];
  398.     };};
  399.   videodot(ll+ii,m++,3);
  400. for (i=0; i<KK; i++) {b[i]=0.0; for (j=0; j<KK; j++) b[i]+=op[i][j];};
  401. videotrdot(TROW,TCOL,b[0],b[1],b[2],1);
  402.   if (e<=0.0001) break;
  403.   if (kbdst()) {kbdin(); break;};
  404.   };
  405. videocursor(0,BROW+8,0);
  406. printf("(2-block) "); 
  407. for (i=0; i<KK; i++) printf("%2d:%5.3f ",i,b[i]);
  408. videotrdot(TROW,TCOL,b[0],b[1],b[2],3);
  409. }
  410.  
  411. /* evaluate the three-block probabilities after one generation */
  412. thrbl(x,a) double x[KK][KK][KK], a[KK][KK][KK]; {
  413. int    i0, i1, i2, i3, i4;
  414. int    j0, j1, j2;
  415. double w, b[KK][KK];
  416.  
  417. for (j0=0; j0<KK; j0++) {
  418. for (j1=0; j1<KK; j1++) {
  419.   for (j2=0; j2<KK; j2++) x[j0][j1][j2]=0.0;
  420.   b[j0][j1]=0.0;
  421.   for (j2=0; j2<KK; j2++) b[j0][j1]+=a[j0][j1][j2];  
  422.   };};
  423.  
  424. for (i0=0; i0<KK; i0++) {
  425. for (i1=0; i1<KK; i1++) {
  426. for (i2=0; i2<KK; i2++) {
  427. for (i3=0; i3<KK; i3++) {
  428. for (i4=0; i4<KK; i4++) {
  429.   j0=ascrule[i0][i1][i2]-'0';
  430.   j1=ascrule[i1][i2][i3]-'0';
  431.   j2=ascrule[i2][i3][i4]-'0';
  432.   w=a[i0][i1][i2]*a[i1][i2][i3]*a[i2][i3][i4];
  433.   if (w!=0.0) w/=b[i1][i2]*b[i2][i3];
  434.   x[j0][j1][j2]+=w;
  435.   };};};};};
  436. }
  437.  
  438. /* iterate the 3-block parameters to find self-consistent values */
  439. /* ll - initial line    */
  440. /* mm - length of line    */
  441. /* nn - number of lines    */
  442. thrblfreq(ll,mm,nn) int ll, mm, nn; {
  443. int    ii, i, j, k, l, m, n;
  444. double op[KK][KK][KK], np[KK][KK][KK];
  445. double b[KK], bb[KK][KK], d, e, f, s;
  446.  
  447. m=0;
  448. f=(double)mm;
  449. s=1.0/((double)(KK*KK*KK));
  450. n=(int)(f*s);
  451. videodot(ll,m++,3);
  452. for (i=0; i<KK; i++) {
  453. for (j=0; j<KK; j++) {
  454. for (k=0; k<KK; k++) {
  455.   op[i][j][k]=s;
  456.   for (l=0; l<n; l++) videodot(ll,m++,k);
  457.   };};};
  458. videodot(ll,m++,3);
  459.  
  460. for (ii=1; ii<=nn; ii++) {
  461.   e=0.0; m=0;
  462.   thrbl(np,op);
  463.   videodot(ll+ii,m++,3);
  464.   for (i=0; i<KK; i++) {
  465.   for (j=0; j<KK; j++) {
  466.   for (k=0; k<KK; k++) {
  467.     n=(int)(f*np[i][j][k]);
  468.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,k);
  469.     d=op[i][j][k]-np[i][j][k];
  470.     e+=d<0.0?-d:d;
  471.     op[i][j][k]=np[i][j][k];
  472.     };};};
  473.   videodot(ll+ii,m++,3);
  474.   if (e<=0.0001) break;
  475.   if (kbdst()) {kbdin(); break;};
  476.   };
  477.  
  478. for (i=0; i<KK; i++) {
  479.   b[i]=0.0;
  480.   for (j=0; j<KK; j++) {
  481.   for (k=0; k<KK; k++) {
  482.     b[i]+=op[i][j][k];
  483.     };}; };
  484. videocursor(0,BROW+8,0);
  485. printf("(3-block) "); 
  486. for (i=0; i<KK; i++) printf("%2d:%5.3f ",i,b[i]);
  487. videotrdot(TROW,TCOL,op[0],op[1],op[2],l);
  488.  
  489. for (i=0; i<KK; i++) {
  490. for (j=0; j<KK; j++) {
  491.   bb[i][j]=0.0;
  492.   for (k=0; k<KK; k++) {
  493.     bb[i][j]+=op[i][j][k];
  494.     };}; };
  495.  
  496. /*videocursor(0,BROW+9,0);*/
  497. /*for (i=0; i<KK; i++) for (j=0; j<KK; j++)*/ 
  498. /*  printf("%1d%1d:%5.3f ",i,j,bb[i][j]);*/
  499. /*videodot(199-(int)(100.0*bb[1][1]),BORG+(int)(100.0*b[1]),2);*/
  500. }
  501.  
  502. /* evaluate the four-block probabilities after one generation */
  503. foubl(x,a) double x[KK][KK][KK][KK], a[KK][KK][KK][KK]; {
  504. int    i0, i1, i2, i3, i4, i5;
  505. int    j0, j1, j2, j3;
  506. double w, b[KK][KK][KK];
  507.  
  508. for (j0=0; j0<KK; j0++) {
  509. for (j1=0; j1<KK; j1++) {
  510. for (j2=0; j2<KK; j2++) {
  511.   for (j3=0; j3<KK; j3++) x[j0][j1][j2][j3]=0.0;
  512.   b[j0][j1][j2]=0.0;
  513.   for (j3=0; j3<KK; j3++) b[j0][j1][j2]+=a[j0][j1][j2][j3];  
  514.   };};};
  515.  
  516. for (i0=0; i0<KK; i0++) {
  517. for (i1=0; i1<KK; i1++) {
  518. for (i2=0; i2<KK; i2++) {
  519. for (i3=0; i3<KK; i3++) {
  520. for (i4=0; i4<KK; i4++) {
  521. for (i5=0; i5<KK; i5++) {
  522.   j0=ascrule[i0][i1][i2]-'0';
  523.   j1=ascrule[i1][i2][i3]-'0';
  524.   j2=ascrule[i2][i3][i4]-'0';
  525.   j3=ascrule[i3][i4][i5]-'0';
  526.   w=a[i0][i1][i2][i3]*a[i1][i2][i3][i4]*a[i2][i3][i4][i5];
  527.   if (w!=0.0) w/=b[i1][i2][i3]*b[i2][i3][i4];
  528.   x[j0][j1][j2][j3]+=w;
  529.   };};};};};};
  530. }
  531.  
  532. /* iterate the 4-block parameters to find self-consistent values */
  533. /* ll - initial line    */
  534. /* mm - length of line    */
  535. /* nn - number of lines    */
  536. foublfreq(ll,mm,nn) int ll, mm, nn; {
  537. int    ii, i0, i1, i2, i3, l, m, n;
  538. double op[KK][KK][KK][KK], np[KK][KK][KK][KK];
  539. double b[KK], bb[KK][KK], d, e, f, s;
  540.  
  541. m=0;
  542. f=(double)mm;
  543. s=1.0/((double)(KK*KK*KK*KK));
  544. n=(int)(f*s);
  545. videodot(ll,m++,3);
  546. for (i0=0; i0<KK; i0++) {
  547. for (i1=0; i1<KK; i1++) {
  548. for (i2=0; i2<KK; i2++) {
  549. for (i3=0; i3<KK; i3++) {
  550.   op[i0][i1][i2][i3]=s;
  551.   for (l=0; l<n; l++) videodot(ll,m++,i3);
  552.   };};};};
  553. videodot(ll,m++,3);
  554.  
  555. for (ii=1; ii<=nn; ii++) {
  556.   e=0.0; m=0;
  557.   foubl(np,op);
  558.   videodot(ll+ii,m++,3);
  559.   for (i0=0; i0<KK; i0++) {
  560.   for (i1=0; i1<KK; i1++) {
  561.   for (i2=0; i2<KK; i2++) {
  562.   for (i3=0; i3<KK; i3++) {
  563.     n=(int)(f*np[i0][i1][i2][i3]);
  564.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,i3);
  565.     d=op[i0][i1][i2][i3]-np[i0][i1][i2][i3];
  566.     e+=d<0.0?-d:d;
  567.     op[i0][i1][i2][i3]=np[i0][i1][i2][i3];
  568.     };};};};
  569.   videodot(ll+ii,m++,3);
  570.   if (e<=0.0001) break;
  571.   if (kbdst()) {kbdin(); break;};
  572.   };
  573.  
  574. for (i0=0; i0<KK; i0++) {
  575.   b[i0]=0.0;
  576.   for (i1=0; i1<KK; i1++) {
  577.   for (i2=0; i2<KK; i2++) {
  578.   for (i3=0; i3<KK; i3++) {
  579.     b[i0]+=op[i0][i1][i2][i3];
  580.     };};}; };
  581.  
  582. videocursor(0,BROW+8,0);
  583. printf("(4-block) "); 
  584. for (i0=0; i0<KK; i0++) printf("%2d:%5.3f ",i0,b[i0]);
  585. videotrdot(TROW,TCOL,op[0],op[1],op[2],l);
  586.  
  587. for (i0=0; i0<KK; i0++) {
  588. for (i1=0; i1<KK; i1++) {
  589.   bb[i0][i1]=0.0;
  590.   for (i2=0; i2<KK; i2++) {
  591.   for (i3=0; i3<KK; i3++) {
  592.     bb[i0][i1]+=op[i0][i1][i2][i3];
  593.     };}; };};
  594.  
  595. /*videocursor(0,BROW+9,0);*/
  596. /*for (i0=0; i0<KK; i0++) for (i1=0; i1<KK; i1++)*/ 
  597. /*  printf("%1d%1d:%5.3f ",i0,i1,bb[i0][i1]);*/
  598. /*videodot(199-(int)(100.0*bb[1][1]),BORG+(int)(100.0*b[1]),2);*/
  599. }
  600.  
  601. /* evaluate the five-block probabilities after one generation */
  602. fivbl(x,a) double x[KK][KK][KK][KK][KK], a[KK][KK][KK][KK][KK]; {
  603. int    i0, i1, i2, i3, i4, i5, i6;
  604. int    j0, j1, j2, j3, j4;
  605. double w, b[KK][KK][KK][KK];
  606.  
  607. for (j0=0; j0<KK; j0++) {
  608. for (j1=0; j1<KK; j1++) {
  609. for (j2=0; j2<KK; j2++) {
  610. for (j3=0; j3<KK; j3++) {
  611.   for (j4=0; j4<KK; j4++) x[j0][j1][j2][j3][j4]=0.0;
  612.   b[j0][j1][j2][j3]=0.0;
  613.   for (j4=0; j4<KK; j4++) b[j0][j1][j2][j3]+=a[j0][j1][j2][j3][j4];  
  614.   };};};};
  615.  
  616. for (i0=0; i0<KK; i0++) {
  617. for (i1=0; i1<KK; i1++) {
  618. for (i2=0; i2<KK; i2++) {
  619. for (i3=0; i3<KK; i3++) {
  620. for (i4=0; i4<KK; i4++) {
  621. for (i5=0; i5<KK; i5++) {
  622. for (i6=0; i6<KK; i6++) {
  623.   j0=ascrule[i0][i1][i2]-'0';
  624.   j1=ascrule[i1][i2][i3]-'0';
  625.   j2=ascrule[i2][i3][i4]-'0';
  626.   j3=ascrule[i3][i4][i5]-'0';
  627.   j4=ascrule[i4][i5][i6]-'0';
  628.   w=a[i0][i1][i2][i3][i4]*a[i1][i2][i3][i4][i5]*a[i2][i3][i4][i5][i6];
  629.   if (w!=0.0) w/=b[i1][i2][i3][i4]*b[i2][i3][i4][i5];
  630.   x[j0][j1][j2][j3][j4]+=w;
  631.   };};};};};};};
  632. }
  633.  
  634. /* iterate the 5-block parameters to find self-consistent values */
  635. /* ll - initial line    */
  636. /* mm - length of line    */
  637. /* nn - number of lines    */
  638. fivblfreq(ll,mm,nn) int ll, mm, nn; {
  639. int    ii, i0, i1, i2, i3, i4, l, m, n;
  640. double op[KK][KK][KK][KK][KK], np[KK][KK][KK][KK][KK];
  641. double b[KK], bb[KK][KK], d, e, f, s;
  642.  
  643. m=0;
  644. f=(double)mm;
  645. s=1.0/((double)(KK*KK*KK*KK*KK));
  646. n=(int)(f*s);
  647. videodot(ll,m++,3);
  648. for (i0=0; i0<KK; i0++) {
  649. for (i1=0; i1<KK; i1++) {
  650. for (i2=0; i2<KK; i2++) {
  651. for (i3=0; i3<KK; i3++) {
  652. for (i4=0; i4<KK; i4++) {
  653.   op[i0][i1][i2][i3][i4]=s;
  654.   for (l=0; l<n; l++) videodot(ll,m++,i4);
  655.   };};};};};
  656. videodot(ll,m++,3);
  657.  
  658. for (ii=1; ii<=nn; ii++) {
  659.   e=0.0; m=0;
  660.   fivbl(np,op);
  661.   videodot(ll+ii,m++,3);
  662.   for (i0=0; i0<KK; i0++) {
  663.   for (i1=0; i1<KK; i1++) {
  664.   for (i2=0; i2<KK; i2++) {
  665.   for (i3=0; i3<KK; i3++) {
  666.   for (i4=0; i4<KK; i4++) {
  667.     n=(int)(f*np[i0][i1][i2][i3][i4]);
  668.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,i4);
  669.     d=op[i0][i1][i2][i3][i4]-np[i0][i1][i2][i3][i4];
  670.     e+=d<0.0?-d:d;
  671.     op[i0][i1][i2][i3][i4]=np[i0][i1][i2][i3][i4];
  672.     };};};};};
  673.   videodot(ll+ii,m++,3);
  674.   if (e<=0.0001) break;
  675.   if (kbdst()) {kbdin(); break;};
  676.   };
  677. for (i0=0; i0<KK; i0++) {
  678.   b[i0]=0.0;
  679.   for (i1=0; i1<KK; i1++) {
  680.   for (i2=0; i2<KK; i2++) {
  681.   for (i3=0; i3<KK; i3++) {
  682.   for (i4=0; i4<KK; i4++) {
  683.     b[i0]+=op[i0][i1][i2][i3][i4];
  684.     };};};}; };
  685. videocursor(0,BROW+8,0);
  686. printf("(5-block) "); 
  687. for (i0=0; i0<KK; i0++) printf("%2d:%5.3f ",i0,b[i0]);
  688. videotrdot(TROW,TCOL,op[0],op[1],op[2],l);
  689.  
  690. for (i0=0; i0<KK; i0++) {
  691. for (i1=0; i1<KK; i1++) {
  692.   bb[i0][i1]=0.0;
  693.   for (i2=0; i2<KK; i2++) {
  694.   for (i3=0; i3<KK; i3++) {
  695.   for (i4=0; i4<KK; i4++) {
  696.     bb[i0][i1]+=op[i0][i1][i2][i3][i4];
  697.     };}; };};};
  698.  
  699. /*videocursor(0,BROW+9,0);*/
  700. /*for (i0=0; i0<KK; i0++) for (i1=0; i1<KK; i1++) */
  701. /*  printf("%1d%1d:%5.3f ",i0,i1,bb[i0][i1]);*/
  702. /*videodot(199-(int)(100.0*bb[1][1]),BORG+(int)(100.0*b[1]),2);*/
  703. }
  704.  
  705. /* end */
  706.