/*This is the computer program Tie version 2.0 called tie_v2_0.txt. (Turing machine Inference Engine)

 Copyright (C) 2018 John Nixon
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    A copy of the GNU General Public License can be obtained from
    <http://www.gnu.org/licenses/>.

This program is based on the translation into D of TIEv1.0 (written in C++). It generates the set of all
irreducible regular rules (IRR) for any Turing Machine (TM) up to any specified
length, written by John Nixon (john.h.nixon1@gmail.com formerly john.h.nixon@spamcop.net).
Comments are welcome.
In the terminology of paper 1, "rule" is to be read as
"regular rule" i.e. RR. Likewise "irreducible_rule" and "reducible_rule" are
to be read as IRR and RRR respectively as defined in paper 1.
The latest version the program and draft papers are available from www.bluesky-home.co.uk under Turing 
machine analysis. The first paper describing the theory behind it and what the
output means has been published (paper 1):

Mathematica Aeterna, (2013) Vol. 3, no. 9, 709 - 738 (Paper 1)

second paper continuing this study has now been published

Mathematica Aeterna, (2017) Vol. 7, no. 1, 19 - 56 (Paper 2)

Program modifications

New in version 2.0:
The procedure for generating the IRR has been changed to include generating the origins of the IRR i.e. the IRR are now
represented as triplets. This gives a more informative result providing the proof of the
reachability of the LHS's of the IRR and the algorithm has been changed as a result and
is now much easier to describe. It also runs about twice as fast as version 1.3 which is partly due to avoiding the
unnecessary copying of IRR through the default method D has for copying objects (pointer
copies) instead of physical copies obtained with .dup. I had concerns about doing this at
first until I realised that it is quite easy to check the program for the absence of  "back doors"
i.e. other variables that point to the same memory that could be used to change such a variable. There
will often be such other variables initially and it must be checked that they go out of scope without modifying
the variable of interest. This mostly deals with my earlier comment about .dup.

New in version 1.3: the IRR are printed with the additional "L" or "R" giving the
pointer position on the LHS of the IRR. This is just a convenience because the pointer
position is included as before, but the 'type' of the rule e.g. LL, LR, RL or RR can now be
read off quickly after the second CS e.g. to determine if the rule has type RL or LR which allows
it to be extended by an additional symbol and not be reducible.

New in version 1.2:
(1) The program has been translated from C++ into D (dlang.org).
[D is much easier to use than C++, has bounds checking making anything like Valgrind
unnecessary as far as I can see, has many more features and has similar syntax.
This otherwise excellent language does unfortunately require the programmer to set up
and use functions called ``dup" that duplicate things that I think should be done 
automatically in assignments and copying objects into arrays etc. If the programmer
forgets to do this correctly a new pointer will be assigned to the object instead
of making a new object copy, so if the object is subsequently changed, errors could
arise in code that appears naively to be correct.]
(2) The symbols are now in the alphabet {a,b,c,...z} and the states are numbered 1,2,3,...
allowing for Turing Machines with many more states and symbols to be analysed.
This also removes ambiguity in symbolic expressions where exponents are used and is in accord
with the notation in Paper 2.
(3) A much improved format for reading Turing machines.
(4) There was a modification to the way halting is handled: halting is now considered
to be a property of a TM instruction rather than a TM state. As such, the TM can either move right (R)
left(L) by one square, or halt (H) at each TM step. This avoids needing a distinguished state, but it means that
a configuration set (CS) reached as a result of a halt must specify this explicitly e.g. by putting (H) after it.
In any case my interest is mainly in describing non-terminating computations so I probably won't
make much use of this in future.

The IRR are essentially the non-redundant set of computational short-cuts that
can be derived by running the Turing Machine(TM) on a finite length of tape.
Every computation of the TM as a sequence of single TM steps can be replaced
by an equivalent shorter sequence of applications of the short-cut rules,
(though the set of such rules may be and usually is infinite). Any computation with a TM
expressed in terms of its IRR is reduced in terms of the number of computation steps needed.

By correctly guessing the iterative formulae or method for obtaining the IRR for a particular TM
with the aid of the list of IRR generated by this program,
and verifying them by an induction argument, it is possible to obtain long-
term information on the behaviour of the TM that is not obvious initially, and
may shed light on the problem which the TM was designed to solve. Because TM's
can be set up to solve any problem that can be solved by an algorithm, the
possibilities of application of this method seem limitless.

The notation for each irreducible regular rule (IRR), including the original
Turing Machine (TM) table, is in the following order: two configuration sets (CS)'s
each defined by a quadruple (state,pointer position,length,tape symbols) representing
the input and output respectively of the rule, followed by the status of the rule 
(L left-moving,R right-moving,H halt), and finally the numbers of steps involved in its
derivation (number of derivation steps from rules of length n-1, number of TM steps).
The lines of the TM machine table must be sorted first by old state increasing, and
then by symbol read in alphabetical order a to z.

This program must be compiled before it can be used. It was written in D and
originally had extension .d, but since webservers can misinterpret such files
I prefer to keep it on the web as a .txt file and rename it before compilation. To compile
it go to dlang.org and download one of the compilers (I use dmd).
then do for example 
mv tie_v2_0.txt tie.d
dmd tie.d
and run it with the command that includes the file name of the TM to be analysed.
./tie <TM_file_name>.txt
*/

