program jones;                                                   

(* This program uses the representation of the braid group B_n on 2^n dimensional space, 
it corresponds with the use of the fundamental 2-dimensional module over the quantum group SU(2).  *)

(* Program written by H.R.Morton, January 1994. 
Modified April 1997, July 1998.
Laurent polynomials now handled in terms of s and not A.
Braid input routine modified May 2014, to use a string ended by 0, for more reliable
string handling on some compilers. File output disabled, again for portability.

The procedures extract and insert give the transformation from pointers to the Laurent polynomials
(which have type 'row'), but the main storage is a collection of 2^n pointers 
giving a vector in 2^n dimensional space with Laurent polynomial entries.

The 'multiply' routine operates on each basis vector in turn by the matrix determined by the braid
and collects its contribution to the weighted trace. The matrix itself has a lot of zero
entries, and it works as a composite of tensor products of identity matrices with 
a very simple 2x2 matrix. Its action is integrated with the numbering of the basis vectors.


*)

const  globaln=12;jmax=300; jmin=-300; globalexpn=4096; maxl=150;

(* globaln is the largest number of strings which will be handled, and 
globalexpn = 2^globaln. These can be altered if using a more or less powerful machine.
If they are increased then the input routine needs amendment to read in higher generators,
and the declaration of the size of entries in the braid array has to be changed. *)

type                                                                            
     row= record coefft:array[jmin..jmax] of integer;                           
                 maxdeg:integer;                                                
                 mindeg:integer                                                 
          end;                                                                  
     ptr1=^term;                                                                
     term=record coefft:integer;                                                
                 expon:jmin..jmax;                                              
                 nexterm:ptr1                                                   
          end;                                                                  
label startpoint;                                                                    
var s,g,gg,k,i,j,l,n,expn,len,


  r,sig:integer;
    algcrt:integer;                                                             
    mistake: boolean;                                                           
    t1:ptr1;                                                               
    z,w,poly,zerovect:row;                                                      
    braid:array[0..maxl] of -11..11;                                              
    braidstr: packed array [1..200] of char;
    anstr: packed array [1..200] of char;
    filename: packed array [1..160] of char;
    weight:array[0..globalexpn] of 0..globaln;                                              
    ch,ans:char;                                                                    
    c:array[0..globalexpn] of ptr1;                                                   
    outfile:text ;                                 

function min(a,b:integer):integer;
        begin if a<b then min:=a else
        min:=b;
        end;  {min}

function max(a,b:integer):integer;
    begin if a>b then max:=a else
    max:=b;
    end;  {max}

 function plus(a,b:row):row;                                                     
var c:row; i:integer;                                                           
   begin  if a.mindeg=-1000 then plus:=b else                                   
   if b.mindeg=-1000 then plus:=a else                                          
   begin                                                                        
      c := zerovect;                                                            
      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 power(num,index:integer):integer;                                      
var i,prod:integer;                                                             
begin prod:=1; for i:=1 to index do                                             
   prod:=prod*num;                                                              
   power:=prod;                                                                 
end;{power}                                                                     
                                                                                
                                                                                
function extract(t:ptr1):row;                                                   
var d		 : integer;
   z		 : row;
   t1		 : ptr1;                                                                     
begin t1	 :=t; z:=zerovect; d:=-1000;                                          
    if not(t1=nil) then z.mindeg :=t1^.expon;                                    
    while not(t1=nil) do                                                        
    begin d:=t1^.expon; z.coefft[d]:=t1^.coefft;                                
    t1:=t1^.nexterm;                                                            
    end; z.maxdeg:=d;                                                          
 extract:=z;
end;{extract}                                                                   

procedure stack(t:ptr1);                                                        
var t1:ptr1;                                                                    
begin while t<> nil do                                                          
   begin t1:=t; t:=t^.nexterm;                                                  
      dispose(t1);                                                            
   end;                                                                       
end;{stack}                                                                     

function  insert(z:row):ptr1;                                                    
var i,d		       : integer;                                                        
   t, t1,t2	       : ptr1;                                                           
begin       d:=-1000;
   if z.mindeg<=-1000 then t:=nil else                                       
   begin                                                                        
      new(t1);                                                                  
      t:=t1; t2:=t1;                                                         
      for i:=z.mindeg to z.maxdeg do if z.coefft[i]<>0 then                     
      begin t1^.coefft:=z.coefft[i];d:=i;    t1^.expon:=i;                      
         t2:=t1;                                                                
         new(t1);                                                               
         t2^.nexterm:=t1;                                                       
      end;                                                                      
      t2^.nexterm:=nil;    
      dispose(t1);                                                            
      if d=-1000 then t:=nil;
   end;
   insert:=t;
