/*   APRIL 4 1991 8pm  *   */ 
/*  rndifs20_c  based on Algorithm March 1991 page 2 
         Random Fractal Iterated Functions                                  */ 

#include <stdio_h> 
int a[21], b[21], c[21], d[21], e[21], f[21], p[21], ps[21];  
int pk[3], fpw[3], fpx[3], fpy[3], fpz[3], sqrt2[3], hfsqrt2[3];  int rnd[3],
ra[3], rd[3], rm[3], rr[3], det[3], tmp[3], coeff[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[];  
int i, j, k, n, t;  

main()   { 
   prompt("choose number");  
   choose_n();  
   initrand();  
   prompt("select random number seed");  
   pickrand();  
   coefpute();  
   prompt("see array values");  
   seearray();  
   prompt("determine raw probabilities");  
   rawprob();  
   prompt("equalize probabilities (sum = 1)");  
   sumprob();  
   prompt("get ranges of graphics values");  
   sizeup();  
   trace("Scaling factors based on 100 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");  
   prompt("see IFS develop, about 5 minutes");  
   forshow(x,y);  
   at(ow,0,0);  
   printf("any key to exit\n");  
   getchar();  
   exit();  
} 

seearray()  { 
   csize(ow,0,0);  
   cls(ow);  printf("%7c%14c%14c%14c%14c%14c\n",'a','b','c','d','e','f');  
   for ( i = 0; i < n; ++i ) { 
          fpica(a,i,tmp);  
          at(ow,( 3 + (i*2)),(0*14));  printf("%13s\n",ftoa(tmp,str));  
          fpica(b,i,tmp);  
          at(ow,( 3 + (i*2)),(1*14));  printf("%13s\n",ftoa(tmp,str));  
          fpica(c,i,tmp);  
          at(ow,( 3 + (i*2)),(2*14));  printf("%13s\n",ftoa(tmp,str));  
          fpica(d,i,tmp);  
          at(ow,( 3 + (i*2)),(3*14));  printf("%13s\n",ftoa(tmp,str));  
          fpica(e,i,tmp);  
          at(ow,( 3 + (i*2)),(4*14));  printf("%13s\n",ftoa(tmp,str));  
          fpica(f,i,tmp);  
          at(ow,( 3 + (i*2)),(5*14));  
          printf("%13s\n",ftoa(tmp,str));  
   } 
   csize(ow,1,0);  
} 

forshow(x,y) int x[3], y[3];  { 
   scale(ow,pscale,xoffset,yoffset);  
   ink(ow,4);  
   cls(ow);  
   for ( i = 0 ; i < 3072 ; i++) { 
       iterate(x,y);  
       point(ow,x,y);  
   } 
} 

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()   { 
int k;  
   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(-6);  
   *maxx = float(-8);  
   *minx = float(8);  
   *miny = float(6);  
   cls(stdout);  
   trace("LOOPING INTO 100 ITERATIONS");  
   for(i=0; i<100; 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,30);  
       printf("%d  \n",i);  
       at(stdout,10,10);  
       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.1");  
   *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");  
   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));  fpu-
       ta(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( ftoa(det,str) == "0" ) { 
          *det = atof(".1");  
          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[21], element, value[3];  { int i;  
   element *= 3;  
   for(i = 0; i < 3 ; i++) { 
       array[i+element] = value[i];  
   } 
} 

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

rand()   { 
   int t, rd[3], rr[3], rm[3];  
       t = 125;  
       *rm = float(t);  
       *rd = atof("2796203");  
       *ra = fmult(ra,rm);  
       *ra = fdiv(ra,rd);  
       t = int(ra);  
       *rr = float(t);  
       *rnd = fsub(ra,rr);  
       *ra = fmult(rnd,rd);  
       return rnd ;  
} 

coef()   { 
   *coeff = fmult(sqrt2,rand());  
   *coeff = fsub(coeff,hfsqrt2);  
   return coeff;  
} 

choose_n()  { 
   cls(stdout);  
   n = 0;  
   while( n<2 || n >8 ) { 
       trace("Choose a number between 2 and 6\n");  fgets(ch,1,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();  
} 

coefpute()  { 
   for ( i = 0; i < n; ++i ) { 
       coef();  
       fputa(a,i,coeff);  
       coef();  
       fputa(b,i,coeff);  
       coef();  
       fputa(c,i,coeff);  
       coef();  
       fputa(d,i,coeff);  
       coef();  
       fputa(e,i,coeff);  
       coef();  
       fputa(f,i,coeff);  
       coef();  
       fputa(p,i,coeff);  
   } 
   trace("Floating point Arrays have been loaded");  
} 

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

initrand()  { 
   *ra = atof("100001");  
   *hfsqrt2 = atof("2");  
   *sqrt2 = sqrt(hfsqrt2);  
   *hfsqrt2 = fdiv(sqrt2,hfsqrt2);  
} 

pickrand()  { 
   cls(stdout);  
   printf("\n\n\n\nTouch   ENTER  to get new random number seeds\n");
   printf("\n\n\n\nTouch   [ SPACE BAR ]  to select seed shown\n");  rand();  
       at (stdout,15,20);  
       ftoa(rnd,str);  
       printf("%s        \n",str);  
   while(getchar() != ' '){ 
       rand();  
       at (stdout,15,20);  
       ftoa(rnd,str);  
       printf("%s        \n",str);  
   } 
   cls(stdout);  
} 

_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);  
}