import std.stdio;
import std.algorithm;
import std.string;
import std.file;
import std.conv;
int nsy;//# symbols
int nst;//# machine states
// Global variable: The set of IRR for the TM, used by main and apply_rules
IRR [][] irreducible_rules;

//********************************** function can_write **************************************
//This function finds a filename to be written to for output and returns it. It ensures
//that no file is overwritten without the consent of the user. "base_name" is usually the input
//datafile name and "extension" is the default file extension for output (including the dot)
//arrays and class objects are not copied
char[] can_write(const char[] base_name,const char[] extension){
  char[] filename = (base_name ~ extension).dup;
  char[] chars;
  char[] question;
  bool write1,fn,overwrite=false;
  do{
    fn = exists(filename);
    if(fn){
      writeln("The file ",filename," already exists.");
      question="Do you want me to overwrite it?".dup;    
      overwrite=yes(question);
      if(!overwrite){
        writeln("The new path and file name is ",base_name,"<extra characters>",extension);
        write("Enter the extra characters: ");
        readln(chars);chars=chomp(chars);
        filename=base_name ~ chars ~ extension;}}
    write1=(!fn) || overwrite;}
  while (!write1);
  writeln("Output will be to the file ",filename);
  return filename;}

//************************************** function yes *****************************************
//This function asks the user the question by dispaying the text 'question' and waits for
//the response y (for yes) or n (for no)
bool yes(const char[] question){
  char ans='*';
  bool ret;
  char[] res;
  while(ans!='y' && ans!='n'){
    std.stdio.write(question,"(y/n): ");
    res=chomp(readln().dup);
    ans=res[0];
    if(res.length==1){
      if(ans=='y')ret = true;
      else if(ans=='n') ret = false;}
    else{ans='*';}}
  return ret;
}

//*********************************** input_int ***************************************
//The purpose of this function is to extract an int from the input stream and handle
//the error that occurs when a non-number is entered, and allow the user to enter it again.
void input_int(ref int i){
  bool error;
  char[] line;
  for(;;){
    line = chomp(readln().dup);
    error=false;
    try{i = to!int(line);} catch (Exception){error=true;}
    if(error == false)break;
    write("An error occurred. Please try again: ");}
}
//************************** struct CS ************************************************
//This class represents a configuration set(CS) as defined in the paper.
//I will use the convention that if p>0, symbols will be added to the left
//of the pointer, and if p=0, symbols will be added to the right.
//This distinguishes the two cases when l=0 because then p=1 and 0 respectively.
struct CS{
  CS dup()const{
    CS cs;
    cs.s = this.s;
    cs.p = this.p;
    cs.l = this.l;
    cs.t = this.t.dup;
    return cs;}
  int s;//state
  int p,l;//pointer position,length of tape given
  char[] t;//part of tape of specified length
}

//******************************** function add_at_x *******************************
//Add symbol sy at the opposite end of the string of symbols from the pointer.
//The routine assumes that the pointer is at one end of the known
//symbols on the tape (other cases will not be needed).
void add_at_x(ref CS cs,char sy){//symbol is added at right. Pointer is at left. ( i.e. cs.p=1)
  ++(cs.l);
  cs.t~=sy;
  if(cs.p==cs.l-1){//the symbol is added at the left, pointer is at right (cs.p=cs.l)
    ++(cs.p);
    rotate(cs.t,1);}}
