/*  rndifs25_c  based on Algorithm March 1991 page 2 
modified to do 20000 pixels for checkerboard 
32 squares * 3 = 96 , make 99 
                  Random Fractal -  Iterated Functions Systems               */


#include <stdio_h> 
int a[99], b[99], c[99], d[99], e[99], f[99], p[99], ps[99];  
int pk[3], fpw[3], fpx[3], fpy[3], fpz[3];  
int rnd[3], ra[3], rd[3], rm[3], rr[3], det[3], tmp[3];  
int w[3], x[3], y[3], z[3], xs[3], ys[3], maxx[3], maxy[3], minx[3], miny[3];
int yscale[3], xscale[3], pscale[3], xoffset[3], yoffset[3];  
char *ow, ch, str[16];  
int i, j, k, n, t, iter;  

main()   { 
   n = 32;  
   initrand();  
   prompt("enter array values");  
   selarray();  
   prompt("determine raw probabilities");  
   rawprob();  
   prompt("equalize probabilities, sum = 1");  
   sumprob();  
   prompt("get ranges of graphics values");  
   sizeup();  
   prompt("see IFS develop, 10 minutes for 10K points");  forshow(x,y);  
   at(ow,0,0);  
   printf("any key to exit");  
   getchar();  
   exit();  
} 

getsit() { 
       printf("?");  
       gets(str);  
       *tmp = atof(str);  
       return tmp;  
} 

selarray()  { 
   csize(ow,0,0);  
   cls(ow);  
   *fpz = atof("8");  
   for ( j = 0; j < n; ++j ) { 
       *tmp = atof(".125");  
       fputa(a,j,tmp);  
       fputa(d,j,tmp);  
       *tmp = atof("0");  
       fputa(b,j,tmp);  
       fputa(c,j,tmp);  
       k = (2*(j%4)+(j/4)%2)%8;  
       *tmp=float(k);  
       *tmp=fdiv(tmp,fpz);  
       printf("value for e @ %u =",j);  
       fputa(e,j,tmp);  
       fpica(e,j,tmp);  
       printf("%s      ",ftoa(tmp,str));  
       k = (j/4);  
       *tmp=float(k);  
       *tmp=fdiv(tmp,fpz);  
       printf("  value for f @ %u =",j);  
       fputa(f,j,tmp);  
       fpica(f,j,tmp);  
       printf("%s      \n",ftoa(tmp,str));  
   } 
   csize(ow,1,0);  
} 

askmore()  { 
   at(ow,24,0);  
   printf("%u0K want another 10K ?",iter);  
   ch = getchar();  
   if(ch != 'n')  { 
       at(ow,24,0);  
       printf("                          ");  
       do10k();  
   } 
} 

forshow(x,y) int x[3], y[3];  { 
   scale(ow,pscale,xoffset,yoffset);  
   ink(ow,6);  
   cls(ow);  
   point(ow,maxx,maxy);  
   point(ow,maxx,miny);  
   point(ow,minx,maxy);  
   point(ow,minx,miny);  
   ink(ow,4);  
   iter = 0;  
   do10k();  
} 

do10k() { 
   for ( i = 0 ; i < 10240 ; i++) { 
       iterate(x,y);  
       point(ow,x,y);
   } 
   iter ++;  
   askmore();  
} 

iterate(x,y) int x[3], y[3];   { 
   k = selectk();  
   fpica(a,k,fpw);  
   fpica(b,k,fpx);  
   fpica(c,k,fpy);  
   fpica(d,k,fpz);  
   *xs = fmult(fpw,x);  
   *tmp = fmult(fpx,y);  
   *xs = fadd(xs,tmp);  
   fpica(e,k,fpw);  
   *xs = fadd(xs,fpw);  
   *ys = fmult(fpy,x);  
   *tmp = fmult(fpz,y);  
   *ys = fadd(ys,tmp);  
   fpica(f,k,fpw);  
   *ys = fadd(ys,fpw);  
   fmove(x,xs);  
   fmove(y,ys);  
   return (x,y);  
} 

selectk()  { 
   fmove(pk,rand());  
   for(k=0; k<n; k++)   { 
       fpica(ps,k,fpw);  
       *fpz = fsub(pk,fpw);  
       if(int(fpz) < 0)   break;  
   } 
   return (k);  
} 

