(* updated the output procedure to work better with Maple, July 1998 *)
(* Original program adapted in 1994 by H.R.Morton from the 1985 Hecke algebra
program of Morton and Short *)

program br9vas;   {Does not use pointers for polynomial handling}
const
   maxl=120; {maximum braid length}
   Vassmax=20;{Vassiliev invariants of at most this degree will be found}                                                                
   Nfact= 362880;      (* 362880 *)                                             
   dNfact= 40320;      (* 40320 *)                                              
   N=9;
type
     polynomial= record coefft:array[0..Vassmax] of integer;
                 maxdeg:integer;
                 mindeg:integer
          end;

     ptr1=polynomial;

var a,aa,b,bb,bs,blocklength,s,g,i,j,len,len1,r,sig:integer;                         
    test,t,tr:integer;        
    vassiliev_type:integer:=Vassmax;
    x,y,poly,zeropoly:polynomial;
    t1,t2:polynomial;                                                  
    n,K,width, nfact,dnfact : integer;
    algcrt  :integer;                                                           
    mistake: boolean;                                                           

    inv: array[1..Nfact] of 1..Nfact;                                         
    turn: array[1..Nfact,1..N] of 1..Nfact;
    calculated: array[1..N] of boolean;	
    braid:array[0..maxl] of -N..N;                                              
    braidstr, braidstr1 : varying[200] of char;
    factor:array[0..N] of 1..Nfact;                                             
    ch, ch1:char;                                                                    
    c:array[1..Nfact] of polynomial;                                               
    d:array[1..dNfact] of polynomial;                                              
    mult:array[1..Nfact,1..N-1] of integer;                                      
    u:array[0..N] of polynomial;                                                   
    v:array[0..2*N] of polynomial;                                               
    infile:text;

               
 function plus(a,b:polynomial):polynomial;                                       
  var c:polynomial; i:integer;
   begin  if a.mindeg<0 then plus:=b else
   if b.mindeg<0 then plus:=a else
   begin
      c := zeropoly;
      c.mindeg:=min(a.mindeg,b.mindeg);
         c.maxdeg:=max(a.maxdeg,b.maxdeg);
         for i:=c.mindeg to c.maxdeg do
         c.coefft[i]:=a.coefft[i]+b.coefft[i];
         plus:=c;
      end;
   end;{plus}

function minus(a:polynomial):polynomial;                            
 var c:polynomial; i:integer;
   begin  if a.mindeg<0 then minus:=zeropoly else

      begin
      c:=zeropoly;   
      c.mindeg:=a.mindeg;
         c.maxdeg:=a.maxdeg;
         for i:=c.mindeg to c.maxdeg do
         c.coefft[i]:=-a.coefft[i];
         minus:=c;
      end;{else}
   end;{minus}

function shift(a:polynomial):polynomial;{multiplies a polynomial in z by the variable z}
var i: integer;
    x:polynomial;
begin
   if a.mindeg<0 then shift:=zeropoly else
   begin
      x:=zeropoly;
      x.mindeg:=a.mindeg+1;
      x.maxdeg:=min(a.maxdeg+1,vassiliev_type);{allows truncation of the invariant at a given degree in z, specified during calculation}

      for i:=x.mindeg to x.maxdeg do x.coefft[i]:=a.coefft[i-1];
      shift:=x;
   end;
end;{shift}

procedure multtab;                                                              
var  g,h,r,cr,cr1,j,k,m,f,ff: integer;                                        
begin                                                                           
   k:=1;                                                                        
   for g:= 1 to nfact do                                                        
   begin                                                                        
      mult[g,1]:=g+k;                                                           
      k:=-k;                                                                    
   end;                                                                         
   for r:=2 to n-1 do                                                            
   begin                                                                        
      f:=factor[r]; ff :=factor[r-1];                                           
      for g:= 1 to nfact do                                                     
      begin                                                                     
         m:= (g-1) div f;                                                       
         cr:= m mod (r+1);                                                      
         m:= (g-1) div ff;                                                      
         cr1:= m mod r;                                                         
         j:= cr1 - cr;                                                          
         if j>= 0 then                                                          
         begin                                                                  
            h:= g + j*ff*(r-1) + f;                                             
            mult[g,r]:= h;                                                      
            mult[h,r]:= g;                                                      
         end;                                                                   
      end;                                                                      
   end;                                                                         