end;{insert}                                                                    

 procedure multiply(r:integer);                                                  
var i,j,k,r1,m,s,h: integer;                                                
    x,y:row;                                                                    
    t1,t2:ptr1;                                                                 
begin r1:=abs(r);                                                               
   for i:=0 to power(2,r1-1)-1 do                                               
   for j:=0 to power(2,n-r1-1)-1 do begin                                       
      s:=weight[g]-weight[i]-weight[j];                                         
      if s=1 then                                                               
      begin                                                                     
      h:=i+power(2,r1+1)*j+power(2,r1-1);                                       
         k:=i+power(2,r1+1)*j+power(2,r1);                                      
         t1:=c[h]; t2:=c[k];                                                    
         if (t1<>nil) or (t2<>nil) then                                         
         begin                                                                  
	    x:= extract(t1); y:= extract(t2);                                  
            if r>0 then                                                         
            begin if t2=nil then  c[h]:=nil                                     
               else                                                             
               begin z.mindeg  :=y.mindeg-1;z.maxdeg:=y.maxdeg-1;               
                  for m:=z.mindeg to z.maxdeg do                                
                  z.coefft[m]:=y.coefft[m+1];                                   
		  c[h]:= insert(z);                                                    
               end;                                                             
               z.maxdeg:=max(x.maxdeg-1,y.maxdeg);                              
               if t2=nil then z.mindeg:=x.mindeg-1 else
               if t1=nil then z.mindeg:=y.mindeg-2 else                         
               z.mindeg:=min(x.mindeg-1,y.mindeg-2);                            
               for m:=z.mindeg to z.maxdeg do                                   
               z.coefft[m]:=x.coefft[m+1]+y.coefft[m]-y.coefft[m+2];            
	       c[k]:= insert(z);                                                       
            end;                                                                
            if r<0 then                                                         
            begin if t1=nil then  c[k]:=nil    else                             
               begin z.mindeg  :=x.mindeg+1;z.maxdeg:=x.maxdeg+1;               
                  for m:=z.mindeg to z.maxdeg do                                
                  z.coefft[m]:=x.coefft[m-1];                                   
		  c[k]:= insert(z);                                                    
               end;                                                             
               z.maxdeg:=max(y.maxdeg+1,x.maxdeg+2);                            
               if t2=nil then z.mindeg:=x.mindeg else                           
               if t1=nil then z.mindeg:=y.mindeg+1 else                         
               z.mindeg:=min(x.mindeg,y.mindeg+1);                              
               for m:=z.mindeg to z.maxdeg do                                   
               z.coefft[m]:=y.coefft[m-1]+x.coefft[m]-x.coefft[m-2];            
	       c[h]:= insert(z);                                                       
	    end;            
                    stack(t1);stack(t2);
         end;                                                                   
      end;                                                                      
   end;                                                                         
end;{multiply} 
                                                                
procedure init;                                                                 
var i,q : integer;                                                            

begin                                                                           
   //rewrite(outfile,'jones.dat');
   weight[0]:=0; weight[1]:=1;
   expn:=power(2,n);
   for i:=1 to n-1 do                                                           
   begin for q:= 1 to power(2,i) do                                             
      weight[q+power(2,i)-1]:=weight[q-1]+1;                                    
   end;                                                                         
   zerovect.mindeg:=-1000;zerovect.maxdeg:=-1000;                               
   for i:=jmin to jmax do zerovect.coefft[i]:=0;                                
   poly:=zerovect;                                                              
   z:=zerovect;                                                                 

end;{init}    
                                                                  
 procedure jonesprint(poly:row;fn:integer);                                      
