home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 275_02 / prob21.c < prev    next >
Text File  |  1979-12-31  |  29KB  |  1,018 lines
  1.  
  2. /* prob21.c                            */
  3. /* program for lcau21 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        3    /* row for bar charts        */
  23. # define EROW        1    /* row for evolution synopsis    */
  24. # define AORG         0    /* x-origin  of contour plot    */
  25. # define BORG       109    /* x-origin  of 2-block param   */
  26. # define CORG       219    /* x-origin  of Bernstein plot  */
  27.  
  28. /* edit the probability screen    */
  29. edtri() {char c;
  30.   videomode(COLGRAF);
  31.   videopalette(YELREGR);
  32.  
  33.   while (0<1) {
  34.   woruno(0,28);
  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': asfreq(3); break;
  47.     case 'e': pevolve(); break;
  48.     case 'g': lifreq(50,2); break;
  49.     case 'G': lifreq(200,1); break;
  50.     case 'm': moncar(1,2); break;
  51.     case 'M': moncar(50,1); break;
  52.     case 'x': pdiff(100); break;
  53.     case 'i': pidiff(100); break;
  54.     case 'j': pjdiff(100); break;
  55.     case 'y': pdiff3(100); break;
  56.     case 'z': pdiff4(100); break;
  57.     case 'w': pdiff5(100); break;
  58.     case 'v': pdiff6(100); break;
  59.     case 't': twoblockfreq(100); break;
  60.     case '1': nblclr(); oneblfreq(8*BROW,300,48); break;
  61.     case '2': nblclr(); twoblfreq(8*BROW,300,48); break;
  62.     case '3': nblclr(); thrblfreq(8*BROW,300,48); break;
  63.     case '4': nblclr(); foublfreq(8*BROW,300,48); break;
  64.     case '5': nblclr(); fivblfreq(8*BROW,300,48); break;
  65.     case '6': nblclr(); sixblfreq(8*BROW,300,48); break;
  66.     case '/': videomode(COLGRAF); videopalette(YELREGR); break;
  67.     case '?': trmenu(); break;
  68.     default: break;
  69.     }; /* end switch */
  70.   };   /* end while  */
  71.   videopalette(WHCYMAG);
  72.   videomode(T80X25);
  73. }      /* end edtri  */
  74.  
  75. /* show menu */
  76. trmenu() {
  77.   videoscroll(BROW,0,BROW+8,40,0,0);
  78.   videocursor(0,BROW,0);
  79.   printf("a       - a priori estimates\n");
  80.   printf("m,M,g,G - sample evolution\n");
  81.   printf("xyzwv   - selfconsistent probabilities\n");
  82.   printf("xij     - iterated s-c probabilities\n");
  83.   printf("t       - graph 2 block probabilities\n");
  84.   printf("123456  - n-block bar charts\n");
  85.   printf("+-      - change color pallette\n");
  86.   printf("e       - 12 lines evolution\n");
  87.   printf("/?(clear screen, show menu), <cr>(exit)\n");
  88. }
  89.  
  90. /* show fourteen lines of evolution at top of screen */
  91. pevolve() {int i, j;
  92.   videoscroll(EROW,0,EROW+1,40,0,0);
  93.   asctobin();
  94.   for (j=8*EROW; j<8*(EROW+2)-2; j++) {
  95.     for (i=0; i<AL; i++) videodot(j,i,arr1[i]);
  96.     onegen(AL);
  97.     };
  98. }
  99.  
  100. /* Clear a space for the n-block frequencies */
  101. nblclr() {videoscroll(BROW,0,BROW+8,40,0,0);}
  102.  
  103. /* make a frame for a graph          */
  104. /* (x,y) = lower left corner; e.g. (0,0) */
  105. /* n     = dimension of frame            */
  106. gfram(x,y,n) int x, y, n; {int i;
  107.  
  108. for (i=0; i<=n; i++) videodot(199-y-i,x,0);
  109. for (i=0; i<=n; i++) videodot(199-y-i,x+n,0);
  110. for (i=0; i<=n; i++) videodot(199-n-y,x+i,0);
  111. for (i=0; i<=n; i++) videodot(199-y,x+i,0);
  112.  
  113. for (i=0; i<=n; i+=10) videodot(199-y-i,x,3);
  114. for (i=0; i<=n; i+=10) videodot(199-y-i,x+n,3);
  115. for (i=0; i<=n; i+=10) videodot(199-n-y,x+i,3);
  116. for (i=0; i<=n; i+=10) videodot(199-y,x+i,3);
  117. }
  118.  
  119. /* put a diagonal in a graph */
  120. gdiag(x,y,n) int x, y, n; {int i;
  121. for (i=0; 2*i<n; i++) videodot(199-y-2*i,x+2*i,3);
  122. }
  123.  
  124. /* graph Bernstein polynomial */
  125. bgraf(x,y,k,n) int x, y, k, n; {int i; double bern(), en, p;
  126. en=(double)(n);
  127. for (i=0; i<n; i++) {
  128.   p=((double)(i))/en;
  129.   videodot(199-y-(int)(en*bern(p,k)),x+i,1);
  130.   };
  131. }
  132.  
  133. /* "Monte Carlo" determination of probabilities */
  134. moncar(n,l) int n, l; {
  135. int i, j, k, b[KK], bb[KK][KK];
  136. double bf[KK], bbf[KK][KK];
  137.  
  138. nblclr();
  139. gfram(BORG,0,100);
  140. asctobin();
  141. for (k=0; k<n; k++) {
  142.   onegen(AL);
  143.   for (i=0; i<KK; i++) b[i]=0;
  144.   for (i=0; i<AL; i++) b[arr1[i]]++;
  145.   for (i=0; i<KK; i++) bf[i]=((double)(b[i]))/((double)(AL));
  146.   for (i=0; i<KK; i++) for (j=0; j<KK; j++) bb[i][j]=0;
  147.   for (i=1; i<AL; i++) bb[arr1[i-1]][arr1[i]]++;
  148.   bb[arr1[AL-1]][arr1[0]]++; 
  149.   for (i=0; i<KK; i++) for (j=0; j<KK; j++) 
  150.     bbf[i][j]=((double)(bb[i][j]))/((double)(AL));
  151.   videodot(199-(int)(100.0*bbf[1][1]),BORG+(int)(100.0*bf[1]),l);
  152.   };
  153.  videocursor(0,BROW+7,0);
  154. printf("(Monte Carlo) "); 
  155. for (i=0; i<KK; i++) printf("%2d:%5.3f ",i,bf[i]);
  156. videocursor(0,BROW+8,0);
  157. for (i=0; i<KK; i++) for (j=0; j<KK; j++) 
  158.   printf("%1d%1d:%5.3f ",i,j,bbf[i][j]);
  159. }
  160.  
  161. /* Generate coefficients of 2nd generation Bernstein Polynomial */
  162. berncoef() {
  163. int i, i0, i1, i2;
  164.  
  165. for (i=0; i<BD; i++) bp[i]=0.0;
  166. for (i0=0; i0<KK; i0++) {
  167. for (i1=0; i1<KK; i1++) {
  168. for (i2=0; i2<KK; i2++) {
  169.   if (ascrule[i0][i1][i2]=='1') bp[i0+i1+i2]+=1.0;
  170.   };};};
  171. }
  172.  
  173. /* Generate coefficients of 3rd generation Bernstein Polynomial */
  174. bernthrd() {
  175. int i, i0, i1, i2, i3, i4;
  176. int j0, j1, j2;
  177.  
  178. for (i=0; i<BD; i++) bp[i]=0.0;
  179. for (i0=0; i0<KK; i0++) {
  180. for (i1=0; i1<KK; i1++) {
  181. for (i2=0; i2<KK; i2++) {
  182. for (i3=0; i3<KK; i3++) {
  183. for (i4=0; i4<KK; i4++) {
  184.   j0=ascrule[i0][i1][i2]-'0';
  185.   j1=ascrule[i1][i2][i3]-'0';
  186.   j2=ascrule[i2][i3][i4]-'0';
  187.   if (ascrule[j0][j1][j2]=='1') bp[i0+i1+i2+i3+i4]+=1.0;
  188.   };};};};};
  189. }
  190.  
  191. /* Generate coefficients of 4th generation Bernstein Polynomial */
  192. bernfrth() {
  193. int i, i0, i1, i2, i3, i4, i5, i6;
  194. int j0, j1, j2, j3, j4;
  195. int k0, k1, k2;
  196.  
  197. for (i=0; i<BD; i++) bp[i]=0.0;
  198. for (i0=0; i0<KK; i0++) {
  199. for (i1=0; i1<KK; i1++) {
  200. for (i2=0; i2<KK; i2++) {
  201. for (i3=0; i3<KK; i3++) {
  202. for (i4=0; i4<KK; i4++) {
  203. for (i5=0; i5<KK; i5++) {
  204. for (i6=0; i6<KK; i6++) {
  205.   j0=ascrule[i0][i1][i2]-'0';
  206.   j1=ascrule[i1][i2][i3]-'0';
  207.   j2=ascrule[i2][i3][i4]-'0';
  208.   j3=ascrule[i3][i4][i5]-'0';
  209.   j4=ascrule[i4][i5][i6]-'0';
  210.   k0=ascrule[j0][j1][j2]-'0';
  211.   k1=ascrule[j1][j2][j3]-'0';
  212.   k2=ascrule[j2][j3][j4]-'0';
  213.   if (ascrule[k0][k1][k2]=='1') bp[i0+i1+i2+i3+i4+i5+i6]+=1.0;
  214.   };};};};};};};
  215. }
  216.  
  217. /* Generate coefficients of 5th generation Bernstein Polynomial */
  218. bernfifth() {
  219. int i, i0, i1, i2, i3, i4, i5, i6, i7, i8;
  220. int j0, j1, j2, j3, j4, j5, j6;
  221. int k0, k1, k2, k3, k4;
  222. int l0, l1, l2;
  223.  
  224. for (i=0; i<BD; i++) bp[i]=0.0;
  225. for (i0=0; i0<KK; i0++) {
  226. for (i1=0; i1<KK; i1++) {
  227. for (i2=0; i2<KK; i2++) {
  228. for (i3=0; i3<KK; i3++) {
  229. for (i4=0; i4<KK; i4++) {
  230. for (i5=0; i5<KK; i5++) {
  231. for (i6=0; i6<KK; i6++) {
  232. for (i7=0; i7<KK; i7++) {
  233. for (i8=0; i8<KK; i8++) {
  234.   j0=ascrule[i0][i1][i2]-'0';
  235.   j1=ascrule[i1][i2][i3]-'0';
  236.   j2=ascrule[i2][i3][i4]-'0';
  237.   j3=ascrule[i3][i4][i5]-'0';
  238.   j4=ascrule[i4][i5][i6]-'0';
  239.   j5=ascrule[i5][i6][i7]-'0';
  240.   j6=ascrule[i6][i7][i8]-'0';
  241.   k0=ascrule[j0][j1][j2]-'0';
  242.   k1=ascrule[j1][j2][j3]-'0';
  243.   k2=ascrule[j2][j3][j4]-'0';
  244.   k3=ascrule[j3][j4][j5]-'0';
  245.   k4=ascrule[j4][j5][j6]-'0';
  246.   l0=ascrule[k0][k1][k2]-'0';
  247.   l1=ascrule[k1][k2][k3]-'0';
  248.   l2=ascrule[k2][k3][k4]-'0';
  249.   if (ascrule[l0][l1][l2]=='1') bp[i0+i1+i2+i3+i4+i5+i6+i7+i8]+=1.0;
  250.   };};};};};};};};};
  251. }
  252.  
  253. /* Generate coefficients of 6th generation Bernstein Polynomial */
  254. bernsixth() {
  255. int i, i0, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10;
  256. int j0, j1, j2, j3, j4, j5, j6, j7, j8;
  257. int k0, k1, k2, k3, k4, k5, k6;
  258. int l0, l1, l2, l3, l4;
  259. int m0, m1, m2;
  260.  
  261. for (i=0; i<BD; i++) bp[i]=0.0;
  262. for (i0=0; i0<KK; i0++) {
  263. for (i1=0; i1<KK; i1++) {
  264. for (i2=0; i2<KK; i2++) {
  265. for (i3=0; i3<KK; i3++) {
  266. for (i4=0; i4<KK; i4++) {
  267. for (i5=0; i5<KK; i5++) {
  268. for (i6=0; i6<KK; i6++) {
  269. for (i7=0; i7<KK; i7++) {
  270. for (i8=0; i8<KK; i8++) {
  271. for (i9=0; i9<KK; i9++) {
  272. for (i10=0; i10<KK; i10++) {
  273.   j0=ascrule[i0][i1][i2]-'0';
  274.   j1=ascrule[i1][i2][i3]-'0';
  275.   j2=ascrule[i2][i3][i4]-'0';
  276.   j3=ascrule[i3][i4][i5]-'0';
  277.   j4=ascrule[i4][i5][i6]-'0';
  278.   j5=ascrule[i5][i6][i7]-'0';
  279.   j6=ascrule[i6][i7][i8]-'0';
  280.   j7=ascrule[i7][i8][i9]-'0';
  281.   j8=ascrule[i8][i9][i10]-'0';
  282.   k0=ascrule[j0][j1][j2]-'0';
  283.   k1=ascrule[j1][j2][j3]-'0';
  284.   k2=ascrule[j2][j3][j4]-'0';
  285.   k3=ascrule[j3][j4][j5]-'0';
  286.   k4=ascrule[j4][j5][j6]-'0';
  287.   k5=ascrule[j5][j6][j7]-'0';
  288.   k6=ascrule[j6][j7][j8]-'0';
  289.   l0=ascrule[k0][k1][k2]-'0';
  290.   l1=ascrule[k1][k2][k3]-'0';
  291.   l2=ascrule[k2][k3][k4]-'0';
  292.   l3=ascrule[k3][k4][k5]-'0';
  293.   l4=ascrule[k4][k5][k6]-'0';
  294.   m0=ascrule[l0][l1][l2]-'0';
  295.   m1=ascrule[l1][l2][l3]-'0';
  296.   m2=ascrule[l2][l3][l4]-'0';
  297.   if (ascrule[m0][m1][m2]=='1') bp[i0+i1+i2+i3+i4+i5+i6+i7+i8+i9+i10]+=1.0;
  298.   };};};};};};};};};};};
  299. }
  300.  
  301. /* evaluate the nth generation Bernstein polynomial at point p */
  302. double bern(p,n) double p; int n; {double q, s, x, r; int i, d;
  303. d=2*n-1;
  304. if (p>0.99) return bp[d];
  305. q=1.0-p; r=p/q; s=0.0; x=1.0;
  306. for (i=0; i<d;  i++) x*=q;
  307. for (i=0; i<=d; i++, x*=r) s+=bp[i]*x;
  308. s*=0.9998;
  309. return(s+0.0001);
  310. }
  311.  
  312. /* graph the probability of the second generation */
  313. pdiff(n) int n; {
  314. berncoef();
  315. gfram(CORG,0,n);
  316. gdiag(CORG,0,n);
  317. bgraf(CORG,0,2,n);
  318. }
  319.  
  320. /* graph the iterated probability of the second generation */
  321. pidiff(n) int n; {int i; double bern(), en, p;
  322. en=1.0/((double)(n));
  323. berncoef();
  324. gfram(CORG,0,n);
  325. gdiag(CORG,0,n);
  326. for (i=0; i<=n; i++) {
  327.   p=((double)(i))*en;
  328.   videodot(199-(int)(100.0*bern(bern(p,2),2)),CORG+(int)(100.0*p),3);
  329.   };
  330. }
  331.  
  332. /* graph the second iterated probability of the second generation */
  333. pjdiff(n) int n; {int i; double bern(), en, p;
  334. en=1.0/((double)(n));
  335. berncoef();
  336. gfram(CORG,0,n);
  337. gdiag(CORG,0,n);
  338. for (i=0; i<=n; i++) {
  339.   p=((double)(i))*en;
  340.   videodot(199-(int)(100.0*bern(bern(bern(p,2),2),2)),CORG+(int)(100.0*p),3);
  341.   };
  342. }
  343.  
  344. /* graph the probability of the third generation */
  345. pdiff3(n) int n; {
  346. bernthrd();
  347. gfram(CORG,0,n);
  348. gdiag(CORG,0,n);
  349. bgraf(CORG,0,3,n);
  350. }
  351.  
  352. /* graph the probability of the fourth generation */
  353. pdiff4(n) int n; {
  354. bernfrth();
  355. gfram(CORG,0,n);
  356. gdiag(CORG,0,n);
  357. bgraf(CORG,0,4,n);
  358. }
  359.  
  360. /* graph the probability of the fifth generation */
  361. pdiff5(n) int n; {
  362. bernfifth();
  363. gfram(CORG,0,n);
  364. gdiag(CORG,0,n);
  365. bgraf(CORG,0,5,n);
  366. }
  367.  
  368. /* graph the probability of the sixth generation */
  369. pdiff6(n) int n; {
  370. bernsixth();
  371. gfram(CORG,0,n);
  372. gdiag(CORG,0,n);
  373. bgraf(CORG,0,6,n);
  374. }
  375.  
  376. /* display frequencies in ascrule in color l */
  377. asfreq(l) int l; {
  378. int    i, j, i0, i1, i2;
  379. int    stat[KK], stal, stac, star;        /* statistic counts */
  380. double staf[KK], pp;
  381.  
  382. videoscroll(BROW,0,BROW+8,40,0,0);
  383. videocursor(0,BROW,0);
  384. pp=1.0/((double)(KK*KK*KK));
  385. stal=0; stac=0; star=0;
  386. for (i=0; i<KK; i++) stat[i]=0;
  387. for (i0=0; i0<KK; i0++) {
  388. for (i1=0; i1<KK; i1++) {
  389. for (i2=0; i2<KK; i2++) {
  390.   j=ascrule[i0][i1][i2]-'0';
  391.   stat[j]++;
  392.   if(j==i0) star++;
  393.   if(j==i1) stac++;
  394.   if(j==i2) stal++;
  395.   };};};
  396. for (i=0; i<KK; i++) printf("%2d    - %5.2f\n",i,((double)stat[i])*pp);
  397. printf("\n");
  398. printf("left  - %5.2f\n",((double)stal)*pp);
  399. printf("still - %5.2f\n",((double)stac)*pp);
  400. printf("right - %5.2f\n",((double)star)*pp);
  401. for (i=0; i<KK; i++) staf[i]=((double)stat[i])*pp;
  402. j=199-(int)(100.0*staf[1]);
  403. i=CORG+(int)(100.0*staf[1]);
  404. if (j>=2) videodot(j-2,i,l);
  405. if (i>=2) videodot(j,i-2,l);
  406. videodot(j,i,1);
  407. if (i<=318) videodot(j,i+2,l);
  408. if (j<=198) videodot(j+2,i,l);
  409. }
  410.  
  411. /* display state frequencies in arr1 as a point on probability graph */
  412. /* n = number of points to plot */
  413. /* l = color of plotted point   */
  414. lifreq(n,l) int n, l; {
  415. int    i, j, os, ns;
  416. double of, nf, s;
  417.  
  418. s=100.0/((double)AL);
  419. asctobin();
  420. for (i=0; i<n; i++) {
  421.   os=0; ns=0;
  422.   for (j=0; j<AL; j++) if (arr1[j]==1) os++;
  423.   onegen(AL);
  424.   for (j=0; j<AL; j++) if (arr1[j]==1) ns++;
  425.   of=(double)os;
  426.   nf=(double)ns;
  427.   videodot(199-(int)(s*nf),CORG+(int)(s*of),l);
  428.   };
  429. }
  430.  
  431. /* evaluate the parameters for two-block probabilities after one generation */
  432. twoblock(z) double *z; {
  433. int    i, i0, i1, i2, i3;
  434. int    j, j0, j1;
  435. double w, x, y, p, q, pp;
  436. double a[KK][KK], b[KK], c[KK][KK];
  437.  
  438. p=z[0]; q=1.0-p; pp=z[1];
  439. a[1][1]=pp; a[1][0]=a[0][1]=p-pp; a[0][0]=1.0-2*p+pp;
  440. b[1]=p; b[0]=q;
  441. for (i=0; i<KK; i++) {for (j=0; j<KK; j++) c[i][j]=0.0;};
  442. for (i0=0; i0<KK; i0++) {
  443. for (i1=0; i1<KK; i1++) {
  444. for (i2=0; i2<KK; i2++) {
  445. for (i3=0; i3<KK; i3++) {
  446.   j0=ascrule[i0][i1][i2]-'0';
  447.   j1=ascrule[i1][i2][i3]-'0';
  448.   x=a[i0][i1]*a[i1][i2]*a[i2][i3];
  449.   y=b[i1]*b[i2];
  450.   w=0.0; if (y>0.00001) w=x/y;
  451.   c[j0][j1]+=w;
  452.   };};};};
  453. z[0]=c[1][0]+c[1][1]; z[1]=c[1][1];
  454. }
  455.  
  456. /* graph the 2-block parameter differences to find self-consistent values */
  457. /* nxn points                                   */
  458. twoblockfreq(n) int n; {
  459. int i, j, l;
  460. double op[2], np[2], s, t, u, v;
  461.  
  462. gfram(0,0,100);
  463.  
  464. s=1.0/((double)(n));
  465. op[1]=s;
  466. for (i=1; i<n; i++) {
  467.   op[0]=s;
  468.   for (j=1; j<n; j++) {
  469.     np[0]=op[0]; np[1]=op[1]; twoblock(np);
  470.     u=np[0]-op[0]; v=np[1]-op[1];
  471.     t=u*u+v*v;
  472.     if (t<0.0025) l=0; else {
  473.     if (t<0.0675) {if((i+j)%2) l=0; else l=3;} else {
  474.     if (t<0.3906) l=3; else {
  475.     if (t<0.7500) {if((i+j)%2) l=3; else l=2;} else {
  476.     if (t<1.5000) l=2; else {
  477.     if (t<3.0000) {if((i+j)%2) l=2; else l=1;} else {
  478.     if (t<6.0000) l=1; else {
  479.                   if((i+j)%2) l=1; else l=0;
  480.        }}}}}}}
  481.     videodot(199-i,j,l);
  482.     if (kbdst()) {kbdin(); return;};
  483.     op[0]+=s;
  484.     };  /* end for-j */
  485.   op[1]+=s;
  486.   };    /* end for-i */
  487. }
  488.  
  489. /* evaluate the one-block probabilities after one generation */
  490. onebl(x,a) double x[KK], a[KK]; {int i0, i1, i2, j0;
  491.  
  492. for (j0=0; j0<KK; j0++) x[j0]=0.0;
  493.  
  494. for (i0=0; i0<KK; i0++) {
  495. for (i1=0; i1<KK; i1++) {
  496. for (i2=0; i2<KK; i2++) {
  497.   j0=ascrule[i0][i1][i2]-'0';
  498.   x[j0]+=a[i0]*a[i1]*a[i2];
  499.   };};};
  500. }
  501.  
  502. /* iterate the 1-block parameters to find self-consistent values */
  503. /* graph the iterative steps in bar-chart form             */
  504. /* ll - initial line    */
  505. /* mm - length of line    */
  506. /* nn - number of lines    */
  507. oneblfreq(ll,mm,nn) int ll, mm, nn; {
  508. int    ii, i, l, m, n;
  509. double op[KK], np[KK];
  510. double d, e, f, s;
  511.  
  512. m=0;
  513. f=(double)mm;
  514. s=1.0/((double)(KK));
  515. n=(int)(f*s);
  516. videodot(ll,m++,3);
  517. for (i=0; i<KK; i++) {op[i]=s; for (l=0; l<n; l++) videodot(ll,m++,i);};
  518.  
  519. for (ii=1; ii<=nn; ii++) {
  520.   e=0.0; m=0;
  521.   onebl(np,op);
  522.   videodot(ll+ii,m++,3);
  523.   for (i=0; i<KK; i++) {
  524.     n=(int)(f*np[i]); if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,i);
  525.     d=op[i]-np[i];
  526.     e+=d*d;
  527.     op[i]=np[i];
  528.     };
  529.   videodot(ll+ii,m++,3);
  530.   if (op[0]<=0.001) break;
  531.   if (op[1]<=0.001) break;
  532.   if (e<=0.0000001) break;
  533.   if (kbdst()) {kbdin(); break;};
  534.   };
  535. videocursor(0,BROW+7,0);
  536. printf("(1-block) "); 
  537. for (i=0; i<KK; i++) printf("%2d:%5.3f ",i,op[i]);
  538. }
  539.  
  540. /* evaluate the two-block probabilities after one generation */
  541. twobl(x,a) double x[KK][KK], a[KK][KK]; {
  542. int    i0, i1, i2, i3;
  543. int    j0, j1;
  544. double w, b[KK];
  545.  
  546. for (j0=0; j0<KK; j0++) {for (j1=0; j1<KK; j1++) x[j0][j1]=0.0;};
  547.  
  548. for (j0=0; j0<KK; j0++) {b[j0]=0.0; for (j1=0; j1<KK; j1++) b[j0]+=a[j0][j1];};
  549.  
  550. for (i0=0; i0<KK; i0++) {
  551. for (i1=0; i1<KK; i1++) {
  552. for (i2=0; i2<KK; i2++) {
  553. for (i3=0; i3<KK; i3++) {
  554.   j0=ascrule[i0][i1][i2]-'0';
  555.   j1=ascrule[i1][i2][i3]-'0';
  556.   w=a[i0][i1]*a[i1][i2]*a[i2][i3];
  557.   if (w!=0.0) w/=b[i1]*b[i2];
  558.   x[j0][j1]+=w;
  559.   };};};};
  560. }
  561.  
  562. /* iterate the 2-block parameters to find self-consistent values */
  563. /* graph the iterative steps in bar-chart form             */
  564. /* ll - initial line    */
  565. /* mm - length of line    */
  566. /* nn - number of lines    */
  567. twoblfreq(ll,mm,nn) int ll, mm, nn; {
  568. int    ii, i, j, l, m, n;
  569. double op[KK][KK], np[KK][KK];
  570. double b[KK], d, e, f, s;
  571.  
  572. m=0;
  573. f=(double)mm;
  574. s=1.0/((double)(KK*KK));
  575. n=(int)(f*s);
  576. videodot(ll,m++,3);
  577. for (i=0; i<KK; i++) {
  578. for (j=0; j<KK; j++) {
  579.   op[i][j]=s;
  580.   for (l=0; l<n; l++) videodot(ll,m++,j);
  581.   };};
  582. videodot(ll,m++,3);
  583.  
  584. for (i=0; i<KK; i++) {b[i]=0.0; for (j=0; j<KK; j++) b[i]+=op[i][j];};
  585. videodot(199-(int)(100.0*op[1][1]),(int)(100.0*b[1]),1);
  586.  
  587. for (ii=1; ii<=nn; ii++) {
  588.   e=0.0; m=0;
  589.   twobl(np,op);
  590.   videodot(ll+ii,m++,3);
  591.   for (i=0; i<KK; i++) {
  592.   for (j=0; j<KK; j++) {
  593.     n=(int)(f*np[i][j]);
  594.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,j);
  595.     d=op[i][j]-np[i][j];
  596.     e+=d<0.0?-d:d;
  597.     op[i][j]=np[i][j];
  598.     };};
  599.   videodot(ll+ii,m++,3);
  600.   for (i=0; i<KK; i++) {b[i]=0.0; for (j=0; j<KK; j++) b[i]+=op[i][j];};
  601.   videodot(199-(int)(100.0*op[1][1]),(int)(100.0*b[1]),2);
  602.   if (e<=0.0001) break;
  603.   if (kbdst()) {kbdin(); break;};
  604.   };
  605. videocursor(0,BROW+7,0);
  606. printf("(2-block) "); 
  607. for (i=0; i<KK; i++) printf("%2d:%5.3f ",i,b[i]);
  608. videocursor(0,11,0);
  609. for (i=0; i<KK; i++) for (j=0; j<KK; j++) 
  610.         printf("%1d%1d:%5.3f ",i,j,op[i][j]);
  611. videodot(199-(int)(100.0*op[1][1]),(int)(100.0*b[1]),3);
  612. videodot(199-(int)(100.0*op[1][1]),BORG+(int)(100.0*b[1]),1);
  613. }
  614.  
  615. /* evaluate the three-block probabilities after one generation */
  616. thrbl(x,a) double x[KK][KK][KK], a[KK][KK][KK]; {
  617. int    i0, i1, i2, i3, i4;
  618. int    j0, j1, j2;
  619. double w, b[KK][KK];
  620.  
  621. for (j0=0; j0<KK; j0++) {
  622. for (j1=0; j1<KK; j1++) {
  623.   for (j2=0; j2<KK; j2++) x[j0][j1][j2]=0.0;
  624.   b[j0][j1]=0.0;
  625.   for (j2=0; j2<KK; j2++) b[j0][j1]+=a[j0][j1][j2];  
  626.   };};
  627.  
  628. for (i0=0; i0<KK; i0++) {
  629. for (i1=0; i1<KK; i1++) {
  630. for (i2=0; i2<KK; i2++) {
  631. for (i3=0; i3<KK; i3++) {
  632. for (i4=0; i4<KK; i4++) {
  633.   j0=ascrule[i0][i1][i2]-'0';
  634.   j1=ascrule[i1][i2][i3]-'0';
  635.   j2=ascrule[i2][i3][i4]-'0';
  636.   w=a[i0][i1][i2]*a[i1][i2][i3]*a[i2][i3][i4];
  637.   if (w!=0.0) w/=b[i1][i2]*b[i2][i3];
  638.   x[j0][j1][j2]+=w;
  639.   };};};};};
  640. }
  641.  
  642. /* iterate the 3-block parameters to find self-consistent values */
  643. /* ll - initial line    */
  644. /* mm - length of line    */
  645. /* nn - number of lines    */
  646. thrblfreq(ll,mm,nn) int ll, mm, nn; {
  647. int    ii, i, j, k, l, m, n;
  648. double op[KK][KK][KK], np[KK][KK][KK];
  649. double b[KK], bb[KK][KK], d, e, f, s;
  650.  
  651. m=0;
  652. f=(double)mm;
  653. s=1.0/((double)(KK*KK*KK));
  654. n=(int)(f*s);
  655. videodot(ll,m++,3);
  656. for (i=0; i<KK; i++) {
  657. for (j=0; j<KK; j++) {
  658. for (k=0; k<KK; k++) {
  659.   op[i][j][k]=s;
  660.   for (l=0; l<n; l++) videodot(ll,m++,k);
  661.   };};};
  662. videodot(ll,m++,3);
  663.  
  664. for (ii=1; ii<=nn; ii++) {
  665.   e=0.0; m=0;
  666.   thrbl(np,op);
  667.   videodot(ll+ii,m++,3);
  668.   for (i=0; i<KK; i++) {
  669.   for (j=0; j<KK; j++) {
  670.   for (k=0; k<KK; k++) {
  671.     n=(int)(f*np[i][j][k]);
  672.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,k);
  673.     d=op[i][j][k]-np[i][j][k];
  674.     e+=d<0.0?-d:d;
  675.     op[i][j][k]=np[i][j][k];
  676.     };};};
  677.   videodot(ll+ii,m++,3);
  678.   if (e<=0.0001) break;
  679.   if (kbdst()) {kbdin(); break;};
  680.   };
  681.  
  682. for (i=0; i<KK; i++) {
  683.   b[i]=0.0;
  684.   for (j=0; j<KK; j++) {
  685.   for (k=0; k<KK; k++) {
  686.     b[i]+=op[i][j][k];
  687.     };}; };
  688. videocursor(0,BROW+7,0);
  689. printf("(3-block) "); 
  690. for (i=0; i<KK; i++) printf("%2d:%5.3f ",i,b[i]);
  691.  
  692. for (i=0; i<KK; i++) {
  693. for (j=0; j<KK; j++) {
  694.   bb[i][j]=0.0;
  695.   for (k=0; k<KK; k++) {
  696.     bb[i][j]+=op[i][j][k];
  697.     };}; };
  698.  
  699. videocursor(0,11,0);
  700. for (i=0; i<KK; i++) for (j=0; j<KK; j++) 
  701.   printf("%1d%1d:%5.3f ",i,j,bb[i][j]);
  702. videodot(199-(int)(100.0*bb[1][1]),BORG+(int)(100.0*b[1]),2);
  703. }
  704.  
  705. /* evaluate the four-block probabilities after one generation */
  706. foubl(x,a) double x[KK][KK][KK][KK], a[KK][KK][KK][KK]; {
  707. int    i0, i1, i2, i3, i4, i5;
  708. int    j0, j1, j2, j3;
  709. double w, b[KK][KK][KK];
  710.  
  711. for (j0=0; j0<KK; j0++) {
  712. for (j1=0; j1<KK; j1++) {
  713. for (j2=0; j2<KK; j2++) {
  714.   for (j3=0; j3<KK; j3++) x[j0][j1][j2][j3]=0.0;
  715.   b[j0][j1][j2]=0.0;
  716.   for (j3=0; j3<KK; j3++) b[j0][j1][j2]+=a[j0][j1][j2][j3];  
  717.   };};};
  718.  
  719. for (i0=0; i0<KK; i0++) {
  720. for (i1=0; i1<KK; i1++) {
  721. for (i2=0; i2<KK; i2++) {
  722. for (i3=0; i3<KK; i3++) {
  723. for (i4=0; i4<KK; i4++) {
  724. for (i5=0; i5<KK; i5++) {
  725.   j0=ascrule[i0][i1][i2]-'0';
  726.   j1=ascrule[i1][i2][i3]-'0';
  727.   j2=ascrule[i2][i3][i4]-'0';
  728.   j3=ascrule[i3][i4][i5]-'0';
  729.   w=a[i0][i1][i2][i3]*a[i1][i2][i3][i4]*a[i2][i3][i4][i5];
  730.   if (w!=0.0) w/=b[i1][i2][i3]*b[i2][i3][i4];
  731.   x[j0][j1][j2][j3]+=w;
  732.   };};};};};};
  733. }
  734.  
  735. /* iterate the 4-block parameters to find self-consistent values */
  736. /* ll - initial line    */
  737. /* mm - length of line    */
  738. /* nn - number of lines    */
  739. foublfreq(ll,mm,nn) int ll, mm, nn; {
  740. int    ii, i0, i1, i2, i3, l, m, n;
  741. double op[KK][KK][KK][KK], np[KK][KK][KK][KK];
  742. double b[KK], bb[KK][KK], d, e, f, s;
  743.  
  744. m=0;
  745. f=(double)mm;
  746. s=1.0/((double)(KK*KK*KK*KK));
  747. n=(int)(f*s);
  748. videodot(ll,m++,3);
  749. for (i0=0; i0<KK; i0++) {
  750. for (i1=0; i1<KK; i1++) {
  751. for (i2=0; i2<KK; i2++) {
  752. for (i3=0; i3<KK; i3++) {
  753.   op[i0][i1][i2][i3]=s;
  754.   for (l=0; l<n; l++) videodot(ll,m++,i3);
  755.   };};};};
  756. videodot(ll,m++,3);
  757.  
  758. for (ii=1; ii<=nn; ii++) {
  759.   e=0.0; m=0;
  760.   foubl(np,op);
  761.   videodot(ll+ii,m++,3);
  762.   for (i0=0; i0<KK; i0++) {
  763.   for (i1=0; i1<KK; i1++) {
  764.   for (i2=0; i2<KK; i2++) {
  765.   for (i3=0; i3<KK; i3++) {
  766.     n=(int)(f*np[i0][i1][i2][i3]);
  767.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,i3);
  768.     d=op[i0][i1][i2][i3]-np[i0][i1][i2][i3];
  769.     e+=d<0.0?-d:d;
  770.     op[i0][i1][i2][i3]=np[i0][i1][i2][i3];
  771.     };};};};
  772.   videodot(ll+ii,m++,3);
  773.   if (e<=0.0001) break;
  774.   if (kbdst()) {kbdin(); break;};
  775.   };
  776.  
  777. for (i0=0; i0<KK; i0++) {
  778.   b[i0]=0.0;
  779.   for (i1=0; i1<KK; i1++) {
  780.   for (i2=0; i2<KK; i2++) {
  781.   for (i3=0; i3<KK; i3++) {
  782.     b[i0]+=op[i0][i1][i2][i3];
  783.     };};}; };
  784.  
  785. videocursor(0,BROW+7,0);
  786. printf("(4-block) "); 
  787. for (i0=0; i0<KK; i0++) printf("%2d:%5.3f ",i0,b[i0]);
  788.  
  789. for (i0=0; i0<KK; i0++) {
  790. for (i1=0; i1<KK; i1++) {
  791.   bb[i0][i1]=0.0;
  792.   for (i2=0; i2<KK; i2++) {
  793.   for (i3=0; i3<KK; i3++) {
  794.     bb[i0][i1]+=op[i0][i1][i2][i3];
  795.     };}; };};
  796.  
  797. videocursor(0,11,0);
  798. for (i0=0; i0<KK; i0++) for (i1=0; i1<KK; i1++) 
  799.   printf("%1d%1d:%5.3f ",i0,i1,bb[i0][i1]);
  800. videodot(199-(int)(100.0*bb[1][1]),BORG+(int)(100.0*b[1]),2);
  801. }
  802.  
  803. /* evaluate the five-block probabilities after one generation */
  804. fivbl(x,a) double x[KK][KK][KK][KK][KK], a[KK][KK][KK][KK][KK]; {
  805. int    i0, i1, i2, i3, i4, i5, i6;
  806. int    j0, j1, j2, j3, j4;
  807. double w, b[KK][KK][KK][KK];
  808.  
  809. for (j0=0; j0<KK; j0++) {
  810. for (j1=0; j1<KK; j1++) {
  811. for (j2=0; j2<KK; j2++) {
  812. for (j3=0; j3<KK; j3++) {
  813.   for (j4=0; j4<KK; j4++) x[j0][j1][j2][j3][j4]=0.0;
  814.   b[j0][j1][j2][j3]=0.0;
  815.   for (j4=0; j4<KK; j4++) b[j0][j1][j2][j3]+=a[j0][j1][j2][j3][j4];  
  816.   };};};};
  817.  
  818. for (i0=0; i0<KK; i0++) {
  819. for (i1=0; i1<KK; i1++) {
  820. for (i2=0; i2<KK; i2++) {
  821. for (i3=0; i3<KK; i3++) {
  822. for (i4=0; i4<KK; i4++) {
  823. for (i5=0; i5<KK; i5++) {
  824. for (i6=0; i6<KK; i6++) {
  825.   j0=ascrule[i0][i1][i2]-'0';
  826.   j1=ascrule[i1][i2][i3]-'0';
  827.   j2=ascrule[i2][i3][i4]-'0';
  828.   j3=ascrule[i3][i4][i5]-'0';
  829.   j4=ascrule[i4][i5][i6]-'0';
  830.   w=a[i0][i1][i2][i3][i4]*a[i1][i2][i3][i4][i5]*a[i2][i3][i4][i5][i6];
  831.   if (w!=0.0) w/=b[i1][i2][i3][i4]*b[i2][i3][i4][i5];
  832.   x[j0][j1][j2][j3][j4]+=w;
  833.   };};};};};};};
  834. }
  835.  
  836. /* iterate the 5-block parameters to find self-consistent values */
  837. /* ll - initial line    */
  838. /* mm - length of line    */
  839. /* nn - number of lines    */
  840. fivblfreq(ll,mm,nn) int ll, mm, nn; {
  841. int    ii, i0, i1, i2, i3, i4, l, m, n;
  842. double op[KK][KK][KK][KK][KK], np[KK][KK][KK][KK][KK];
  843. double b[KK], bb[KK][KK], d, e, f, s;
  844.  
  845. m=0;
  846. f=(double)mm;
  847. s=1.0/((double)(KK*KK*KK*KK*KK));
  848. n=(int)(f*s);
  849. videodot(ll,m++,3);
  850. for (i0=0; i0<KK; i0++) {
  851. for (i1=0; i1<KK; i1++) {
  852. for (i2=0; i2<KK; i2++) {
  853. for (i3=0; i3<KK; i3++) {
  854. for (i4=0; i4<KK; i4++) {
  855.   op[i0][i1][i2][i3][i4]=s;
  856.   for (l=0; l<n; l++) videodot(ll,m++,i4);
  857.   };};};};};
  858. videodot(ll,m++,3);
  859.  
  860. for (ii=1; ii<=nn; ii++) {
  861.   e=0.0; m=0;
  862.   fivbl(np,op);
  863.   videodot(ll+ii,m++,3);
  864.   for (i0=0; i0<KK; i0++) {
  865.   for (i1=0; i1<KK; i1++) {
  866.   for (i2=0; i2<KK; i2++) {
  867.   for (i3=0; i3<KK; i3++) {
  868.   for (i4=0; i4<KK; i4++) {
  869.     n=(int)(f*np[i0][i1][i2][i3][i4]);
  870.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,i4);
  871.     d=op[i0][i1][i2][i3][i4]-np[i0][i1][i2][i3][i4];
  872.     e+=d<0.0?-d:d;
  873.     op[i0][i1][i2][i3][i4]=np[i0][i1][i2][i3][i4];
  874.     };};};};};
  875.   videodot(ll+ii,m++,3);
  876.   if (e<=0.0001) break;
  877.   if (kbdst()) {kbdin(); break;};
  878.   };
  879. for (i0=0; i0<KK; i0++) {
  880.   b[i0]=0.0;
  881.   for (i1=0; i1<KK; i1++) {
  882.   for (i2=0; i2<KK; i2++) {
  883.   for (i3=0; i3<KK; i3++) {
  884.   for (i4=0; i4<KK; i4++) {
  885.     b[i0]+=op[i0][i1][i2][i3][i4];
  886.     };};};}; };
  887. videocursor(0,BROW+7,0);
  888. printf("(5-block) "); 
  889. for (i0=0; i0<KK; i0++) printf("%2d:%5.3f ",i0,b[i0]);
  890.  
  891. for (i0=0; i0<KK; i0++) {
  892. for (i1=0; i1<KK; i1++) {
  893.   bb[i0][i1]=0.0;
  894.   for (i2=0; i2<KK; i2++) {
  895.   for (i3=0; i3<KK; i3++) {
  896.   for (i4=0; i4<KK; i4++) {
  897.     bb[i0][i1]+=op[i0][i1][i2][i3][i4];
  898.     };}; };};};
  899.  
  900. videocursor(0,11,0);
  901. for (i0=0; i0<KK; i0++) for (i1=0; i1<KK; i1++) 
  902.   printf("%1d%1d:%5.3f ",i0,i1,bb[i0][i1]);
  903. videodot(199-(int)(100.0*bb[1][1]),BORG+(int)(100.0*b[1]),2);
  904. }
  905.  
  906. /* evaluate the six-block probabilities after one generation */
  907. sixbl(x,a) double x[KK][KK][KK][KK][KK][KK], a[KK][KK][KK][KK][KK][KK]; {
  908. int    i0, i1, i2, i3, i4, i5, i6, i7;
  909. int    j0, j1, j2, j3, j4, j5;
  910. double w, b[KK][KK][KK][KK][KK];
  911.  
  912. for (j0=0; j0<KK; j0++) {
  913. for (j1=0; j1<KK; j1++) {
  914. for (j2=0; j2<KK; j2++) {
  915. for (j3=0; j3<KK; j3++) {
  916. for (j4=0; j4<KK; j4++) {
  917.   for (j5=0; j5<KK; j5++) x[j0][j1][j2][j3][j4][j5]=0.0;
  918.   b[j0][j1][j2][j3][j4]=0.0;
  919.   for (j5=0; j5<KK; j5++) b[j0][j1][j2][j3][j4]+=a[j0][j1][j2][j3][j4][j5];  
  920.   };};};};};
  921.  
  922. for (i0=0; i0<KK; i0++) {
  923. for (i1=0; i1<KK; i1++) {
  924. for (i2=0; i2<KK; i2++) {
  925. for (i3=0; i3<KK; i3++) {
  926. for (i4=0; i4<KK; i4++) {
  927. for (i5=0; i5<KK; i5++) {
  928. for (i6=0; i6<KK; i6++) {
  929. for (i7=0; i7<KK; i7++) {
  930.   j0=ascrule[i0][i1][i2]-'0';
  931.   j1=ascrule[i1][i2][i3]-'0';
  932.   j2=ascrule[i2][i3][i4]-'0';
  933.   j3=ascrule[i3][i4][i5]-'0';
  934.   j4=ascrule[i4][i5][i6]-'0';
  935.   j5=ascrule[i5][i6][i7]-'0';
  936.   w=a[i0][i1][i2][i3][i4][i5]*a[i1][i2][i3][i4][i5][i6]*a[i2][i3][i4][i5][i6][i7];
  937.   if (w!=0.0) w/=b[i1][i2][i3][i4][i5]*b[i2][i3][i4][i5][i6];
  938.   x[j0][j1][j2][j3][j4][j5]+=w;
  939.   };};};};};};};};
  940. }
  941.  
  942. /* iterate the 6-block parameters to find self-consistent values */
  943. /* ll - initial line    */
  944. /* mm - length of line    */
  945. /* nn - number of lines    */
  946. sixblfreq(ll,mm,nn) int ll, mm, nn; {
  947. int    ii, i0, i1, i2, i3, i4, i5, l, m, n;
  948. double op[KK][KK][KK][KK][KK][KK];
  949. double np[KK][KK][KK][KK][KK][KK];
  950. double b[KK], bb[KK][KK], d, e, f, s;
  951.  
  952. m=0;
  953. f=(double)mm;
  954. s=1.0/((double)(KK*KK*KK*KK*KK*KK));
  955. n=(int)(f*s);
  956. videodot(ll,m++,3);
  957. for (i0=0; i0<KK; i0++) {
  958. for (i1=0; i1<KK; i1++) {
  959. for (i2=0; i2<KK; i2++) {
  960. for (i3=0; i3<KK; i3++) {
  961. for (i4=0; i4<KK; i4++) {
  962. for (i5=0; i5<KK; i5++) {
  963.   op[i0][i1][i2][i3][i4][i5]=s;
  964.   for (l=0; l<n; l++) videodot(ll,m++,i5);
  965.   };};};};};};
  966. videodot(ll,m++,3);
  967.  
  968. for (ii=1; ii<=nn; ii++) {
  969.   e=0.0; m=0;
  970.   sixbl(np,op);
  971.   videodot(ll+ii,m++,3);
  972.   for (i0=0; i0<KK; i0++) {
  973.   for (i1=0; i1<KK; i1++) {
  974.   for (i2=0; i2<KK; i2++) {
  975.   for (i3=0; i3<KK; i3++) {
  976.   for (i4=0; i4<KK; i4++) {
  977.   for (i5=0; i5<KK; i5++) {
  978.     n=(int)(f*np[i0][i1][i2][i3][i4][i5]);
  979.     if (n>0) for (l=0; l<n; l++) videodot(ll+ii,m++,i5);
  980.     d=op[i0][i1][i2][i3][i4][i5]-np[i0][i1][i2][i3][i4][i5];
  981.     e+=d<0.0?-d:d;
  982.     op[i0][i1][i2][i3][i4][i5]=np[i0][i1][i2][i3][i4][i5];
  983.     };};};};};};
  984.   videodot(ll+ii,m++,3);
  985.   if (e<=0.0001) break;
  986.   if (kbdst()) {kbdin(); break;};
  987.   };
  988. for (i0=0; i0<KK; i0++) {
  989.   b[i0]=0.0;
  990.   for (i1=0; i1<KK; i1++) {
  991.   for (i2=0; i2<KK; i2++) {
  992.   for (i3=0; i3<KK; i3++) {
  993.   for (i4=0; i4<KK; i4++) {
  994.   for (i5=0; i5<KK; i5++) {
  995.     b[i0]+=op[i0][i1][i2][i3][i4][i5];
  996.     };};};};}; };
  997. videocursor(0,BROW+7,0);
  998. printf("(6-block) "); 
  999. for (i0=0; i0<KK; i0++) printf("%2d:%5.3f ",i0,b[i0]);
  1000.  
  1001. for (i0=0; i0<KK; i0++) {
  1002. for (i1=0; i1<KK; i1++) {
  1003.   bb[i0][i1]=0.0;
  1004.   for (i2=0; i2<KK; i2++) {
  1005.   for (i3=0; i3<KK; i3++) {
  1006.   for (i4=0; i4<KK; i4++) {
  1007.   for (i5=0; i5<KK; i5++) {
  1008.     bb[i0][i1]+=op[i0][i1][i2][i3][i4][i5];
  1009.     };}; };};};};
  1010.  
  1011. videocursor(0,11,0);
  1012. for (i0=0; i0<KK; i0++) for (i1=0; i1<KK; i1++) 
  1013.   printf("%1d%1d:%5.3f ",i0,i1,bb[i0][i1]);
  1014. videodot(199-(int)(100.0*bb[1][1]),BORG+(int)(100.0*b[1]),2);
  1015. }
  1016.  
  1017. /* end */
  1018.