//********************************* function add_at_ptr ****************************
//Add symbol sy to cs at the pointer position.
//The routine assumes that the pointer is just off one end of the known
//symbols on the tape.
void add_at_ptr(ref CS cs,char sy){
  ++(cs.l);
  cs.t~=sy;
  if(!cs.p){
    cs.p=1;
    rotate(cs.t,1);
  }}
//********************************* function add_at_end ****************************
//Add symbol sy to cs at the end of the string closest to the pointer position
//This assumes that the pointer is at one end of the string.
void add_at_end(ref CS cs,char sy){
  ++(cs.l);
  cs.t~=sy;
  if(cs.p==1){
    cs.p=2;
    rotate(cs.t,1);
  }}
//********************************* function add_at_left ****************************
//Add symbol sy to cs at the end of the string closest to the pointer position
//This assumes that the pointer is at one end of the string.
void add_at_left(ref CS cs,char sy){
  ++(cs.l);
  ++(cs.p);
  cs.t~=sy;
  rotate(cs.t,1);
 }
//********************************* function add_at_right ****************************
//Add symbol sy to cs at the end of the string closest to the pointer position
//This assumes that the pointer is at one end of the string.
void add_at_right(ref CS cs,char sy){
    ++(cs.l);
    cs.t~=sy;}
//******************************** function rotate *********************************
//Rotates a char[] array of length l by r places to the right. 

void rotate(char[] str,int r){
  char[] temp = str.dup;
  int l=to!int(str.length);
  temp.length = 2*l-r;
  foreach(i;0 .. l-r)temp[i+l] = temp[i];
  foreach(i;0 .. l)str[i] = temp[i+l-r];}
//******************************* function subset **********************************
//Returns whether cs1 is a subset of cs2 (includes equality)
bool subset(const ref CS cs1,const ref CS cs2){
  int ll;
  if(cs1.s!=cs2.s)return 0;//false if states differ
  ll=cs1.p-cs2.p;//difference between left end of tape represented
  if(ll<0)return 0;
  if((cs2.l-cs2.p)>(cs1.l-cs1.p))return 0;
  char[] sub;
  sub.length=cs2.l;
  foreach(i;0 .. cs2.l)sub[i]=cs1.t[i+ll];
  return sub==cs2.t;}
//******************************** struct IRR **************************************
//relation between two CS's c1 and c2 defined by the action of the TM. In the paper
//special types of rule were defined: regular rules i.e. rules defined by a simple 
//recursive definition, and the subset of these that are non-redundant i.e. derived
//in more than one step. These the the irreducible regular rules (IRR). Later the
//array of CS's c0 was added to the definition. These are all the CS's that lead by
//one or more TM steps to c1.
struct IRR{
  CS[] c0;
  CS c1,c2;
  char lsig;//pointer position in c1: 'L' for 1 and 'R' for c1.l, blank if c1.l=1
  char sig;//L,R,H, or C i.e. left moving, right moving, halting, or cycling resp.
//See the paper for details.
  int n_derivation;//# derivation steps
  int n_machine;//# machine steps
  IRR dup(){
    IRR r;
    r.c0 = this.c0.dup;
    r.c1 = this.c1.dup;
    r.c2 = this.c2.dup;
    r.lsig = this.lsig;
    r.sig = this.sig;
    r.n_derivation = this.n_derivation;
    r.n_machine = this.n_machine;
    return r;}
}

//****************************** function compare **********************************
//Order relation for sorting rules: Compare first by machine states of the LHS's. If
//these are the same compare by pointer positions of the LHS's, if these are equal
//and the pointer is at 1, compare by symbol strings on the "tape", and if the pointer
//is not at 1, compare by the reversed symbol strings on the "tape".
bool compare(const ref IRR rule1,const ref IRR rule2){
  if(rule1.c1.s!=rule2.c1.s){
    return rule1.c1.s< rule2.c1.s;}
  else if(rule1.c1.p != rule2.c1.p){
    return rule1.c1.p < rule2.c1.p;}
  else if(rule1.c1.p==1){
    return rule1.c1.t < rule2.c1.t;}
  else{
    char[] str1=rule1.c1.t.dup;
    char[] str2=rule2.c1.t.dup;
    reverse(str1);
    reverse(str2);
    return str1<str2;}}