end; {multtab}                                                                  
                                                                                
                                                                                
procedure reverse(r:integer);                                                   
var g,s : integer;                                                          
    p:polynomial;
begin                                                                           
   for g:= 1 to factor[r] do                                                    
   begin                                                                        
      if g < inv[g] then                                                        
      begin                                                                     
         p:= c[g]; s:= inv[g];                                                  
         c[g]:= c[s]; c[s]:= p;                                                 
      end;                                                                      
   end;                                                                         
end;{reverse}                                                                   
                                                                                
procedure listinverse;                                                          
var p,m,cr,g,s,r,f: integer;                                                  
begin                                                                           
   inv[1]:=1;                                                                   
   for s:= 1 to n-1 do                                                         
   begin                                                                        
      f:=factor[s];                                                             
      for g:=f+1 to factor[s+1] do                                              
      begin                                                                     
         p:= inv[g-f];                                                          
         m:= (g-1) div f;                                                       
         cr:= m mod (s+1);                                                      
         r:=s+1-cr;                                                             
         inv[g]:= mult[p,r];                                                    
      end;                                                                      
   end;                                                                         
end; {listinverse}                     

procedure turnover(n:integer);
var g,i: integer;
begin
   if calculated[n]=false then
   begin                                         
	writeln;
	writeln('The table for turned-over braids on',n:3,' strings is being prepared.');
	turn[1,n]:=1;
	for g:=1 to factor[n] do	
	begin
	   for i:=1 to n do
	   begin
		if mult[g,i]>g then
		turn[mult[g,i],n]:=mult[turn[g,n],n-i];
	   end; {i}
	end; {g}
	calculated[n]:=true;
   end;
end; {turnover}
		                                                                                
procedure multiply(a,b:integer);                                                
var sign,g,h: integer;                                                          
    x,y:polynomial;                              
begin                                                                           
   if r>0 then sign:= 1 else sign:= -1;                                         
   r:=abs(r);                                                                   
   for g:=a to b do                                                             
   begin                                                                        
      h:= mult[g,r];                                                            
      if h > g then                                                             
      begin                                                                     
         if (c[g].mindeg>-1) or (c[h].mindeg>-1) then
         begin                                                        
            x:=c[g]; y:=c[h]; 
            if sign>0 then                                                         
            begin                                                                  
               c[h]:=plus(x,shift(y));
               c[g]:=y;
            end else                                                               
            begin                                                                  
               c[g]:=plus(y,minus(shift(x)));
               c[h]:=x;
            end;  
         end;                                                                   
      end;                                                                      
   end;                                                                         
 end; {multiply}                                                                
                                                                                
procedure wrapup;                                                               
var i,j,k,s,a,b,g,rr,f,ff:integer;                                        
begin                                                                           
   for s:=n-1 downto 1  do                                                      
   begin                                                                        
      b:=0; f:=factor[s];                                                       
      for i:=0 to n-s-1 do                                                      
      begin                                                                     
         b:=b+factor[s+i]; a:=b+1-f;                                            
         for rr:=1 to s-1 do                                                    
         begin                                                                  
            k:=(s-rr+1)*f;r:=rr;                                                
            multiply(a+k,b+k);                                                  
            for g:=a+k to b+k do                                                
            begin
               if c[g].mindeg>-1 then
               begin
                  c[g-f]:=plus(c[g],c[g-f]);
               end;
            end;{g}
         end;{rr}                                                                   
         ff:=f-factor[s+i];                                                     
         for g:=a to b do if (i>0)  then                                        
            begin
               if c[g].mindeg>-1 then
               begin
                  c[g+ff]:=plus(c[g],c[g+ff]);
               end;
            end;
      end;{i}                                                                      
   end;{s}                                                                         
   f:=0;                                                                        
   for j:=0 to n-1 do                                                           
   begin                                                                        
      f:=f+factor[j];                                                           
      u[j]:=c[f];                                                               
   end;                                                                         
