Program preim;
uses Graph,Crt;
label 0,1,2,3,4;
const
  a=6;         {alphabet size}
  d=2;         {diameter = 2*radius}
  nprod=24;    {number of productions}
  nforb=8;     {number of forbidden words}
  PROD:array[0..nprod-1] of string =  {local rule}
       ('x100','x110','x203','x213','34x3','44x3','35x0','45x0',
        'x142','14x5','1240','1543','1zz2','yy45','12x1','x544',
        'x1x0','x4x3','xyz2','yzx5','x022','53x5','z2x1','x5y4');
  FORB:array[0..nforb-1] of string =     {forbidden words}
       ('10','11','20','21','34','44','35','45');
var X,Y:array[-d..1000] of byte;  {list of images and preimages}
    iforb,         {pointer to forbidden word}
    lx,ly,         {length of preimage and image}
    nx,ny,         {number of preimages in X and images in Y}
    numx,numy,     {number of preimages found and to be found}
    len,i,j,l,n,j0,j1:longint;
    out:text;


function rule(i:integer):byte;   {local rule}
label 1,2,3;
var iprod,j,j0,j1:integer;
    st:string;
begin
  rule:=X[i+1];
  for iprod:=0 to nprod-1 do begin
      for j:=0 to d do begin
          st:=copy(PROD[iprod],j+1,1);
          if st='x' then goto 1;
          if (st='y') then if (X[i+j]<3) then goto 1 else goto 2;
          if (st='z') then if (X[i+j]>2) then goto 1 else goto 2;
          val(st,j0,j1);
          if j0=X[i+j] then goto 1;
          goto 2;
        1:end;
      val(copy(PROD[iprod],d+2,1),j0,j1); rule:=j0;
      goto 3;
    2:end;
3:end;

begin
 assign(out,'preim.txt');
 rewrite(out);
 writeln(out,'Cellular automaton with non-sofic mixing limit set');
 write(out,'local rule:');
 for i:=0 to nprod-1 do begin
     if (i mod 8)=0 then begin writeln(out); write(out,'   '); end;
     write(out,copy(PROD[i],1,3),':',copy(PROD[i],4,1),', ');
     end;
 writeln(out);
 writeln(out,'   x in (0,1,2,3,4,5) y in (0,1,2), z in (3,4,5)');
 writeln(out);
 write(out,'forbidden words:');
 for i:=0 to nforb-1 do write(out,' ',FORB[i],',');
 writeln(out); writeln(out);
 for iforb:=0 to nforb-1 do begin
     ly:=2; numy:=1;
     for i:=0 to ly-1 do begin
         val(copy(FORB[iforb],i+1,1),j0,j1); Y[i]:=j0;
         end;
   0:nx:=0; ny:=0; lx:=ly+d;
   1:write(out,'preimages of ');
     for i:=0 to ly-1 do write(out,Y[ny*ly+i]); write(out,':');
     numx:=0;
     for j:=0 to d do X[nx*lx+j]:=0;
     i:=0;
   2:if (rule(nx*lx+i)<>Y[ny*ly+i]) and (i>=0) then begin
      3:if X[nx*lx+i+d]<a-1 then begin
           X[nx*lx+i+d]:=X[nx*lx+i+d]+1;
           goto 2;
           end;
        if i>-d then begin i:=i-1; goto 3; end;
        goto 4;
        end;
     if (rule(nx*lx+i)=Y[ny*ly+i]) or (i<0) then begin
        if i<ly-1 then begin
           i:=i+1; X[nx*lx+i+d]:=0; goto 2; end;
        numx:=numx+1;
        if numx>1 then write(out,',');
        if (numx mod 7=1) then begin
           writeln(out); write(out,'   ');  end;
        for j:=0 to lx-1 do begin
            write(out,X[nx*lx+j]); X[nx*lx+lx+j]:=X[nx*lx+j]; end;
        nx:=nx+1;
        goto 3;
        end;
   4:ny:=ny+1;
     if numx=0 then write(out,' none');
     writeln(out);
     if ny<numy then goto 1;
     for n:=0 to nx-1 do begin
         for l:=0 to lx-1 do Y[n*lx+l]:=X[n*lx+l]
         end;
     numy:=nx;
     if numy>0 then begin ly:=lx; goto 0; end;
     end;
  close(out);
  end.