//************************* function new_origins ************************************
//New_origins starts with the CS cs and finds all single backward TM steps from there and their
//resulting CS's. This is repeated for each branch recursively, until
//one of 4 conditions occurs in each branch. Condition 1: the pointer finishes at the opposite end from 
//where it started in cs. Condition 2: the pointer finishes at the same end to where it started in cs.
//Condition 3: the computation can go no further because backward steps from there are impossible.
//Finally, condition 4: a CS is reached that has been entered before.
//Of course the computation on a branch halts in any of these cases. All the CS's arrived at as a result of
//branches ending in condition 1 constitute the array returned.
//
//The algorithm is based on a depth-first search of the tree. An array of the cs's obtained on
//the way (L) must be maintained in order to check for a loop (condition 4). This has to be updated at each step of
//the search. The first time new_origins is called for a CS cs, L is [cs]. Otherwise L is the current list of all
//previous CS's in the backward search for origins, excluding the one that has just been or being found.
//The new element of L is added before each call to new_origins.
//
//The CS cs is initally obtained by adding a symbol at one end or the other (the opposite end to where
//the pointer is) in IRR.c1 which is the same end to where the pointer is for the first call to
//new_origins. The initial pointer position in IRR.c1 is passed down through the variable
//"right_at_start" meaning the pointer is at the right in IRR.c1. This is used to distinguish
//case 1 and case 2 when the pointer ends up at the opposite end or same end
//respectively.
void new_origins(ref CS cs,bool right_at_start,ref CS[] L,ref CS[] ori){
  bool cycle=canFind(L,cs);
  bool right_at_end = cs.p == cs.l;
  bool left_at_end = cs.p == 1;
  if(cycle || right_at_end || left_at_end){//conditions under which new_origins must stop.
    if(!cycle && ((right_at_start && left_at_end) || (!right_at_start && right_at_end)))ori~=cs;//update ori. no dup because cs is not altered
//before it goes out of scope.
    if(L.length)--L.length;
}//remove last item in L if there are any.
  else{
//get array(one_step_ori) of all cs2 that go to cs in 1 TM step
    CS[] one_step_ori;
    int j;
    int ind;
    foreach(k1; 0 .. nsy){
      foreach(k2; 1 .. nst+1){
        j=(k2-1)*nsy+k1;
        if(irreducible_rules[0][j].c2.s==cs.s){
          ind = cs.p-irreducible_rules[0][j].c2.p;//index of symbol matched (cs.p-2 if TM goes right,cs.p if TM goes left)
          if(ind>=0 && ind < cs.l){
            if(irreducible_rules[0][j].c2.t[0]==cs.t[ind]){
              CS cs2;//cs2 is obtained by a backward TM step from cs.
              cs2.s=irreducible_rules[0][j].c1.s;
              if(irreducible_rules[0][j].sig=='L')cs2.p=cs.p+1;else if(irreducible_rules[0][j].sig=='R')cs2.p=cs.p-1;
              cs2.l=cs.l;
              cs2.t=cs.t.dup;//without the .dup, the next line will change cs.t that changes cs2 that could change one_step_ori in the
//calling foreach loop below.(check).
              cs2.t[ind]=irreducible_rules[0][j].c1.t[0];
              one_step_ori~=cs2;//dup not needed: cs2 goes out of scope next, so it cannot provide a "back door" to one_step_ori
}}}}}
//This is the recursive bit
    foreach(cs2;one_step_ori){
      L~=cs;
      new_origins(cs2,right_at_start,L,ori);}}}

//************************* function apply_rules ************************************
//The function "apply_rules" applies the previously found IRR up to length j to carry out substitution
//steps for the TM computation starting from the LHS named cs.
//Do this repeatedly until one of the following
//results occurs: (1)pointer moves outside the length of tape with known symbols,
//to the right (R) or left(L).
//(2) a halt condition occurs (H), (3) an infinitely repeating cycle occurs(C).
//The result is a new rule with the stopping condition (LRHC) indicated and
//the numbers of derivation (i.e. substitution) and machine steps involved. In the case C there
//is no final CS. In the output there must be one, but C is
//indicated showing that it is not really final.
//The new rule is added to IRR, because it must be irreducible
//as the result of the correct functioning of this computer program.
//Each computation will end in state L,R,H or C, meaning the pointer is
//at the left,right, the computation halts or does not terminate, respectively.
//On output, cs now is the final result of the substitutions. The new rules get
//inserted into the array of arrays of IRR "irreducible_rules" in the main program.

