# Maple procedure for calculating the multivariable Alexander polynomial
# of a link, presented as a braid on n strings.
# H.R.Morton and Julian Hodgson, Liverpool University, January 1996.

# Output procedure amended slightly August 1999, HRM, to ensure that
# knots are properly distinguished from links.

with(linalg):
readlib(factors):

braid := proc()

local m,n,r,funct,string_permutation,burau_sub,burau_matrix,burau_current,burau_product,for_subtraction,i,j,k,i1,x,y,count,plus_minus,i_th_string,offset,one_string_links,lone_number,trash,trash_count,trash_pointer,orbit_book,long_link,strings,flag,flag1,flag2,temp,long_link_pos,div,poly,result,result_factors,result_compare,result_final;
global X,T,characpoly;

# To use this procedure for calculating the Alexander polynomial of
# a closed braid, start Maple and then read in this file. 
# You can calculate the polynomial for the closure of a braid on
# n strings with braid word sigma_1 sigma_2^{-1} sigma_3, for example,
# by typing 
# braid(n,1,-2,3);



#       Note that the n-array T is global so that when the final 
#  expression  is  returned to  the user it  can be manipulated.
# When the procedure has finished, besides the information on screen,
# the array T[1]..T[n] of variables in the eventual output can be altered
# to allow for substitutions of your choice.

# The characteristic polynomial of the Burau matrix is also global,
# and is stored as characpoly, in terms of the variables in T, and a 
# variable X.



	if args[1]<2 
	then ERROR(`Not enough strings: Must be two or more!`)
	fi;

# n is number of strings

n:=args[1];

	for i from 2 to nargs do

		if abs(args[i]) > (n-1) 
		then ERROR(`String number out of range!`)
		fi

	od;

# m is number of permutations

m:=nargs-1:

	if m=0 
	then ERROR(`No permutations specified!`) 
	fi;

funct:=(i) -> i:
string_permutation:=vector(n,funct):

# Works out order of strings at other end of braid:

	for count from 1 to m do

	i:=args[count+1];
	i_th_string:=abs(i);
	
i1:=string_permutation[i_th_string];string_permutation[i_th_string]:=string_permutation[i_th_string+1];string_permutation[i_th_string+1]:=i1

	od;

# Create meridians T[1],...,T[n]

T:=array(1..n):

# Works out relations on the meridians T[i]:


# one_string_links stores `lone strings`

one_string_links:=array(1..n):lone_number:=0:

# trash stores strings processed (not lone)

trash:=array(1..n):trash_count:=0:
	
# orbit_book is a list of strings with orbits together with their orbits

orbit_book:=array(sparse,1..n,1..n+1):long_link:=0:           


strings:=0:

	while strings<>n do
	
	strings:=strings+1:

	flag:=0;
	flag1:=0;
	trash_pointer:=0;

		if trash_count<>0

		then    while flag=0 do

			trash_pointer:=trash_pointer+1;

				if trash_pointer=trash_count then flag:=1 fi;
				if trash[trash_pointer]=strings then 
flag:=1;flag1:=1 fi

			od;
		fi;

		if flag1=0 then


			if string_permutation[strings]=strings

			then    
lone_number:=lone_number+1;one_string_links[lone_number]:=strings;
				print(T[strings],`lone string`)

			else    flag:=0;
				temp:=strings;
				long_link_pos:=1;
				long_link:=long_link+1;
				orbit_book[long_link,long_link_pos]:=strings;

					while flag=0 do

					temp:=string_permutation[temp];

						if temp<>strings 

						then    
long_link_pos:=long_link_pos+1:
							
orbit_book[long_link,long_link_pos]:=temp:
							
trash_count:=trash_count+1:
							
trash[trash_count]:=temp

						else    flag:=1

						fi
					od
			fi

		fi
	
	od;

	for i from 1 to long_link do

	long_link_pos:=1;
	flag:=0;

#		lprint(`We have relations:`);

		while flag=0 do

		long_link_pos:=long_link_pos+1;

			if orbit_book[i,long_link_pos]=0

			then flag:=1

			else    
#   print(T[orbit_book[i,long_link_pos]]=T[orbit_book[i,1]]);
				
T[orbit_book[i,long_link_pos]]:=T[orbit_book[i,1]]
			fi
		od
	od;

# Defines 3 by 3 submatrix used forming burau matrices

burau_sub:=array(1..2):
burau_sub[2]:=matrix([[1,0,0],[r,-r,1],[0,0,1]]):
burau_sub[1]:=inverse(burau_sub[2]):