sizeup() { 
   trace("BEGINNING OF SIZEUP");  
   *x = float(0);  
   *y = float(0);  
   *xs = float(0);  
   *ys = float(0);  
   *maxy = float(-1000);  
   *maxx = float(-1000);  
   *minx = float(1000);  
   *miny = float(1000);  
   cls(stdout);  
   trace("LOOPING INTO 200 ITERATIONS");  
   for(i=0; i<200; i++) { 
       iterate(x,y);  
       if (fcmp(y,maxy)  >0)  {  fmove(maxy,y); } 
       if (fcmp(x,maxx)  >0)  {  fmove(maxx,x); } 
       if (fcmp(y,miny)  <0)  {  fmove(miny,y); } 
       if (fcmp(x,minx)  <0)  {  fmove(minx,x); } 
       at(stdout,10,31);  
       printf("%u  \n",i);  
       at(stdout,10,8);  
       printf("X min = %s   ",  ftoa(minx,str));  
       at(stdout,10,36);  
       printf("X max = %s   \n",ftoa(maxx,str));  
       at(stdout,16,25);  
       printf("Y min = %s   ",  ftoa(miny,str));  
       at(stdout,4,25);  
       printf("Y max = %s   \n",ftoa(maxy,str));  
   } 
   trace("BEGINNING GRAPHICS SCALING");  
   *w = atof("1.2"); /* 10% margin  */ 
   *tmp = fsub(maxy,miny);  
   *tmp = fabs(tmp);  
   *yscale = fmult(tmp,w);  
   *tmp = fsub(maxx,minx);  
   *tmp = fabs(tmp);  
   *xscale = fmult(tmp,w);  
   if(fcmp(xscale,yscale) >0 )
      fmove(pscale,xscale);  
   else
      fmove(pscale,yscale);
   fmove(xoffset,minx);  
   fmove(yoffset,miny);  
   *fpx = fsub(maxx,minx);  
   *fpx = fabs(fpx);  
   *fpy = fsub(maxy,miny);  
   *fpy = fabs(fpy);  
   *w = atof("1.3");  
   *fpw = fmult(pscale,w);  
   *fpw = fsub(fpw,fpx);  
   *fpz = fsub(pscale,fpy);  
   *w = atof(".5");  
   *fpw = fmult(w,fpw);  
   *fpz = fmult(w,fpz);  
   *xoffset = fsub(minx,fpw);  
   *yoffset = fsub(miny,fpz);  
   trace("END OF SCALING\n\n");  
   trace("Scaling factors based on 200 iterations:");
   printf("pscale =  %s\n",ftoa(pscale,str));
   printf("xoffset = %s  \n",ftoa(xoffset,str));  
   printf("yoffset = %s  \n",ftoa(yoffset,str));  
   printf("x = %s  \n",ftoa(x,str));  
   printf("y = %s  \n",ftoa(y,str));
   scale(stdout,pscale,xoffset,yoffset);  
   trace("Ready to do a IFS\n");  
   return(pscale,xoffset,yoffset);  
} 

sumprob()   { 
   z[0] = float(0);  
   for( i = 0; i < n; i++) { 
       fpica(p,i,fpw);  
       *tmp = fdiv(fpw,pk);  
       *z = fadd(z,tmp);  
       printf("cumulative probability at %u is %s \n",i,ftoa(z,str));
       fputa(ps,i,z);
       fputa(p,i,tmp);  
   } 
   printf("sum of probabilities = %s\n",ftoa(z,str));  
} 

rawprob()   { 
   *pk = float(0);  
   for( i=0; i<n; i++)  { 
       fpica(a,i,fpw);  
       fpica(d,i,fpx);  
       fpica(b,i,fpy);  
       fpica(c,i,fpz);  
       *det = fmult(fpw,fpx);  
       *tmp = fmult(fpy,fpz);  
       *det = fsub(det,tmp);  
       *det = fabs(det);  
       if( *det == float(0) ) { 
          *det = atof("+.01");  
          trace("ZERO DETERMINANT !");  
       } 
       *pk = fadd(pk,det);  
       fputa(p,i,det);  
       printf("det(%u) = %s\n",i,ftoa(det,str));  
   } 
   printf("Sum of raw probabilities = %s \n",ftoa(pk,str));  
} 

fputa(array,element,value) int  array[99], element, value[3];  { int i;  
   element *= 3;  
   for(i = 0; i < 3 ; i++) { 
       array[i+element] = value[i];  
   } 
} 

fpica(array,element,value) int  array[99], element, value[3];  { int i;  
   element *= 3;  
   for(i = 0; i < 3 ; i++) { 
       value[i] = array[i+element];  
   } 
} 

rand()   { 
       *ra = fmult(ra,rm);  
       *ra = fdiv(ra,rd);  
       t = int(ra);  
       *rr = float(t);  
       *rnd = fsub(ra,rr);  
       *ra = fmult(rnd,rd);  
       return rnd ;  
} 

choose_n()  { 
   cls(stdout);  
   n = 0;  
   while( n<2 || n >33 ) { 
       trace("Choose a number between 2 and 33\n");
       fgets(ch,3,stdin);  
       n = atoi(ch);  
       printf("%2u was chosen \n",n);  
   } 
   return(n);  
} 

prompt(str) char str[];  { 
   printf("\nTouch [ SPACE BAR ] to %s\n",str);  
   getchar();  
} 

trace(str) char str[];   { 
   printf("\ntrace %s\n",str);  
} 

initrand()  { 
   *rm = float(125);  
   *rd = atof("2796203");  
   *ra = atof("100001");  
   for( i = 0; i<99; i++)  { 
       a[i]=0;  
       b[i]=0;  
       c[i]=0;  
       d[i]=0;  
       e[i]=0;  
       f[i]=0;  
       p[i]=0;  
       ps[i]=0;  
   } 
} 

_console() { 
   mode(4);  
   fopen(ow,"con","w");  
   window(ow,512,256,0,0);  
   paper(ow,0);  
   ink(ow,4);  
   border(ow,1,2);  
   cls(ow);  
}