void apply_rules(ref CS cs,ref int n_steps,ref int nm,ref char d,ref char lsig){//cslist is working space, d is "sig" L,R,C, or H
  if(cs.p==1)lsig='L'; else lsig='R';
  n_steps=0;
  CS[] cs_array;
  cs_array~=cs.dup;//dup is needed because cs_array is a record of the changing values of cs.
  n_steps=0;//# derivation steps
  nm=0;//# machine steps
  bool match;
  int j=cs.l;
  d='0';
  IRR irr2;
  for(;;){//repeatedly apply appropriate rules to continue computation
    if(cs.p==0){d='L';break;}
    if(cs.p>cs.l){d='R';break;}
    match=0;
    for(int k=j-2;k>=0;--k){//find the longest applicable rule
      foreach(irr3;irreducible_rules[k]){
        if(subset(cs,irr3.c1)){match=1;
          irr2 = irr3;//irr3 next goes out of scope so dup not needed provided irr2 is not altered bolow.
          break;}}
      if(match==1)break;}
//match now indicates if a matching rule has been found. This should be true always
    if(!match){
      throw new Exception("Error: no matching rule found");}
//Do substitution. The first position is element 0 as far as replace is concerned
    if(cs.p==1){//pointer at left in cs and irr2.
      foreach(i;0 .. irr2.c2.l){
        cs.t[i]=irr2.c2.t[i];}}
    else if(cs.p>1){//pointer at right in cs and irr2.
      foreach(i;0 .. irr2.c2.l){
        cs.t[j-irr2.c2.l+i]=irr2.c2.t[i];}}
    cs.s=irr2.c2.s;
    cs.p+=irr2.c2.p-irr2.c1.p;
    cs.l=j;//cs is now the result of applying rule irr2 to previous cs
    ++n_steps;
    nm+=irr2.n_machine;
    if(count(cs_array[],cs)>0 || irr2.sig=='C')d='C';
    cs_array~=cs.dup;//dup needed
    if(irr2.sig=='H'){d='H';break;}
    if(d=='C')break;//if this CS was found before, there is an endless stationary loop so go to next case.
    }}