# Creates burau matrices

burau_matrix:=array(1..2,1..n-1):

	for j from 1 to 2 do
	
		for k from 1 to n-1 do

		burau_matrix[j,k]:=array(sparse,1..n,1..n);

			for x from 1 to n do
			burau_matrix[j,k][x,x]:=1
			od;

			if k=1 
			then offset:=2 
			else offset:=1 
			fi;
	
			for x from offset to 3 do

				for y from offset to 3 do
				
burau_matrix[j,k][k+x-2,k+y-2]:=burau_sub[j][x,y]
				od
	
			od
		od
	od;

#       Here we work out each individual burau matrix with the 
# corresponding meridian  as its variable  and keep a  running 
# total of the product.

burau_product:=array(identity,1..n-1,1..n-1):
for_subtraction:=array(identity,1..n-1,1..n-1):
string_permutation:=vector(n,funct):

	for count from 1 to m do

	i:=args[count+1];
	i_th_string:=abs(i);
	plus_minus := (3+i_th_string/i)/2 ;

		if signum(i)=1 

		then r:=T[string_permutation[i_th_string+1]]
		else r:=T[string_permutation[i_th_string]]

		fi;

	
	burau_current:=burau_matrix[plus_minus,i_th_string]:
	burau_current:=submatrix(burau_current,1..n-1,1..n-1):

# At this stage A is a matrix whose values are integers or functions of r.

	burau_current:=map(eval,burau_current):

# Now "map" has forced A to take on board that r:=T[i] for some i.

		if signum(i)=-1

		then    burau_current:=evalm(r*burau_current);
			for_subtraction:=evalm(r*for_subtraction)

		fi;

#       Here we multiply the matrix by the meridian whose value 
# it takes only if it is the inverse form. This keeps fractions
# out of the  calculation  rendering  it quicker. We also  keep 
# track of all the factors we have used for later.

#	print(sigma[i_th_string]);

#		if signum(i)=-1 
#		then print(inverse) 
#		fi;

#	print(burau_current);

	burau_product:=evalm(burau_current &* burau_product):

#	print(pi*sigma['i']);
#	print(burau_product);
	
	
i1:=string_permutation[i_th_string];string_permutation[i_th_string]:=string_permutation[i_th_string+1];string_permutation[i_th_string+1]:=i1
	
	od;



lprint(`This is the final composite matrix upon realisation of the above 
relations (if any):`);
print(burau_product);

div:=1;

	for i from 1 to n do

	div:=div*T[i]

	od;

div:=1-div;

lprint(`Subtracting the identity and taking the determinant gives:`);

poly:=evalm(burau_product-for_subtraction);
poly:=det(poly);


print(poly);
lprint(`which factorises to:`);

poly:=factor(poly);

print(poly);


lprint(`We divide by 1-T[1]...T[n] in factorised form:`);
div:=factor(div);
print(div);

result:=poly/div;
result:=factor(result);

knot:=1;
for i from 1 to n do if T[i]<>T[1] then knot:=0; fi;od;
	
        if knot=0


	then    lprint(`to give:`);
		
	else    lprint(`to give (here we give the polynomial having got rid of the factor 1-t which remains for a knot):`);
		result:=result*(1-T[1]);
		result:=simplify(result);
		
	fi;

print(result);

# This routine removes any factors of T[i] and should be removed if the
# extra time needed for it is significant

result_factors:=factors(result);
result_compare:=result_factors[1];
result_final:=result_factors[1];

	if result_compare<>result

	then 	i:=1;
		flag2:=0;

			while flag2=0 do

			# This routine returns flag=1 if 
			# the current factor is a power of T[i]

			flag:=0;
			flag1:=0;
			j:=0;

				while flag1=0 do

				j:=j+1;

					if result_factors[2][i][1]=T[j]

					then 	flag:=1;
						flag1:=1

					fi;

					if j=n
					then 	flag1:=1
					fi
 
				od;

				if flag=0 
				then 	
result_final:=result_final*(result_factors[2][i][1])^result_factors[2][i][2]
				fi;

			
result_compare:=result_compare*(result_factors[2][i][1])^result_factors[2][i][2];

				if expand(result_compare)=expand(result)
				then	flag2:=1
				fi;

			i:=i+1
			
			od

	fi;
characpoly:=det(evalm(X*burau_product-for_subtraction)):

lprint(`or without multiples of T[i]s:`);

RETURN(factor(expand(result_final)))


end;