end; {wrapup}                                                                   
                                                                                
procedure init;                                                                 
var i,j : integer;                                                              
begin                                                                           
    writeln;
    writeln('This program will handle braids of up to ',N:1,' strings');
    writeln('Give the number of strings in the largest braid to be calculated');

    readln(n);
    if n>N then writeln('String index too large for the program');
    
    for j:=1 to n do calculated[j]:= false; (* waiting for a table to be calculated *)
  
    zeropoly.mindeg:=-1000;zeropoly.maxdeg:=-1000;
    for i:=0  to Vassmax do zeropoly.coefft[i]:=0;

    for j:=0 to n do  u[j]:=zeropoly;
    for j:=0 to 2*n do  v[j]:=zeropoly;                                              
    factor[0]:=1;                                                               
    for j:= 1 to n do factor[j]:= j*factor[j-1] ;                              
  
    for i:=2 to factor[n] do  c[i]:=zeropoly;                                              
    c[1].mindeg:=0;
    c[1].coefft[0]:=1;
    c[1].maxdeg:=0;
    for j:=1 to Vassmax do c[1].coefft[j]:=0;
    nfact:=factor[n];
    dnfact:=factor[n-1];



    writeln('Please wait while multiplication table is calculated');
    multtab;                                                                    
    writeln('Now calculating the table of inverses');
    writeln;
    writeln;
    listinverse;                                                                

end; {init}                                                                     
                                                                                
procedure reinit;                                                               
var i,j : integer;                                                              
begin                                                                           
    for j:=0 to n do  u[j]:=zeropoly;                                                  
    for j:=0 to 2*n do  v[j]:=zeropoly;                                              
    for i:=2 to nfact do  c[i]:=zeropoly;
    c[1].mindeg:=0;
    c[1].coefft[0]:=1;
    c[1].maxdeg:=0;
    for j:=1 to Vassmax do c[1].coefft[j]:=0;

end; {reinit}                                                                   
                                                                                
procedure symmetrise;                                                           
var t,g,s: integer;      (* b = factor[n-K], t gives degree of block *)           
begin               
   if ch1='t' then turnover(n-K);                                                            
   for t:= 0 to K do                                                          
   begin                                                                        
      for g:= 1 to b do                                                         
      begin                      
         if ch1='r' then  s:= inv[g]
         else if ch1='t' then s:=turn[g,n-K];      
         if s >= g then                                                         
         begin                                                                  
            c[g+t*b]:=plus(c[g+t*b],minus(c[s+t*b]));
            c[s+t*b]:=minus(c[g+t*b]);
         end;                                                                   
      end; {g}                                                                     
   end; {t}                                                                        
end;{symmetrise}                                                                
                                                                                
                                                                                
procedure fileout;                                                              
var i,k:integer;
    test, first:boolean;
    poly:polynomial;                                                                  
begin                                                                           
      writeln('Link constructed from the braid ',braidstr1);
      writeln('with ',K:1,' string(s) closed off');
      writeln('and composed with the braid ',braidstr);  
      writeln;                                                         
      writeln('algebraic crossing number =',algcrt:3);                           
      writeln;                                                         
      writeln('The Homfly polynomial P(v,z) is the following polynomial in u=(1-v^2)/z and z,');
      writeln('multiplied by v^',algcrt-n-K+1:3);
      writeln;                                                      
      begin test:=false;
         for i:=0 to vassiliev_type do
         begin  first:=true;
         for k:=0 to width-1  do
            begin
                poly:=v[width-1-k];
               if poly.coefft[i]<>0 then
               begin
                   if first=true then
                        begin
                        if test = false then
                            begin
                            write('  (',poly.coefft[i]:2,' *u^',k:2);                    
                            end
                        else write(' +(',poly.coefft[i]:2,' *u^',k:2);
                        first:=false;
                        end
                      else
                        begin
                        if poly.coefft[i]>0 then
                            begin
                            write(' +',poly.coefft[i]:2,' *u^',k:2);                     
                            end
                        else write('  ',poly.coefft[i]:2,' *u^',k:2);
                        end;
                   test:=true;
               end;
            end; {k}
         if first=false then
         writeln(') *z^',i:2);
         end; {i}
         if test= false then writeln('0');
      end;
      writeln(' + O(z^',vassiliev_type+1:2,')');
      writeln;