var sdeg,k,t:integer;                                                           
begin                                
writeln;
writeln('The Jones polynomial, as a function of s = t^(1/2), when multiplied by s+s^(-1),');
writeln('is the following Laurent polynomial P in s;');
writeln;
write(' P = ');                                           
if poly.mindeg=-1000 then write(0:4)                             
 else for k:=poly.mindeg to poly.maxdeg do                                    

    begin t:=poly.coefft[k]; if t>0 then                                       
          begin sdeg:=k-algcrt;                                               
          if fn=1 then writeln(outfile,'+',t,'*s^(',sdeg:2, ')')                        
          else 
	  write(' +',t:2,'*s^(',sdeg:2,')');                                
          end else if t<0 then
          begin sdeg:=k-algcrt;
          if fn=1 then writeln(outfile,t,'*s^(',sdeg:2, ')')
          else
          write(' ',t:2,'*s^(',sdeg:2, ')');
          end;                                                                 
    end;                                                                        
    writeln;                                                            
end; {jonesprint}

                                                                           
procedure clearscreen;                                                          
begin  writeln('press return');                         
      readln(anstr);                                                      
end;                                                                            
                                                    
begin; {main program}   
     //filename:='braidout';rewrite(outfile,filename);
   n:=globaln;expn:=globalexpn;
writeln;   
writeln('This program calculates the Jones polynomial of a closed braid on at most twelve strings');
writeln('It uses a 2^n dimensional representation of the n-string braid group.');

    
    repeat writeln;                                                     
    
startpoint: 
 writeln('Give braid to be calculated, on at most', globaln:3,' strings');
   writeln('Use 1,2,.. for the standard braid generators and -1,-2,... for their inverses.');
   writeln('Write the braid word without any separators, and finish with 0');
   writeln;
   writeln('for example 1-232-3-10');
   writeln('(For generators 10 and 11 use a and b.)');
   writeln;
   writeln('Type * to exit the program');

             
    readln(braidstr);                                                     

       if not(braidstr[1]='*') then                                                   
       begin

	    writeln('The number of strings to be used in this braid is currently',n:3);
	  writeln;
       writeln('If you want to change this number type y');
       writeln('else press any key, and return, to continue with' ,n:3, ' strings');
       ans:='N';
       readln(ans);
       while ans='y' do
       begin writeln('Give the number of strings in the braid');
	  readln(n);
	  if (n>globaln) or (n<2) then
	     writeln('The number of strings should be >1 and at most',globaln:3)
	  else ans:='N';
       end;
     
        init;


	  i:=1;
       while braidstr[i]<>'0' do
       begin    i:=i+1;end;
       len:=i-1;
(* This is  a way to convert the braidstring input to an array of integers *)

      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 i:=i+1;                                                         
              if ch='a' then braid[i]:=sig*10 else
              if ch='b' then braid[i]:=sig*11 else
              braid[i]:=sig*(ord(ch)-ord('0'));                                 
              algcrt:=algcrt+sig; sig:=1;                                       
               if abs(braid[i])>n-1 then mistake:=true;                         
          end;                                                                  
      end; len:=i;                                                              
      if (len>maxl-15) then writeln('DANGER..braid too long!!');writeln;         
      if mistake=true then                                                      
      begin writeln('wrong elements in braid');writeln;goto startpoint;               
      end;                                                                      
      writeln('Starting to calculate');
      for g:=0 to expn-1 do                                                     
      begin z := zerovect;

     if (g+1) mod 100 =0 then  writeln('reached row  ',g+1:1, ' of ', expn:1);
     //This gives a check on progress if the computation is slow. Uncomment if needed.

       for j :=0 to expn-1 do c[j]:=nil;                                    
         new(c[g]);                                                             
         (* Check(Poly) *) ;                                                           
         c[g]^.coefft:=1;                                                       
         c[g]^.expon:=0;                                                        
         c[g]^.nexterm:=nil;                                                    
         for i:=1 to len do begin                                               
             multiply(braid[i]);(* Check(Poly) *)  end;                                
         if c[g]<>nil then                                                      
         begin                                                                  
	    z:= extract(c[g]);                                                         
            w:=zerovect;                                                        
            l:=2*weight[g]-n;                                                   
            for k:=z.mindeg to z.maxdeg do                                      
                w.coefft[k+l]:=z.coefft[k];                                   
            w.mindeg:=z.mindeg+l;                                             
            w.maxdeg:=z.maxdeg+l;                                             
            poly:=plus(poly,w);                                                 

            for gg:= 0 to expn-1 do
            stack(c[gg]);    
(* This clears the memory after each row of the matrix is dealt with *)
         end;                                                                   
                                                                                
      end;                                                                      

   jonesprint(poly,0);                                                          
   end;                                                                         
   until braidstr[1]='*';                                                          
end.