//********************************* main program **********************************************
void main(string[] args){
  writeln("This is TIEv2.0 copyright (C) 2018 John Nixon.\nThis program comes with ABSOLUTELY NO WARRANTY;\nThis is free software, and you are welcome
  to redistribute it\nunder certain conditions; for more information see the source code file.\n");
  if(args.length!=2){
    writeln("Usage: <directory and program name>  <input file for Turing Machine>");
    writeln("For example from the same directory as the program:");
    writeln(" ./<executable program name> <input file for Turing Machine>");
    writeln("See the initial comments in the program source file for more information.");
    writeln("The information in the file of the Turing Machine must be presented in the following");
    writeln("order with each element separated from the next by any number (including zero) of");
    writeln("whitespace characters (tabs spaces newlines etc).");
    writeln("to allow great flexibility of presentation:");
    writeln("Number of TM states (nst), Number of TM symbols(nsy), then for each combination of these:");
    writeln("Old state,symbol read,new state,symbol printed,direction of movement (L or R) or halt (H).");
    writeln("The states are 1,2,... nst, and the symbols are a,b,c,...");
    writeln("The lines of the TM machine table must be sorted first by old state increasing, and then by");
    writeln("symbol read in alphabetical order a to z.");
    writeln("This last requirement simplifies the algorithm for searching for reverse TM steps");
    writeln("For example\n3 \n2 \n1a2bR\n1b1bH\n2a3bR\n2b3aL\n3a1aL\n3b1aR\n");
    writeln("For verification, the program's interpretation of the user specified TM is first shown.");
    return;}
  string inputfilename=args[1];
  File infile = File(inputfilename,"r");
  if(!exists(inputfilename)){writeln("I cannot find file ",inputfilename);return;}
  char[] outputfilename=can_write(inputfilename,"_out.txt");
  File of;//output file;
  of.open(outputfilename.idup,"w");
  infile.readf(" %s",&nst);infile.readf(" %s",&nsy);
  int c=nsy*nst;
  irreducible_rules.length = 1;
  write("See the documentation in the program source file for the notation.\n");
  of.write("See the documentation in the program source tie.d file for the notation.\n");
  writeln("The Turing Machine table expressed as IRR of length 1");
  IRR r1;
  foreach(oc;0 .. c){
    infile.readf(" %d",&r1.c1.s);
    r1.c1.l=1;
    r1.c1.t.length=1;
    infile.readf(" %s",&r1.c1.t[0]);
    r1.c1.p=1;
    infile.readf(" %d",&r1.c2.s);
    r1.c2.l=1;
    r1.c2.t.length=1;
    infile.readf(" %s",&r1.c2.t[0]);
    infile.readf(" %s",&r1.sig);
    if(r1.sig=='L')r1.c2.p=0;else if(r1.sig=='R')r1.c2.p=2;else r1.c2.p=1;
    r1.n_derivation=1;r1.n_machine=1;
    r1.lsig=' ';
    irreducible_rules[0]~=r1.dup;//dup is needed here because there is only one r1
    writeln(r1);}
  write("\nEnter the maximum length of derived computation rules: ");
  int n;
  input_int(n);
  write("\nEnter 1 for showing all IRR, 2 for showing only IRR of type LL, or 3 for showing only IRR of type RR: ");
  int opt;
  input_int(opt);
  irreducible_rules.length = n;
  write("\nThe numbers of irreducible regular (computation) rules");
  of.write("\nThe numbers of irreducible regular (computation) rules");
  write("\n(IRR) derived from this Turing Machine for different tape lengths are as follows:");
  of.write("\n(IRR) derived from this Turing Machine for different tape lengths are as follows:\n");
  writeln("\nSee the output file "~outputfilename~" for the list of IRR obtained.");
  int n_machine;//# machine steps
  int n_derivation;//# derivation steps
  char lsig;//pointer position in c1: 'L' for 1 and 'R' for c1.l, blank if c1.l=1
  char sig;//L,R,H, or C i.e. left moving, right moving, halting, or cycling resp.
  int count;
  CS csi;
  CS[][CS] rtm;//reversed Turing Machine steps represented as an associative array. This structure is filled in next.
//Do this first for the TM steps going left
  int index;
foreach(st;1 .. nst+1)foreach(sy;0 .. nsy){
    CS cs1;
    cs1.s=st;
    cs1.p=0;
    cs1.l=1;
    cs1.t.length=1;
    cs1.t[0]=to!char(sy+97);
    rtm[cs1]=[];
    foreach(st1;1 .. nst+1)foreach(sy1;0 .. nsy){
      index=(st1-1)*nsy+sy1;
      if(irreducible_rules[0][index].c2.p==0 && irreducible_rules[0][index].c2.s == st && irreducible_rules[0][index].c2.t[0] == to!char(sy+97)){
        CS cs2;
        cs2.s=irreducible_rules[0][index].c1.s;
        cs2.p=1;
        cs2.l=1;
        cs2.t.length=1;
        cs2.t[0]=irreducible_rules[0][index].c1.t[0];
        rtm[cs1]~=cs2;}}}//no dup needed because cs2 goes out of scope next.
//then for the TM steps going right
foreach(st;1 .. nst+1)foreach(sy;0 .. nsy){
    CS cs1;
    cs1.s=st;
    cs1.p=2;
    cs1.l=1;
    cs1.t.length=1;
    cs1.t[0]=to!char(sy+97);
    rtm[cs1]=[];
    foreach(st1;1 .. nst+1)foreach(sy1;0 .. nsy){
      index=(st1-1)*nsy+sy1;
      if(irreducible_rules[0][index].c2.p==2 && irreducible_rules[0][index].c2.s == st && irreducible_rules[0][index].c2.t[0] == to!char(sy+97)){
        CS cs2;
        cs2.s=irreducible_rules[0][index].c1.s;
        cs2.p=1;
        cs2.l=1;
        cs2.t.length=1;
        cs2.t[0]=irreducible_rules[0][index].c1.t[0];
        rtm[cs1]~=cs2;}}}//no dup needed
//generate the IRR(2)
//add symbols at the pointer s.t. forward TM step takes the pointer to opposite end of tape (l=2).
//alpha_n is the integer representation of the symbol called alpha in the papers: 0 for a ,1 for b etc. to give an irreducible rule
  foreach(element;rtm.byKeyValue){
    if(element.value==[])continue;
    foreach(alpha_n;0 .. nsy){
      index=(element.key.s-1)*nsy+alpha_n;//use ordering of TM to compute index
      bool left = element.key.p==0 && irreducible_rules[0][index].c2.p==2;//add symbol at left
      bool right = element.key.p==2 && irreducible_rules[0][index].c2.p==0;//add symbol at right
      if(left || right){
        IRR r;
        foreach(cs0;element.value){
//construct the IRR triplet of length 2 and add it to the IRR
          CS cs1=cs0.dup;//cs0 is not modified by following statements
          if(left)add_at_left(cs1,to!char(alpha_n+97));
          if(right)add_at_right(cs1,to!char(alpha_n+97));
          r.c0~=cs1;}//ok because cs1 then goes out of scope.
//r.c1=key with alpha added on left or right
        CS cs=element.key.dup;//this protects element.key from changes via cs
        add_at_ptr(cs,to!char(alpha_n+97));
        r.c1=cs;//ok because cs is not changed before it goes out of scope
        CS cs2=cs.dup;
        apply_rules(cs2,n_derivation,n_machine,sig,lsig);
        r.n_derivation=n_derivation;
        r.n_machine=n_machine;
        r.sig=sig;
        if(r.c1.p==1)r.lsig='L';else r.lsig='R';
        r.c2=cs2;//ok
        irreducible_rules[1]~=r;
}}}
//write IRR(2)
  writeln("length: ",2," # IRR = ",irreducible_rules[1].length);
  of.writeln("length: ",2," # IRR = ",irreducible_rules[1].length);
  CS[] L;//array of CS's in the path. kept to check for cycling in new_origins. Reused in every call.
  CS[] ori;//list of new origins on output.
  CS[] oris;//combined array of origins.
  bool right_at_start;
//carry out analysis for the required number of iterations i.e.
  foreach(i;2 .. n){//given all irreducible rules (IRR) of length i, find all extensions of them to length i+1.
    count=0;//initialise count of the number of IRR found of length i.
    foreach(r2;irreducible_rules[i-1]){
      if(!((r2.lsig=='L' && r2.sig=='R') || (r2.lsig=='R' && r2.sig=='L')))continue;//ignore IRR that are not of type LR or RL
//for each origin in r2, find all alpha values and the set of new origins when this alpha is added, giving a matrix of bool meaning
//the symbol alpha is used, and the array of new origin(s) if the symbol alpha is added.
//Find all the new origins, pooling the results for each previous origin in r2.
//then get new IRR. This involves getting the RHS by apply_rules.
//writeln("r2 = ",r2);
      foreach(i1;0 .. nsy){
        oris=[];
        foreach(cs;r2.c0){
          csi=cs.dup;
          add_at_end(csi,to!char(i1+97));
          right_at_start = cs.p == 1;//i.e pointer is at left in the cs that new origins is first called with.
	  L=[];
          ori=[];
          new_origins(csi,right_at_start,L,ori);
          oris~=ori;}
        if(!oris.length)continue;
        IRR r;
        r.c0=oris;//nothing changes oris before it goes out of scope so r.c0 will not be altered, so dup is not needed
        r.c1=r2.c1.dup;//r2.c1 must not be altered by change to r.c1 so dup is needed here.
        add_at_x(r.c1,to!char(i1+97));
        r.c2=r2.c2.dup;//likewise dup is needed here.
        add_at_ptr(r.c2,to!char(i1+97));
        apply_rules(r.c2,n_derivation,n_machine,sig,lsig);
        r.n_derivation=n_derivation+1;//the 1 came from the original rule
        if(n_derivation>i+1)writeln("n_derivation = ", n_derivation,"r2.c2 = ",r2.c2);
        r.n_machine=n_machine+r2.n_machine;
        r.sig=sig;
        if(r.c1.p==1)r.lsig='L';else r.lsig='R';
        irreducible_rules[i]~=r;}}//needs no dup here 
//Write out results
    writeln("length: ",i+1," # IRR = ",irreducible_rules[i].length);
    of.writeln("length: ",i+1," # IRR = ",irreducible_rules[i].length);}
  if(opt==1)of.write("All IRR are shown\n");
  if(opt==2)of.write("Only IRR of type LL are shown\n");
  if(opt==3)of.write("Only IRR of type RR are shown\n");
  foreach(i; 0 .. n){
    if(irreducible_rules[i].length!=0){
      if(i)of.writeln("IRR of length ",i+1);
      else of.writeln("The Turing Machine table expressed as IRR of length 1");
      sort!(compare)(irreducible_rules[i]);
      foreach(it;irreducible_rules[i]){
        if((opt==1)||(opt==2 && it.c1.p==1 && it.c2.p==0)||(opt==3 && it.c1.p==it.c1.l && it.c2.p==it.c2.l+1))of.writeln(it);}}}
}//modified to only print IRR of type LL or RR or all