end; {fileout}
                                                                                
begin; {main program}                
   init;
   repeat
   writeln('This program allows you to specify a ',n:1,'-string braid,');
   writeln('and then close off some of the strings to get a tangle on fewer strings,');
   writeln('before multiplying by a further ',n:1,'-string braid, and closing.');
   writeln('The output is the Homfly polynomial,');

   writeln;
   writeln;
   writeln('ENTER FIRST BRAID, using 1,2,.. for the standard braid generators ');
   writeln('and -1,-2,... for their inverses.');                                        
   readln(braidstr);                                                     

   writeln('The invariant will be truncated at terms of degree', vassiliev_type:3);
   writeln('If you want to alter this degree then type y.');
   writeln('Otherwise press RETURN twice to continue.');
   read(ch);
   if (ch ='y') then 
   begin
      writeln('The present  degree for truncation is ',vassiliev_type:1,' Enter the new value.');
      readln(vassiliev_type);
   end;

   writeln('The braid is to have some strings closed off at the right-hand side,');
   writeln('before the next part is added.');
   writeln;
   writeln('Enter the number of strings, anything from 0 to ',n-1:1);
   writeln('to be closed off.');
   readln(K);
   if (K>n-1) then writeln('You are trying to close too many strings.');
   
   width:=n+K;
   len:=length(braidstr);                                                 
   sig:=1; algcrt:=0; i:=0; mistake:=false;                                     
   for j:=1 to len do                                                           
   begin ch:=braidstr[j];                                                       
      if ch='-' then sig:=-1                                                    
      else                                                                      
      begin if ch <>' ' then                                                    
         begin  i:=i+1;                                                         
            braid[i]:=sig*(ord(ch)-ord('0'));                                   
            algcrt:=algcrt+sig; sig:=1;                                         
            if abs(braid[i])>n-1 then mistake:=true;                            
         end;                                                                   
      end;                                                                      
   end; len:=i;                                                                 
   writeln('braid length = ',i:3);                                              
   for i:=1 to len do                                                           
   begin                                                                        
      r:=braid[i]; multiply(1,nfact);                                           
      writeln('started to multiply at position ',i:3);                          
   end;                                                                         
(* do partial wrapup and storage to get c right for next part *)                
  if K>0 then
  begin
   for s:= n-1 downto n-K do                                                      
   begin                                                                        
      blocklength:= factor[s];                                                  
      for g:= 1 to blocklength do d[g]:= c[g];                                  
      for t:= 0 to n-s-1 do                                                     
      begin                                                                     
         a:= t*factor[s+1]; b:= t*factor[s];                                    
         aa:=a+factor[s];  bb:= b+factor[s];                                    
         for g:= 1 to blocklength do                                            
         begin                                                                  
            if t>0 then 
            begin
               if c[g+a].mindeg>-1 then
               begin
                  c[g+b]:=plus(c[g+a],c[g+b]);
               end;
            end;
            c[g+bb]:= c[g+aa];                                                  
         end;{g}
         for j:= 2 to s do                                                      
         begin                                                                  
            bs:= a+j*factor[s];                                                 
            for g:= 1 to blocklength do                                         
            begin                                                               
               c[g]:= c[g+bs];                                                  
            end;                                                                
            reverse(s);                                                         
            for r:= s-1 downto s-j+1 do multiply(1,blocklength);                
            reverse(s);                                                         
            for g:= 1 to blocklength do                                         
            begin
               if c[g].mindeg>-1 then
               begin
                  c[g+bb]:=plus(c[g],c[g+bb]);
               end;
            end;{g}
         end; {s}                                                                  
      end;{t}                                                                      
      for g:= 1 to blocklength do                                               
      c[g]:= d[g];                                                              
   end;{s}
  end;                                                                         
   b:= factor[n-K];


   writeln;
   writeln('If you want to use the difference between this tangle and');
   writeln('the same tangle after a half-turn in the subsequent calculation, type');
   writeln('r   to subtract the reversed tangle,');
   writeln('t   to subtract the tangle conjugated by a half-twist braid.');
   writeln;
   writeln('Press RETURN  twice  to continue the calculation without alteration.');

   readln(ch1);
   if (ch1='r') or (ch1='t') then symmetrise;

   writeln;                                                            
