program commutators;
(* version for filtered lists of Hall elements;
   differs only in file extensions ('f' added) *)
(* input file: hall_bas.N (Hall generators as commutators in A,B)
   -- file provided by 3rd party (necklaces.pas?);
   the union of hall_bas.N for n<=12 is contained in hall1-12.txt  *)
(* output file: bas_comm.N (Hall generators expanded over the 
   monomials in A,B and written as vectors *)
(* the same Hall generators as polynomials can be written on standard
   output *)
(* Global array VarPres initialized once and forever   ;
   function "Store" disabled *)
(* version of 1.11.03: operators T and R added;
   T is swapping A<->B; R is reversal of monomials;
   results written on "bas_comm.Nt" and "bas_comm.Nr", resp. *)
const
    MaxNumVar = 8200;    (* maximal number of variables         *)
    namelength = 12;   (* number of characters in identifiers *)
    nulname = '';
type
    nametype =  string[12];
    string = string[255];
    varsettype= array[1..MaxNumVar] of nametype;
    reltype = (succeeds,precedes,similar,cancels);
    rational = record
                 num: longint;
                 den: longint;
               end;
    exptype = array[1..MaxNumVar] of shortint;
    monom = record
              coef: rational;
              exp: exptype;
            end;
    polynom = ^node;
    node = record
                term: monom;
                next: polynom
              end;
    listnamedpol = ^lnnode;
    lnnode = record
                 name: nametype;
                 pol: polynom;
                 next: listnamedpol;
             end;
    opertype = (unplus,unminus,postarplus,postarmin,
                biplus,biminus,aster,slash,commut,arrow,
                leftpar,rightpar,nonop,beginop);
    lextype = (oper,ident,number,unary,binary);
    lexem = record
              case which: lextype of
                oper: (whato: opertype);
                ident: (whats: nametype);
                number: (whatn: longint);
              end;
    lexlist = ^lexnode;
    lexnode = record
                item: lexem;
                next: lexlist
              end;
    expr = ^enode;
    enode = record
               case tag: lextype of
                  binary: (biop: opertype;
                          arg1, arg2: expr);
                  unary:  (unop: opertype;
                          arg: expr);
                  number: (num: longint);
                  ident:  (nam: nametype);
                end;
    vector = array[1..maxnumvar] of integer;
var
    NumVar: integer; (* number of currently present variables *)
    ExtVars: listnamedpol;    (* definitions for external variables *)
    VarPres: varsettype; (* variables involved at any time (the rest *)
                         (* of the array filled with blanks)         *)
    delimiters : set of char ;
    letters : set of char ;
    digits : set of char ;
    p,p1: polynom; 
    s,s1: string;
    n: integer; (* degree of the monomials *)
    i: integer; (* index *)
    infile,outfile,outft,outfr: text;
    last: boolean;
    offset: integer;
    v: vector;
    ext: string;
    ch: char;
    count: integer;

procedure Error (msg: string);
begin
    writeln;
    writeln('*** ERROR  ',msg);
    readln; halt
end;


function SumNum (a,b: longint): longint;
var s: longint;
begin
    s:= a+b;
    if ( (a>0) and (b>0) and (s<0) )
       or ( (a<0) and (b<0) and (s>0) )
    then Error  (' Integer addition overflow ');
    sumnum:= s
end;

function ProNum (a,b: longint): longint;
begin
   pronum:= a*b;
end;

function GCD (a,b: longint): longint;
begin
   if a=0 then gcd:= b
   else if b=0 then gcd:= a
   else if a<b then gcd:= gcd (a, b mod a)
        else gcd:= gcd (a mod b, b);
end;

procedure DivideByGCD (var a,b: longint);
var g: longint;
begin
    if b<0 then begin a:=-a; b:= -b end;
    g:= gcd (abs(a), abs(b));
    a:= a div g;
    b:= b div g;
end;

procedure AddRat (a,b: rational; var s: rational);
var g,a1,b1: longint;
begin
    if (a.den=1) and (b.den=1) then begin
                                      s.num:= sumnum(a.num,b.num);
                                      s.den:=1;
                                    end
    else begin
           g:= GCD(a.den,b.den);
           a1:= a.den div g;
           b1:= b.den div g;
           s.num:= sumnum(pronum(a.num,b1), pronum(b.num,a1));
           s.den:= pronum(a.den,b1);
           DivideByGCD (s.num,s.den);
         end;
end;

procedure MulRat (a,b: rational; var p: rational);
begin
    if (a.den=1) and (b.den=1) then begin
                                      p.num:= pronum(a.num,b.num);
                                      p.den:=1
                                    end
    else begin
           DivideByGCD (a.num,b.den);
           DivideByGCD (b.num,a.den);
           p.num:= pronum(a.num,b.num);
           p.den:= pronum(a.den,b.den);
         end;
end;

procedure InvRat (r: rational; var i: rational);
begin
    if r.num = 0 then Error('inverse of 0')
    else if r.num > 0 then begin i.den:= r.num; i.num:= r.den end
    else begin i.den:= -r.num; i.num:= -r.den end
end;

function digit(ch:char):integer;
begin
   if ch='A' then digit:=0
   else if ch='B' then digit:=1
   else begin writeln('wrong character'); halt end;
end; 

function power(a,b:integer):integer;
var i,r:integer;
begin
   r:=1;  
   for i:=1 to b do r:=r*a;
   power:=r
end;

function nvar (v:nametype): integer;
   (* ad hoc hack for the list of variables AAA, AAB, ABB etc *)
var i,l,r: integer;
begin
   l:=length(v);
   r:=0;
   for i:=1 to l do r:=r*2+digit(v[i]);
   r:=r+power(2,l)-1;
   nvar:=r;
{   nvar:=numvar+1-r (* to comply with the reversion done in InitVarsAB *)}
end;

function IsConstant (m:monom) : boolean;
var i: integer;
begin
    IsConstant:= true;
    if m.coef.num <> 0 then
    for i:= 1 to NumVar do
        if m.exp[i] <> 0 then IsConstant:= false
end;

function divides (t1,t2: exptype; var ratio: exptype): Boolean;
var i: integer;
begin
    divides:= true;
    for i:= 1 to numvar do
    begin
       ratio[i]:= t2[i]-t1[i];
       if ratio[i]<0 then divides:= false;
    end;
end;

function Relation (m,n: monom): reltype;
label 21;
var
    i:integer; s: rational;
begin
    i:=1;
    while (i<=NumVar) do
       if m.exp[i]=n.exp[i] then i:=i+1 else goto 21 ;
21: if i > NumVar then begin
                          AddRat (m.coef,n.coef,s);
                          if S.num = 0 then relation:= cancels
                          else relation:= similar
                        end
    else if m.exp[i] > n.exp[i] then relation:= precedes
         else relation:= succeeds
end;

function CopyPol (p: polynom): polynom;
var q, last: polynom;
begin
    if p=nil then q:= nil else
    begin
        new(q);
        with q^ do
        begin
           term.coef:= p^.term.coef;
           term.exp:= p^.term.exp;
           next:= nil;
        end;
        last:= q;
        p:= p^.next;
        while p<>nil do
        begin
           new(last^.next);
           last:=last^.next;
           with last^ do
           begin
                 term.coef:= p^.term.coef;
                 term.exp:= p^.term.exp;
                 next:=nil
           end;
           p:= p^.next;
        end;
    end;
    CopyPol:= q
end;

procedure FreePol (var p: polynom);
var q: polynom;
begin
    while p<>nil do
    begin
       q:=p;
       p:=p^.next;
       dispose(q)
    end
end;

function ConstPol (n,d: longint): polynom;  (* rational constant n/d *)
var
   p: polynom; i: integer;
begin
   if n=0 then ConstPol:= nil else
   begin
      new (p);
      p^.next:= nil;
      p^.term.coef.num:= n;
      p^.term.coef.den:= d;
      for i:= 1 to
      MaxNumVar do p^.term.exp[i]:= 0;
      ConstPol:= p;
   end;
end;

function OneTerm (m: monom): polynom;
                 (* polynomial consisting of one monomial m *)
var
   p: polynom;
begin
   new (p);
   p^.next:= nil;
   p^.term:= m;
   OneTerm:= p
end;

function Extract (np: nametype; lp: listnamedpol): polynom;
label 21;
begin
    while (lp<>nil) do if (lp^.name=np)
        then goto 21 else lp:= lp^.next;
21: if lp=nil then Error('no definition for ' + np)
       else extract:= CopyPol(lp^.pol)
end;

function VarPol (v: nametype): polynom;
var
    m: monom;
    numv,i: integer;
begin
    numv:= nvar(v);
    if numv<>0 then
    begin
       m.coef.num:= 1;
       m.coef.den:= 1;
       for i:= 1 to MaxNumVar do m.exp[i]:=0;
       m.exp[numv]:= 1;
       VarPol:= OneTerm(m);
    end
    else VarPol:= Extract (v, ExtVars);
end;


function SumPol (p,q: polynom): polynom;
                    (*  creates the sum of two polynomials
                    referenced by p and q on new space  *)
var
   r: polynom; m: monom;
begin
   if p=nil then r:=CopyPol(q)
   else if q=nil then r:=CopyPol(p)
   else case relation (p^.term, q^.term) of
           cancels: r:= SumPol (p^.next, q^.next);
           similar: begin
                       m:= p^.term;
                       AddRat(m.coef,q^.term.coef,m.coef);
                       new(r);
                       r^.term:= m;
                       r^.next:= SumPol (p^.next, q^.next);
                    end;
           precedes: begin
                       new(r);
                       r^.term:= p^.term;
                       r^.next:= SumPol (p^.next, q);
                     end;
           succeeds:  begin
                       new(r);
                       r^.term:= q^.term;
                       r^.next:= SumPol (p, q^.next);
                     end;
        end;
   SumPol:= r
end;

procedure Add1to2 (p: polynom; var q: polynom);
var r: polynom;
begin
   r:= SumPol (p,q);
   freepol(q);
   q:= r
end;


procedure Negate (p: polynom);
                 (*  creates -p on the same space  *)
var r:polynom;
begin
    r:= p;
    while r<>nil do
    begin
        r^.term.coef.num:= - r^.term.coef.num;
        r:=r^.next
    end;
end;

procedure Add1Times2to3 (m: monom; p: polynom; var q: polynom);

          (*  creates m*p+q on the space occupied by q  *)
          (*  and, if need be, some new space           *)
          (*  while preserving p                        *)
   var
     r: polynom; n: monom;

      procedure MulMon (m,n: monom; var res: monom);
      var i:integer;
      begin
          MulRat (m.coef,n.coef, res.coef);
          for i:=1 to NumVar do
              res.exp[i]:= m.exp[i] + n.exp[i];
          for i:=NumVar+1 to MaxNumVar do res.exp[i]:=0;
      end;

begin            (*  Add1Times2to3  *)
     if p<>nil then
     if q=nil then begin
                    q:= CopyPol(p); r:=q;
                    while r<>nil do
                    begin
                       MulMon (m,r^.term,r^.term);
                       r:= r^.next
                    end;
                 end
     else begin
       MulMon (m,p^.term,n);
       case relation (n, q^.term) of
           cancels:  begin
                        r:= q; q:=q^.next; dispose (r);
                        Add1Times2to3 (m, p^.next, q);
                     end;
           similar:  begin
                        AddRat (n.coef, q^.term.coef,q^.term.coef);
                        Add1Times2to3 (m, p^.next, q^.next);
                     end;
           precedes: begin
                        Add1Times2to3 (m,p^.next,q);
                        new(r);
                        r^.term:= n;
                        r^.next:= q;
                        q:= r;
                     end;
           succeeds: Add1Times2to3 (m, p, q^.next);
       end;    (* case *)
     end;      (* else *)
end;         (* Add1Times2to3 *)

function ProMonPol (m: monom; p: polynom): polynom;
var q: polynom;
begin
    q:= nil;
    Add1Times2to3 (m,p,q);
    ProMonPol:= q;
end;

function ProPol (p1,p2: polynom): polynom;

            (* returns product of p and q on new space  *)

var res: polynom;

begin                               (*  ProPol  *)
    res:= nil;
    if p2<>nil then
    while p1<>nil do
    begin
        Add1Times2to3 (p1^.term, p2, res);
        p1:= p1^.next
    end;
    ProPol:= res
end;                                 (*  ProPol  *)

procedure MulPol (p: polynom; var q: polynom);
             (*  replaces q by the product p*q  *)
             (*  leaving p unchanged            *)
var
    r: polynom;
begin
    r:= ProPol (p,q);
    FreePol (q);
    q:= r
end;

procedure MulCoef (c: rational; p: polynom);
                 (*  creates c*p on the same space  *)
var r:polynom;
begin
    r:= p;
    while r<>nil do
    begin
        MulRat (c,r^.term.coef,r^.term.coef);
        r:=r^.next
    end;
end;

function PowerPol (p: polynom; n: integer): polynom;
             (*  creates p^n on new space  *)
var r,res: polynom; c: rational; m: monom; i: integer;
begin
    if p=nil then if n>0 then res:=nil
                  else if n=0 then Error('exponent 0')
                       else        Error('exponent < 0 ')
    else if p^.next = nil  (* monomial, negative exponents allowed *)
         then if n=0 then res:= ConstPol (1,1)
              else begin
                     m:= p^.term;
                     for i:= 1 to numvar do m.exp[i]:= n*m.exp[i];
                     if n<0 then begin
                                   invrat (m.coef, c);
                                   m.coef:= c;
                                   n:= -n;
                                 end
                     else c:=m.coef;
                     for i:= 1 to n-1 do
                          MulRat (m.coef, c, m.coef);
                     res:= OneTerm (m);
                   end
         else if n<0 then Error( 'exponent < 0 ')
           else if n=0 then res:= ConstPol(1,1)
             else if n=1 then res:= CopyPol(p)
               else if odd(n) then
                        begin
                           r:= PowerPol (p, n div 2);
                           res:= ProPol (r,r);
                           MulPol (p,res);
                           FreePol (r);
                        end
                     else begin
                            r:= PowerPol (p, n div 2);
                            res:= ProPol (r,r);
                            FreePol (r);
                          end;
    PowerPol:= res
end;

procedure fwritenum (var f: text; n: longint);
begin
   write (f,n:0);
end;

procedure writenum (n: longint);
begin
   fwritenum (output, n);
end;
    
function VarLength (v: nametype): integer;
var n: integer;
begin
{    n:= 0;
    while (n+1 <= namelength) and (v[n+1]<>' ') do  n:= n+1;
    VarLength:= n }
    varlength:=length(v)
end;

function MonLength (m: monom) : integer;

    (* returns the number of characters in the standard
                        presentation of m             *)
const
    ln10=2.302585;
var
    i,l: integer;

    function NumLength (n:longint) : integer;
    begin
       if n=0 then Error('zero coefficient') else
       NumLength:= trunc(ln(abs(n))/ln10)+1
    end;

begin
    if m.coef.num=m.coef.den then l:= 3 else
       l:= NumLength(m.coef.num)+3; (* leading sign and 2 spaces incl.*)
    if m.coef.den<>1 then l:= l+ NumLength (m.coef.den)+1;    (* / *)
    for i:=1 to NumVar do if m.exp[i] <> 0 then
    begin
       l:= l + VarLength(VarPres[i]) +1;   (*  '*' and variable *)
       if m.exp[i]<>1 then l:= l + NumLength(m.exp[i])+1;
                                           (* '^' and exponent  *)
       if m.exp[i]<0 then l:= l+1;    (*  '-' sign  *)
    end;
    MonLength:= l;
end;    (* monlength *)

procedure FWriteVar (var outf: text; v: nametype);
var i: integer;
begin
    for i:= 1 to VarLength (v) do
        write(outf,v[i]);
end;

procedure WriteVar (v: nametype);
begin
    FWriteVar (output, v);
end;

procedure FWriteMon (var outf: text; m: monom; first: boolean );
   var i, fnv: integer;  (* first non-vanishing factor *)
   begin
      if not first and (m.coef.num>0) then write(outf,' + ');
      if m.coef.num < 0 then write(outf,' - ');
      m.coef.num:= abs(m.coef.num);
      if (m.coef.num<>m.coef.den) or IsConstant(m) then
      begin
          fwritenum (outf,m.coef.num);
          if m.coef.den<>1 then begin
                                   write(outf,'/');
                                   fwritenum(outf,m.coef.den);
                                end;
      end;
      fnv:= 0; repeat fnv:= fnv+1
               until (fnv>NumVar) or (m.exp[fnv]<>0);
      if fnv <= NumVar then           (* non-constant monomial *)
      for i:= fnv to NumVar do if m.exp[i]<>0 then
      begin
          if (i>fnv) or (m.coef.num<>m.coef.den) then write(outf,'*');
                     (* do not write first *, if m has coefficient +-1 *)
          fwritevar (outf,VarPres[i]);
          if m.exp[i]<>1 then begin
                                 Write(outf,'^');
                                 fwritenum(outf,m.exp[i])
                              end;
      end;
   end;

procedure FWritePol (var outf: text; p: polynom; var offset: integer );
                        (* writes to an abstract text file always
                        keeping track of the change of margin offset *)
var i,monlen: integer;
begin
    if p=nil then write (outf,' 0') else
    begin i:=1;
       repeat
           monlen:= MonLength (p^.term);
           if monlen > (78 - offset) then
           begin
              writeln (outf);
              offset:=1;
           end;
           FWriteMon (outf,p^.term, i=1);
           offset:= offset + monlen;
           p:=p^.next;
           i:= i+1
       until p=nil;
    end;
end;

procedure FNameWritePol (filename: string; p: polynom);
var f: text; unity: integer;
begin
    assign (f, filename);
    rewrite (f);
    unity:= 1;
    FWritePol (f, p, unity);
    close (f);
end;

procedure writemon (m: monom);
begin
   fwritemon (output,m,true);
end;

procedure WritePol (p: polynom );
                                (* writes to console *)
var dummy: integer;
begin
    dummy:= 1;
    FWritePol (output, p, dummy);
end;

procedure KillBlanks (var s: string);
var  i: integer; snew: string;
begin
   snew:='';
   for i:= 1 to length(s) do if s[i]<>' ' then
        snew:= snew + s[i];
   s:= snew;
end;

function char2int (c:char): integer;
begin
    char2int:=ord(c)-ord('0')
end;

procedure ReadLInt (s: string; var i: integer; var a: longint);

     (* reads part of a line representing an integer
      and shifts i to the first non-digit  *)

  begin
     if length(s)=0  then a:=1 else
        if s[i]='+' then
        begin
           i:= succ(i); ReadLInt(s,i,a)
        end
        else if s[i]='-' then
        begin
           i:= succ(i); ReadLInt(s,i,a); a:=-a;
        end
        else if s[i] in digits then
        begin
           a:=char2int(s[i]); i:= succ(i);
               while (s[i] in digits) and (i<=length(s)) do
               begin
                  a:=a*10 + char2int(s[i]) ;
                  i:= succ(i);
               end;
        end
        else a:=1
  end;         (*   ReadLInt   *)

procedure ReadIdent (s: string; var i: integer; var id: nametype);
         (* read an identifier out of a string s beginning from position i
         and shift the index i correspondingly *)
begin
   id:='';
   while (i<=length(s)) and (s[i] in digits+letters) do
   begin
      id:= id + s[i];
      i:= succ(i);
   end;
end;

procedure getlex (s: string; var i: integer; predop: opertype; var l: lexem);
var ch: char; newname: nametype;
    j,len: integer;
begin
   len:= length(s);
   ch:= s[i];
   if ch in delimiters then
   begin
      l.which:= oper;
      case ch of
        '+': if predop in [nonop,rightpar] then l.whato:= biplus
             else if predop = arrow then l.whato:= postarplus
                 else l.whato:= unplus;
        '-': if predop in [nonop,rightpar] then l.whato:= biminus
             else if predop = arrow then l.whato:= postarmin
                  else l.whato:= unminus;
        '*': l.whato:= aster;
        '/': l.whato:= slash;
        '^': l.whato:= arrow;
        '(': l.whato:= leftpar;
        ')': l.whato:= rightpar;
        ',': l.whato:= commut;
      end;
      i:= i+1;
   end
   else if ch in letters then
   begin
      l.which:= ident;
      ReadIdent (s,i,newname);    (* as a side-effect, this
                               shifts i to the position after ident*)
      if length (newname) > namelength then Error('namelength exeeded')
      else begin
              l.whats:= newname;
{              store(newname); (* in the array VarPres *)  }
           end;
   end
   else if ch in digits then
   begin
      l.which:= number;
      ReadLInt (s,i,l.whatn);
   end
   else Error ('Illegal character '+ch);
end;

procedure addtolexlist (var f: lexlist; l: lexem);
var
   g: lexlist;
begin
   new (g);
   g^.item:= l;
   g^.next:= f;
   f:= g;
end;

function Invert (p: lexlist): lexlist; (* I.Nesterov *)
var pb,pf: lexlist;
begin
   pb:= nil;
   while p<>nil do
   begin
      pf:= p^.next;
      p^.next:= pb;
      pb:= p;
      p:= pf;
   end;
   invert:= pb;
end;

procedure s2f (s: string; var f: lexlist);
var i: integer; l: lexem; predop: opertype;
begin
   i:= 1;
   KillBlanks (s);
   f:= nil;
   predop:= beginop;
   while i<= length(s) do
   begin
      getlex (s,i,predop,l);
      addtolexlist (f,l);
      if l.which = oper then predop:= l.whato
          else predop:= nonop;
   end;
   f:= Invert (f);
end;

function charimg (op: opertype): string;
   (* auxiliary function; for debugging only *)
begin
   case op of
       unplus: charimg:= 'u+';
       unminus : charimg:= 'u-';
       postarplus: charimg:= 'p+';
       postarmin  : charimg:= 'p-';
       biplus   : charimg:=  'b+';
       biminus : charimg:= 'b-';
       aster   : charimg:=  '*';
       slash  : charimg:=  '/';
       arrow  : charimg:=  '^';
       leftpar  : charimg:= '(';
       rightpar : charimg:= ')';
       commut  : charimg:= ',';
       nonop : charimg:=   'non';
       beginop : charimg:= 'beg';
    end;
end;

procedure displlexlist (l: lexlist);
begin
   while l<>nil do
   begin
     case l^.item.which of
       oper: writeln ('oper: ',charimg(l^.item.whato));
       ident: writeln ('iden: ',l^.item.whats);
       number: begin
                  write ('numb: ');
                  writenum (l^.item.whatn);
                  writeln;
               end;
     end;
     l:= l^.next;
   end;  
end;  

procedure displtree (t:expr);
begin
   case t^.tag of
       binary: begin
                 writeln('binary ',charimg(t^.biop));
                 displtree(t^.arg1);
                 displtree(t^.arg2);
               end;
       unary:  begin
                 writeln('unary ',charimg(t^.unop));
                 displtree(t^.arg);
               end;
       number: write(' ',t^.num,' ');
       ident:  write(' ',t^.nam,' ');
   end;
end;

procedure transpostfix (inlist: lexlist; var outlist: lexlist);
type
   stack = ^stacknode;
   stacknode = record
                  op:  opertype;
                  next: stack;
               end;
var
  shop: stack;
  lex: lexem;
  op: opertype;

    procedure displstack (s: stack);
    begin
       while s<>nil do
       begin
          write (charimg(s^.op),' ');
          s:= s^.next;
       end;
    end;

    function prior (d: opertype): integer;
    begin
       case d of
          leftpar:         prior:= 0;
          biplus,biminus,
               rightpar:   prior:= 1;
          commut:          prior:=2;
          aster, slash:          prior:= 3;
          unplus, unminus:       prior:= 4;
          arrow:                 prior:= 5;
          postarplus,postarmin:  prior:= 6;
          else  Error('priority undefined');
        end;
    end;

    procedure push (nextop: opertype; var st: stack);
            (* put new item on top of the stack *)
    var s: stack;
    begin
       new(s);
       s^.op:= nextop;
       s^.next:= st;
       st:= s;
    end;

           procedure RemoveTop (var st: stack);
           var s: stack;
           begin
              s:= st^.next;
              dispose (st);
              st:= s;
           end;

    procedure pop (var st: stack; newop: opertype; var f: lexlist);
            (* move all items on top of the stack
              whose priority is >= prior(newop) to output;
              leftpar pops nothing;  rightpar pops everything
              to nearest leftpar and annihilates with the latter *)
    label 21,22;
    var l: lexem; pr: integer;

    begin
{writeln('entering pop');}
       if newop= leftpar then   (* do nothing *)
       else if newop= rightpar then
       begin
         while st<>nil  do 
         begin
           if st^.op=leftpar then goto 21; 
           l.which:= oper;
           l.whato:= st^.op;
           if l.whato = rightpar then Error('excessive right parenthesis')
              else addtolexlist (f,l);
           RemoveTop (st);
         end;
21:      if (st=nil) then Error('left parenthesis missing')
             else if st^.op<>leftpar then Error('Left parenthesis missing')
                 else RemoveTop (st);
       end else
begin{}
{if st=nil then writeln('st=nil') else 
    begin
       writeln('st<>nil');
       writeln('st^.op=', charimg(st^.op));
    end;}   
       while (st<>nil)  do
       begin
          if prior (st^.op) < prior(newop)  then goto 22;
          l.which:= oper;
          l.whato:= st^.op;
          addtolexlist (f,l);
{writeln('going to remove top of stack');
write('current state of the stack: '); displstack (st);
writeln;}
          RemoveTop (st);
       end;
22:
end;{}
    end;

    procedure clearstack (st: stack; var outl: lexlist);
    var l: lexem;
    begin
         while st<>nil do
         begin
           l.which:= oper;
           l.whato:= st^.op;
           if l.whato <> leftpar then addtolexlist (outl,l);
           RemoveTop (st);
         end;
    end;

begin
   outlist:= nil;
   shop:= nil;
   while inlist<>nil do
   begin
      lex:= inlist^.item;
      case lex.which of
         oper: begin
                  op:= lex.whato;
                  pop (shop, op, outlist);
                  if op <> rightpar then push (op,shop);
               end;
         ident,number:  addtolexlist (outlist, lex);
      end; (* case *)
      inlist:= inlist^.next
   end;
   clearstack (shop, outlist);
end;     (* transpostfix *)    


procedure lex2expr (lex: lexem; var en: enode);
         (* translate a single lexem *)
begin
    with en do
    begin
      case lex.which of
         oper: begin
                 case lex.whato of
                   biplus..arrow: begin
                                     biop:= lex.whato;
                                     tag:= binary;
                                     arg1:= nil;
                                     arg2:= nil;
                                  end;
                   unplus..postarmin: begin
                               unop:= lex.whato;
                               tag:= unary;
                               arg:= nil;
                            end;
                   else Error('error in operation');
                 end;
               end;
         ident:  begin
                    tag:=ident;
                    nam:= lex.whats;
                 end;
         number: begin
                    tag:=number;
                    num:= lex.whatn;
                 end;
       end;
    end;
end;

function list2expr (l: lexlist): expr;  (* I.Nesterov *)
label 25;
const depth = 30;  (* stack depth  *)
var r, newr, head: expr;
    stack: array[1..depth] of expr;
    sp: integer;      (* stack pointer *)
begin
   if l=nil then head:= nil else new (head);
   r:= head;
   sp:=0;
   while l<>nil do
   begin
     lex2expr (l^.item,r^);
25:  case r^.tag of
        number,ident: if sp<>0 then begin
                         r:= stack[sp];
                         sp:= sp-1; goto 25;
                       end;
        unary: if r^.arg = nil then
                  begin
                    sp:= sp+1;
                    if sp>depth then Error(' stack too shallow');
                    stack[sp]:= r;
                    new (newr);
                    r^.arg:= newr;
                    r:= newr;
                  end
                  else if sp<>0 then begin
                         r:= stack[sp];
                         sp:= sp-1; goto 25;
                       end;
        binary: if r^.arg1 = nil then
                  begin
                    sp:= sp+1;
                    if sp>depth then Error(' stack too shallow');
                    stack[sp]:= r;
                    new (newr);
                    r^.arg1:= newr;
                    r:= newr;
                  end
                  else if r^.arg2 = nil then
                  begin
                    sp:= sp+1;
                    if sp>depth then Error(' stack too shallow');
                    stack[sp]:= r;
                    new (newr);
                    r^.arg2:= newr;
                    r:= newr;
                  end
                  else if sp<>0 then begin
                         r:= stack[sp];
                         sp:= sp-1; goto 25;
                       end;
     end;  (* case *)
     l:= l^.next;
   end;   (* while *)
   list2expr:= head ;
end;

procedure cmultmon (m,n:monom; var mn: monom);
  (* e.g., 2*AB, 3*BA -->  6*ABBA *)
var
   i,j:integer;
   c1,c2,c: rational;
   id1,id2,id: nametype;
begin
   i:=0; repeat inc(i) until m.exp[i]<>0;
   id1:= varpres[i]; c1:=m.coef;
   j:=0; repeat inc(j) until n.exp[j]<>0;
   id2:= varpres[j]; c2:=n.coef;
   id:=id1+id2;
   mulrat(c1,c2,mn.coef);
   for i:=1 to numvar do mn.exp[i]:=0;
   mn.exp[nvar(id)]:=1;
end;

function CCommut (p,q:polynom): polynom;
   (* concatenation commutator of polynomials;
      result on new space *)
var
   qr,r,s: polynom;
   m,n,k: monom;
begin
   r:=nil;
   while p<>nil do
   begin
      m:=p^.term;
      qr:=q;
      while qr<>nil do
      begin
         n:=qr^.term;
         cmultmon(n,m,k);
         s:=OneTerm(k);
         Add1To2 (s,r);
         freepol (s);
         cmultmon(m,n,k);
         s:=OneTerm(k);
         negate(s);
         Add1To2 (s,r);
         freepol(s);
         qr:=qr^.next
      end;
      p:=p^.next
   end;
   CCommut:=r;
end;

function Cproduct (p,q:polynom): polynom;
   (* concatenation product of polynomials;
      result on new space *)
var
   qr,r,s: polynom;
   m,n,k: monom;
begin
   r:=nil;
   while p<>nil do
   begin
      m:=p^.term;
      qr:=q;
      while qr<>nil do
      begin
         n:=qr^.term;
         cmultmon(m,n,k);
         s:=OneTerm(k);
         Add1To2 (s,r);
         freepol(s);
         qr:=qr^.next
      end;
      p:=p^.next
   end;
   Cproduct:=r;
end;

procedure rev_str (s: nametype; var r: nametype);  (* reverse a string *)
var 
   len,i:integer;
begin
   r:='';
   len:=length(s);
   for i:=1 to len do r:=r+s[len+1-i];
end;

procedure rev_mon (var m:monom);
        (* reverse all variables in a monomial, e.g. AABA -> ABAA *)
var 
   n:monom;
   i,e:integer;
   v: nametype;
begin
   n.coef:=m.coef;
   for i:=1 to maxnumvar do n.exp[i]:=0;
   for i:=1 to maxnumvar do
   begin
      e:=m.exp[i];
      if e<>0 then 
      begin
         rev_str(varpres[i],v);
         n.exp[nvar(v)]:=e
      end;
   end;
   m:=n;
end;      

procedure rev_pol (var p: polynom); 
        (* reverse all variables in a polynomial, e.g. AABA -> ABAA *)
var 
   q,r:polynom;
   m: monom;
begin
   q:=nil;
   while p<>nil do
   begin
      m:= p^.term;
      rev_mon (m);
      r:=oneterm (m);
      add1to2 (r,q);
      freepol(r);
      p:=p^.next
   end;
   freepol(p);
   p:=q;
end;

function tree2pol (t:expr): polynom;
var p,q,r: polynom; n: longint;
begin
   if t=nil then p:= nil
   else case t^.tag of
      unary:  case t^.unop of
                unminus:  begin
                             p:= tree2pol(t^.arg);
                             negate (p);
                          end;
                unplus:   p:= tree2pol(t^.arg);
                postarmin: begin
                             p:= tree2pol(t^.arg);
                             negate (p);
                           end;
                postarplus:  p:= tree2pol(t^.arg);
              end;
      binary: case t^.biop of
                biplus: begin
                           q:= tree2pol(t^.arg1);
                           r:= tree2pol(t^.arg2);
                           p:= sumpol(q,r);
                           freepol(q); freepol(r);
                        end;
                biminus: begin
                            q:= tree2pol (t^.arg2);
                            r:= tree2pol (t^.arg1);
                            negate (r);
                            p:= sumpol (q,r);
                            freepol(q); freepol(r);
                         end;
                arrow: begin
                          r:= tree2pol (t^.arg1);
                          if r=nil then n:=0 else
                          if (r^.next=nil) and (isconstant (r^.term))
                                and (r^.term.coef.den=1)
                                then n:= r^.term.coef.num
                             else Error ('non-numeric exponent');
                          freepol (r);
                          q:= tree2pol (t^.arg2);
                          p:= PowerPol (q,n);
                          freepol(q);
                       end;
                slash: begin
                          q:= tree2pol (t^.arg1);
                          if q=nil then Error('zero denominator') else
                          if q^.next=nil  then
                          begin
                             r:= PowerPol (q,-1); freepol(q);
                          end  else Error ('non-monomial divisor');
                          q:= tree2pol (t^.arg2);
                          p:= ProPol (q,r);
                          freepol(q); freepol(r);
                       end;
                aster: begin
                           q:= tree2pol(t^.arg1);
                           r:= tree2pol(t^.arg2);
                           p:= propol(q,r);
                           freepol(q); freepol(r);
                        end;
                commut: begin
                           q:= tree2pol(t^.arg1);
                           r:= tree2pol(t^.arg2);
                           p:= ccommut(q,r);
                           freepol(q); freepol(r);
                        end;
              end;
      ident:  p:= VarPol (t^.nam);
      number: p:= ConstPol (t^.num,1);
  end;
  tree2pol:= p;
end;

procedure freelist (l: lexlist);
var r: lexlist;
begin
   while l<>nil do
   begin
      r:= l^.next;
      dispose (l);
      l:= r;
   end;
end;

procedure freeexpr (ex: expr);
begin
   if ex<>nil then
   case ex^.tag of
      ident, number: dispose (ex);
      unary: begin
               freeexpr (ex^.arg);
               dispose (ex);
             end;
      binary: begin
               freeexpr (ex^.arg1);
               freeexpr (ex^.arg2);
               dispose (ex);
             end;
     end;
end;

function s2p (s: string): polynom;
              (* transforms string into polynomial *)
var
   infix, postfix: lexlist;
   ex: expr;
   p:polynom;
begin
{writeln('s2p entered');}
   s2f (s,infix);
{writeln('transformed to infix');
displlexlist (infix);  writeln;}
   transpostfix (infix, postfix);   freelist (infix);
{writeln('transformed to postfix');
displlexlist (postfix);  writeln;}
   ex:= list2expr (postfix);        freelist (postfix);
{writeln('transformed to tree');
displtree(ex); readln;}
   p:= tree2pol (ex);             freeexpr (ex);
{writeln('transformed to polynom');
writepol(p); readln;}
   s2p:=p;
end;

function copy (s: string; st,le: integer) : string;
var
   i: integer;
   c: string;
begin
   c:='';
   for i:= st to st+le-1 do  
          c:= c+ s[i];
   copy:= c;
end;


procedure FReadPol (var inpf: text; var p: polynom);
         (* reading from an abstract text file, provided that every line
            has correct syntax and the whole polynomial
            is the sum of these lines; input is terminated either by eof
            or by semicolon *)
var
   buf: string;
   temp: polynom;
   nosemicolon: boolean;
begin
   p:= nil;
   nosemicolon:= true;
   while nosemicolon and not eof(inpf) do
   begin
      readln(inpf,buf);
      if buf='' then temp:=nil
      else if buf[1]='%' then temp:= nil
                          (*  %  signifies a line of comment *)
      else if buf[length(buf)] <> ';' then
         begin
            nosemicolon:= true;
            temp:= s2p(buf);
         end
         else   (* last semicolon *)
         begin
            nosemicolon:= false;
            temp:= s2p(copy(buf,1,length(buf)-1));
         end;
      Add1to2(temp,p);
      freepol(temp);
   end;
end;

procedure ReadPol (var p: polynom);
            (* reading from keyboard, terminated with ; *)
begin
    FReadPol (input, p);
end;

procedure FNameReadPol (fname: string; var p: polynom);
var inpf: text;
begin
    assign (inpf, fname);
    reset (inpf);
    FReadPol (inpf,p);
    close (inpf);
end;

procedure FReadOrdPol (var inpf: text; var p: polynom);
            (* reading from an abstract text file, provided that its syntax
               is such as if it were output by FWritePol:
               monomials are arranged by descend of order;
               termination by semicolon  *)
var
   buf: string;
   temp: polynom;
   nosemicolon: boolean;

   procedure  Append1to2 (t: polynom; var q: polynom);
   var r: polynom;
   begin
      if q=nil then q:= t else
      begin
         r:= q;
         while r^.next<>nil do r:= r^.next;
         r^.next:= t;
      end;
   end;

begin
   p:= nil;
   nosemicolon:= true;
   while nosemicolon and not eof(inpf) do
   begin
      readln(inpf,buf);
      if buf='' then temp:=nil
      else if buf[1]='%' then temp:= nil
                          (*  %  signifies a line of comment *)
      else if buf[length(buf)] <> ';' then
         begin
            nosemicolon:= true;
            temp:= s2p(buf);  (* develop to SOrd2p ! *)
         end
         else   (* last semicolon *)
         begin
            nosemicolon:= false;
            temp:= s2p(copy(buf,1,length(buf)-1)); (* develop to SOrd2p ! *)
         end;
      Append1to2(temp,p);
   end;
end;

procedure FNameReadOrdPol (fname: string; var p: polynom);
var inpf: text;
begin
    assign (inpf, fname);
    reset (inpf);
    FReadOrdPol (inpf,p);
    close (inpf);
end;

procedure defsets;
begin
    delimiters := ['(', ')', ',', '+', '-', '*', '/', '^'];
    letters := ['a'..'z','A'..'Z'];
    digits := ['0'..'9']; 
end;

procedure initstr (var s:string; n: integer);
var i: integer;
begin
   s:='';
   for i:=1 to n do s:=s+'A';
end;

procedure nextstr (var s:string; var last:boolean);
var 
   i,r,l:integer;
begin
   last:=false;
   l:=length(s);
   r:=l;
   while (r>0) and (s[r]='B') do dec(r);
   if r=0 then last:=true else
   begin
      s[r]:='B';
      for i:=r+1 to l do s[i]:='A'
   end;
end;

procedure InitVarsAB(n:integer);
var
   s:string;
   last:boolean;
   k:integer;
   vpr:varsettype;
begin
   NumVar:=0;
   for k:=1 to n do
   begin
      initstr(s,k);
      repeat
         inc(NumVar);
         varpres[NumVar]:=s;
         nextstr(s,last)
      until last
   end;
{   for k:=1 to numvar do  (* reverse the sequence so that B becomes   *)
                          (* higher than A  *)
   begin
      vpr[k]:=varpres[numvar+1-k];
   end;
   varpres:=vpr;}
end;

procedure initvec (var v:vector);
var i: integer;
begin
   for i:=1 to numvar do v[i]:=0;
end;

procedure pol2vec (p:polynom; var v:vector);
var k:integer;
begin
   initvec (v);
   while p<>nil do
   begin
      k:=0;
      repeat inc(k) until p^.term.exp[k]=1;
      v[k]:=v[k]+p^.term.coef.num;
      p:=p^.next
   end;
end;

procedure fwritevec (var t:text; v:vector);
begin
   write(t,'[');
   for i:= (numvar div 2)+1 to numvar-2 do write(t,v[i],',');
   write(t,v[numvar-1],']');
end;

procedure swapAB (s: string; var s1:string);
var
   i: integer;
   c: char;
begin
   s1:='';
   for i:=1 to length(s) do
   begin
      c:=s[i];
      if c='A' then c:='B' else if c='B' then c:='A';
      s1:=s1+c;
   end;
end;

begin
   defsets;
{ edit the curly braces below to get the required version }
   write('n='); readln(n);
   InitVarsAB(n);
   if n<10 then ext:=chr(ord('0')+n)
      else if n=10 then ext:='10'
      else if n=11 then ext:='11'
      else if n=12 then ext:='12';
   assign(infile,'hall_bas.'+ext+'f');
   reset(infile);
   assign(outfile,'bas_comm.'+ext+'f');
   rewrite(outfile);
   assign(outft,'bas_comm.'+ext+'tf');
   rewrite(outft);
   assign(outfr,'bas_comm.'+ext+'rf');
   rewrite(outfr);
   writeln(outfile,'matr'+ext+' := [');
   count:=0;
   repeat
      inc(count);
      readln(infile,s);
      p:=s2p(s);
      writeln(count);
{      write('p[',n,count,'] := ');
      offset:=3+n;
      writepol(p); writeln(';');}
      pol2vec (p,v);
      fwritevec (outfile,v);
      writeln(outfile,',');
      rev_pol (p);
      pol2vec (p,v);
      write (outfr,'vr[',count,']:=');
      fwritevec (outfr,v);
      writeln(outfr,';');
      freepol(p);
      swapAB (s,s1);
      p1:=s2p(s1);
      pol2vec (p1,v);
      write (outft,'vt[',count,']:=');
      fwritevec (outft,v);
      writeln(outft,';');
      freepol(p1);
   until eof(infile);
   writeln(outfile,'];');
   writeln(outft,'size:=',count,';');
   writeln(outfr,'size:=',count,';');
   close(outfile);
   close(outft);
   close(outfr);
   close(infile);
end.