(*   for g:= 1 to b do                                                            
   begin                                                                        
      if c[g]<>0 then writeln('entry ',g:3,' in c = ',c[g]:3);          
   end;                                                                         
*)
   for g:= 1+b to (K+1)*b do                                                  
   begin                                                                        
      d[g-b]:= c[g];                                                            
(*      if c[g]<>0 then writeln('entry ',g:3,' in c = ',c[g]:3);          
*)
   end;{g}                                                                         
   for g:= 1+b to nfact do c[g]:=zeropoly;                                             
   len1:=len;(* Remember first braid and its length *)
   braidstr1:=braidstr;

 (* now read in second braid and multiply by partial closure of 
the first, where the blocks holding the different degree information in T                     
   are taken one after another, controlled by the variable t  *)                
          
   writeln('ENTER SECOND BRAID, again as a braid on ',n:1,' strings. ');
                                                                
   readln(braidstr);                                                     
   len:=length(braidstr);                                                       
   sig:=1;  i:=0; mistake:=false;                                               
   for j:=1 to len do                                                           
   begin ch:=braidstr[j];                                                       
      if ch='-' then sig:=-1                                                    
      else                                                                      
      begin if ch <>' ' then                                                    
         begin  i:=i+1;                                                         
            braid[i]:=sig*(ord(ch)-ord('0'));                                   
            algcrt:=algcrt+sig; sig:=1;                                         
            if abs(braid[i])>n-1 then mistake:=true;                            
         end;                                                                   
      end;                                                                      
   end;                                                                         
   len:=i;                                                                      
                                                                                
   for t:= 0 to K do                                                          
   begin                                                                        
      writeln('Starting to deal with terms of degree ',t:1,' in T');
      writeln;
      test:= 1;                                                                 
      if t>0 then                                                               
      begin                                                                     
         test:= 0;                                                              
         for g:= 1 to b do    (* b= factor[n-K] *)                                
         begin                                                                  
            c[g]:= d[g+(t-1)*b];                                                
            if c[g].mindeg > -1 then test:= 1;                                         
         end;{g}                                                                   
         for g:= b+1 to nfact do c[g]:= zeropoly;
      end;                                                                      
      if test = 1 then                                                          
      begin                                                                     
         for   j:= 1 to len do                                                  
         begin                                                                  
            r:= braid[j];                                                       
            multiply(1,nfact);                                                  
            writeln('Multiplication done for crossing ',j+len1:1);
         end;                                                                   
         writeln;
                                                                                
      (* for g:= 1 to nfact do                                                  
         begin                                                                  
            if c[g] <>0 then writeln(t:3,g:5,c[g]:3);                   
         end;  *)                                                               
         wrapup;                                                                
         for j:= 0 to n-1 do  v[t+j]:=plus( v[t+j] , u[j]);
      end;                                                                      
   end;                                                                         
   fileout;                                              
   writeln('Do you want to calculate another example? Enter y/n.');
   readln(ch);
   if (ch='n') then halt;                 
   reinit;
   until false;
end.