# galois
# by Ravi Vakil
# Fri. July 5, 2002, 11 am

#
#   Frank Sottile is altering this mainly for more informative output
#    circa Jan 2018

#
# (See notes from roughly June 5.)
#
# To use:  read galois
#		biggalois(3,7)
#
# The goal of this program is:
# (i)   to implement the geolr in string notation
# (ii)  to calculate solutions to Schubert problems 
# (iii) to make a fancier way of asking Schubert questions 
# (iv)  to determine when I know that the monodromy group
#       is full or alternating  
# (v)   to check the monodromy for an entire Grassmannian 
#
#----------------------------------------------------------------------
# TABLE OF CONTENTS
#----------------------------------------------------------------------
#
# Table of contents
# Notes
# Variables
# Initial commands
# Useful subroutines
# 	initpos
#	wempty
#	subtoseq
#	seqtosub
#	start
#	checkergame
#	pcheckergame
#	finishgame
#	gcheckergame
#	galoisstart
#	gal
#	biggalois
#	prprint
#	prprint3
#	addtotlist
#	addtolist2
#	addtolist3
#	checkzerocycle
#	btournamentstart
#	bbournament
#	tournament
#	findinlist
#	filsub
#	compare
#	sortorder
#	trim
#	trim3
#	checktlist
#	lr
#
# Next:  schubert
#
#----------------------------------------------------------------------
# NOTES
#----------------------------------------------------------------------
# Programs to write:
#
# codimension(ul,vl,w,vr,ur)
#
#
#
#
#----------------------------------------------------------------------
# VARIABLES
#----------------------------------------------------------------------
# n and k are fixed
# T is -1
# [ul (vl { w }  vr)  ur ]
# 	ul and ur are lists of 0's and 1's
# 	vl and vr are lists of 0's, 1's, and T's
# 	w is a list of 0, 1, or nothing (later change to a number)
# slist1 = [alpha,beta] for "better tournament"
# slist2 = corresponding Schubert number

#----------------------------------------------------------------------
# INITIAL COMMANDS
#----------------------------------------------------------------------
#restart:
with(combinat):
slist1 := []: slist2 := []:  # list 1 is Schubert problems; 2 is answer
tlist1 := []:    # list of Schubert problems mid-tournament
tlist2 := []:    # the answer
tlist3 := []:    # +1 if known to be A_n or S_n, +2 if can't tell
		# because of n/n issue, +3 if can't tell because of 6/6
		# issue
		# Change!  0 if known to be A_n or S_n,
		# positive integer k if there are k issues en route
		# (included, because we hope this answer is 1 in general,
		# i.e. once I find an issue, I know it is the only one)
		# infinity if there is a Mathieu issue!
# Fix n and k here.

#
#
#

#----------------------------------------------------------------------


#----------------------------------------------------------------------
# SUBROUTINES
#----------------------------------------------------------------------


#----------------------------------------------------------------------
# initpos(aa,bb)
#----------------------------------------------------------------------
# starts with two sequences of 0's and 1's (k 1's, n-k 0's)
# and produces an input list ul,vl,w,vr, ur
initpos := proc(aa,bb)
local output;
	return[ aa[1..(nops(aa)-1)],[],[aa[nops(aa)]], [],bb];
end:

#----------------------------------------------------------------------
# wempty(ul,vl,vr,ur)
#----------------------------------------------------------------------
# produces acceptable output for lr

wempty := proc(ul, vl, vr, ur)
local output, w, ul1, vl1, vr1, ur1;

	if ul = [] then
		output := [true, vl];
		return output;
	else
#	vr := vl;vl := [];w := [ul[nops(ul)]];ul := ul[1..(nops(ul)-1)];
#	output := [false, [ul,vl,w,vr,ur]]:
		output := [false, [ul[1..(nops(ul)-1)], [], 
			[ul[nops(ul)]], vl, ur]];
		return output;
	fi;
end:

#
#----------------------------------------------------------------------
# subtoseq(a,n)
#----------------------------------------------------------------------
# turns a subset a of {1,...,n} into a sequence of 0's and 1's
subtoseq := proc(a,n)
local output,i:
output := []:
for i from 1 to n do output:= [op(output),0]: od:
for i from 1 to nops(a) do output[a[i]] := 1: od:
return copy(output):
end:

#----------------------------------------------------------------------
# seqtosub(a)
#----------------------------------------------------------------------
# turns a sequence of 0's and 1's into a subset of {1, ..., n}
seqtosub := proc(a)
local output, i:
output := {}:
for i from 1 to nops(a) do
	if a[i] = 1 then output := output union {i}: fi:
od: 
return convert(output,list):
end:

#----------------------------------------------------------------------
# start(i1,i2)
#----------------------------------------------------------------------
# input:  2 sequences
start := proc(i1,i2)
local a1,a2,flag,i,nn,kk:
nn := nops(i1):
a1 := [op( seqtosub(i1))]:
a2 := [op( seqtosub(i2))]:
kk := nops(a1):
flag := true:
if kk <> nops(a2) or nn <> nops(i2) then 
	print(`ERROR:  not in same Grassmannian!`):
	return false:
fi:
for i from 1 to kk do
	if a1[i] + a2[kk+1-i] < nn+1 then flag := false: fi:
od:
if flag=false then return false: fi:
return initpos(i1,i2):
end:

#----------------------------------------------------------------------
# checkergame(i1,i2)
#----------------------------------------------------------------------
# input:  2 sequences; output a sequence of sequences
checkergame := proc(i1,i2)
local board:
board := start(i1,i2):
if board = false then
	return false:
else
	return finishgame(board):
fi:
end:

#----------------------------------------------------------------------
# pcheckergame(i1,i2) = "pretty checkergame"
#----------------------------------------------------------------------
# input:  2 sequences; output a sequence of sequences
pcheckergame := proc(i1,i2)
local board:
board := start(i1,i2):
if board = false then
	print(`No solution.`):
else
	print(finishgame(board)):
fi:
end:


#----------------------------------------------------------------------
# finishgame(input)
#----------------------------------------------------------------------
#
finishgame := proc(input)
local nx:
nx := lr(input):
if nx[1] = true then 
	return [nx[2]]:
elif nops(nx)=2 then
	return finishgame(nx[2]):
else
	return [op(finishgame(nx[2])),op(finishgame(nx[3]))]
fi:
end:

#----------------------------------------------------------------------
# gcheckergame(input)
# This is a variant on finishgame
# It returns either (i) a sequence (output)
# or (ii) a triplet:  state, and two possible next states
gcheckergame := proc(input)
local nx:
nx := lr(input):
if nx[1] = true then
	return [nx[2]]:
elif nops(nx)=2 then
	return gcheckergame(nx[2]):
else
#print(`In gcheckergame, with input`,input,` returning `,input,nx[2],nx[3]):
	return [input, nx[2],nx[3]]
fi:
end:


#----------------------------------------------------------------------
# galoisstart(al,be)
#----------------------------------------------------------------------
# This is the right way to start the galois calculation,
# where we have only Schubert conditions, and 
# where we don't worry about sorting
galoisstart := proc(al,be)
local alpha, beta,i:
i := sortorder(al,be):
alpha := i[1]: beta := i[2]:
i:= gal(alpha,beta,[],[],[]):
print(`Number = `,i[1]):
if i[2]=0 then 
	print(`Galois group is at least alternating.`):
elif i[2]>0 then
	print(`We have a problem!`):
	if i[2] = infinity then 
		print(`Possible Mathieu issue!`):
	fi:
fi:
end:


#----------------------------------------------------------------------
# gal(al,be,ga,de,ep)
# see btournament for explanation of alpha, beta
# gamma are multiples
# delta are the mid-game checker calculation
# epsilon are the lists of two next boards
# returns [number, answer to the galois question]
##
##  This is where a given state is attempted to degenerate into all possible substates
##
gal := proc(al,be,ga,de,ep)
global tlist1, tlist2, tlist3:  
local i, answer, j, answer1, answer2, k1, k2, temp,
 gt,dt,et,at,bt, number, q, r, gcg, a1, b1, g1, d1, e1, a2, b2, g2, d2, e2:

answer := infinity:

#print(`Starting gal al=`,al,`be=`,be,`ga=`,ga,`de=`,de,`ep=`,ep):

if nops(al)=1 and al[1]=1 and nops(ga)=0 then 
#@print(`Gal input`,al,be,ga,de,ep,`output`, [checkzerocycle(be[1]),1]):
	return [checkzerocycle(be[1]),0]:
fi:


i := findinlist([al,be,ga,de,ep],tlist1):
if i[1] = true then
#@print(`Gal input`,al,be,ga,de,ep,`output`, [tlist2[i[2]],tlist3[i[2]]]):

	return [tlist2[i[2]], tlist3[i[2]]]:
fi:

# First, try degenerating one of the games-in-progress
for i from 1 to nops(ga) do
	gt := copy(ga):  dt := copy(de):  et := copy(ep):
	gt[i] := gt[i]-1:
	j := trim3(gt,dt,et):  gt := j[1]:  dt :=j[2]: et := j[3]:

	k1 := ep[i][1]:
	k2 := ep[i][2]:

	# Do first degeneration, put the answer in answer1
	g1 := copy(gt): d1 := copy(dt): e1 := copy(et):

	gcg := gcheckergame(k1):
	if nops(gcg)=1 then
		# add to (a1, b1)

		j := addtolist2(al, be, gcg[1]):
		a1 := j[1]: b1 := j[2]:

		g1 := gt: d1 := dt : e1 := et:
	else
		# add to (g1, d1, e1)
		j := addtolist3(g1, d1, e1, gcg[1], [gcg[2],gcg[3]]):
		g1 := j[1]:  d1 := j[2]:  e1 := j[3]:
		a1 := al: b1 := be:
	fi:

	answer1 := gal(a1, b1, g1, d1, e1):

	# Do second degeneration, put the answer in answer2
	g2 := copy(gt): d2 := copy(dt): e2 := copy(et):

	gcg := gcheckergame(k2):
	if nops(gcg)=1 then
		# add to (a2, b2)
		j := addtolist2(al, be, gcg[1]):
		a2 := j[1]: b2 := j[2]:
		g2 := gt: d2 := dt : e2 := et:
	else
		# add to (g2, d2, e2)

		j := addtolist3(g2, d2, e2, gcg[1], [gcg[2],gcg[3]]):
		g2 := j[1]:  d2 := j[2]:  e2 := j[3]:
		a2 := al: b2 := be:
	fi:

	answer2 := gal(a2, b2, g2, d2, e2):

	number := answer1[1] + answer2[1]:
	if answer1[1] <> answer2[1] then
		answer := min(answer, answer1[2]+ answer2[2]):
	else
		if answer1[1]=0 or answer1[1]=1 then
			answer := 0:
		elif answer1[1]=6 then
			answer := min(answer, infinity):
		else
			answer := min(answer, answer1[2]+answer2[2]+1):
		fi:
	fi:

	if answer = 0 then 
		temp := addtotlist(al,be,ga,de,ep,number,0):
		#@print(`Gal input`,al,be,ga,de,ep,`output`,[number,1]):
		return [number, 0]:
	fi:

od:



# Second, try intersecting two Schubert classes
for i from 1 to nops(al) do for j from i to nops(al) do
if i <> j or al[i]>1 then

	at := copy(al):  bt := copy(be):  
	gt := copy(ga):  dt := copy(de): et := copy(ep):




	at[i] := at[i]-1:  at[j] := at[j]-1:
	r := start(bt[i], bt[j]):
	if r = false then 
		#@print(`Gal input`,al,be,ga,de,ep,`output`,[0,1]):
		return [0,0]: 
	fi:

	q := trim(at, bt): at := q[1]: bt := q[2]:

	gcg := gcheckergame(r):



	if nops(gcg)=1 then
		# add to (at, bt)
		q := addtolist2(at, bt, gcg[1]):
		at := q[1]: bt := q[2]:
	else
		# add to (gt, dt, et)

		q := addtolist3(gt, dt, et, gcg[1], [gcg[2],gcg[3]]):

		gt := q[1]:  dt := q[2]:  et := q[3]:
	fi:


	answer1 := gal(at,bt,gt,dt,et):
	number := answer1[1]:
	answer := min(answer, answer1[2]):

	if answer = 0 then 
		#@print(`Gal input`,al,be,ga,de,ep,`output`,[number,1]):
		temp := addtotlist(al,be,ga,de,ep,number,0):
		return [number, 0]: 
	fi:

fi:
od: od:

temp := addtotlist(al,be,ga,de,ep,number,answer):
#@print(`Gal input`,al,be,ga,de,ep,`output`,[number,answer]):
return [number, answer]:

end:


#----------------------------------------------------------------------
# biggalois(k,n)
#----------------------------------------------------------------------
# verifies that the monodromy is at least alternating for G(k,n)
# checks only partitions not containing the first row or first column.
# In Vakil's original code, this was called biggaloisnew
#
# goal:  make a list of k-subsets of [n], along with
# codimensions
# sum through partitions of n
# for each partition, go through possible lists
# trim list of 0's
# then check gal
biggalois := proc(k,n)
local i, templ, tt,  codim, t, j, st, en, nu, p,  ap, bb, 
 doneflag, newt,  h, b, c, s, gal1, gal2, tgal1, tgal2, tgal, booger,
	booger2, booger3, numProbs:

 numProbs :=[0,0,0]:  # Num zero, one, > one

#print(`Got here1`):
 
tt := choose(n-2,k-1):				# change from original
booger := true:  # Answers question:  do I know that all groups are
	# at least alternating?
booger2 := []:booger3 := []:  # The list of counterexamples
templ := []: codim := []:

# change from original (next paragraph)
for i from 1 to binomial(n-2,k-1) do		
#	ap := tt[i]:   	
	ap := []:
	for j from 1 to k-1 do:
		ap := [ op(ap), tt[i][j]+1]:
	od:
	ap := [ op(ap), n]:
	tt[i] := ap:

	templ := [op(templ),subtoseq( { op(tt[i])},n)]:
	t := n*k - k*(k-1)/2:
	for j from 1 to k do
		t := t - tt[i][j]:
	od:
	codim := [op(codim), t]:

od:

#print(`Got here3`):


for i from 1 to binomial(n-2,k-1)-1 do		# change from original
	for j from i+1 to binomial(n-2,k-1) do	# change from original
		if codim[i] > codim[j] then
			t:= codim[i]: codim[i] := codim[j] : codim[j] :=t:
			t:= templ[i]: templ[i] := templ[j]: templ[j] := t:
		fi:
	od:
#print( templ[i],codim[i]):
od:


st := array(1..(k*(n-k))):  
	# This says when a certain codimension starts in the list.
en := array(1..(k*(n-k))):  # This says when it ends.
nu := array(1..(k*(n-k))):  # This is en-st+1


#print(`Got here2`):

for i from 2 to binomial(n-2,k-1) do	# change from original


	if codim[i]>codim[i-1] then 

		st[codim[i]]:= i: 
		if codim[i-1]>0 then
			en[codim[i-1]] := i-1: 
		fi:
	fi:

od:

# change from original [entire paragraph]
en[ (k-1)*(n-k-1)] := binomial(n-2,k-1):
for i from (k-1)*(n-k-1)+1 to k*(n-k) do
	en[ i] := binomial(n-2,k-1):
	st[ i] := en[i] +1:
od:
# en[ k*(n-k)] := binomial(n,k):


for i from  1 to k*(n-k) do
	nu[i] := en[i] - st[i] + 1:
od:

#print(`Got here7 nu=`,nu):

bb := array(1..k*(n-k)):

for p from numbpart(k*(n-k)) by -1 to 1  do
	ap := partition(k*(n-k))[p]:

if ap[ nops(ap)] <= (k-1)*(n-k-1) then  # change from original
	
#print(`Got here6`,ap):

	# We write ap in more convenient notation:
	# bb[1] 1's etc.
	for i from 1 to k*(n-k) do
		bb[i] := 0:
	od:
	for i from 1 to nops(ap) do
		bb[ap[i]] := bb[ap[i]]+1:
	od:

# begin last.tmp
# (This should later be indented a bit.)


# Given:  some arrays of size k(n-k):
# a[i] - the number of i's in the partition (so \sum i a[i] = k(n-k))
# 	(now bb)
# b[i] - list of arrays (each running from st[i] to en[i])

# We have a list template of length n choose k (of sequences)
# the first one will never get used!
# newt is ntemplate

b := array[1..(k*(n-k))]:  # This goes before partition loop too.

newt := []:
for i from 1 to binomial(n,k) do
	newt := [op(newt),0]:
od:

for i from 1 to (k-1)*(n-k-1) do  				# change from original

#print(`Got here4`,i,nu[i],bb[i]):

	s := choose( nu[i]+bb[i]-1, nu[i]-1):
#print(`Got here5`):

	b[i] := array[1..nops(s)]:

	# recall:  b[ codimension in question ] [number of multiset] [class of this codimension]
	# variables		i			j			h
	for j from 1 to nops(s) do

		b[i][j] := array( st[i]..en[i] ):

		for h from 2 to (nu[i]-1) do
			b[i][j][h+st[i]-1] := s[j][h]-s[j][h-1]-1:
		od:


		if nu[i]=1 then 
			b[i][j][st[i]] := bb[i]:
		fi:


		if nu[i]>1 then 
			b[i][j][st[i]] := s[j][1]-1:  
			b[i][j][en[i]] := nu[i]+bb[i] - s[j][nu[i]-1]-1:
		fi:


	od:




od:


#print(`Analyzing`,ap):

c := array[1..(k*(n-k))]:

newt := [0]:
for i from 1 to k*(n-k) do
	c[i] := 1:   # 1 <= c[i] <= nops(b[i])
	for h from st[i] to en[i] do
		newt := [op(newt), b[i][1][h]]:
	od:
od:

doneflag := false:
while doneflag = false do

	# Here, I can call gal
	# Instead, I do the following temporary test.
#	print(newt):

	#tgal = trimmed galois input
	gal1 := copy(newt):  gal2 := copy(templ):
	tgal := trim(gal1,gal2):
	tgal1 := tgal[1]: tgal2 := tgal[2]:
#	print(`In particular`,tgal[1],tgal[2]):


# start of section grafted from galoisstart

tgal := sortorder(tgal1,tgal2):
tgal1 := tgal[1]: tgal2 := tgal[2]:
tgal:= gal(tgal1,tgal2,[],[],[]):


     if tgal[1]=0 then numProbs[1]:=numProbs[1]+1: end if:
     if tgal[1]=1 then numProbs[2]:=numProbs[2]+1: end if:
     if tgal[1]>1 then numProbs[3]:=numProbs[3]+1: end if:
if tgal[2]>0 and tgal[2] < infinity then 
#	printf("We have a problem: %A %A %A\n",tgal1,tgal2,tgal[2]):
	booger := false: 
	booger2 := [op(booger2), [ tgal1,tgal2,tgal[2],tgal[1]]]:
		# booger2 is a _triple_
elif tgal[2] =infinity then
#	print(`Possible Mathieu issue!,tgal1,tgal2`):
	booger := false:
	booger3 := [op(booger3), [ tgal1,tgal2]]:
fi:

#	end of piece grafted from galoisstart


	# Next, find the first c from the end we can update
	j := 0:
	for i from 1 to k*(n-k) do
		#nops(b[i]) doesn't work because b[i] is an array
		if c[i] < binomial( nu[i]+bb[i]-1, nu[i]-1) then 
			j := i:
		fi:
	od:
	if j = 0 then
		doneflag := true:
	else

		c[j] := c[j]+1:
		for i from j to k*(n-k) do
			if i>j then c[i] := 1: fi:
			for h from st[i] to en[i] do
				newt[h]  := b[i][c[i]][h] :
			od:
		od:
	fi:
od:


#end last.tmp

fi:  # change from original

	
od: # end loop over number of partitions

if booger = false then
	printf("There were %1d Schubert problems that I",nops(booger2)+nops(booger3)):
	printf(" could not tell were at least alternating:\n"):
	prprint3(booger2): 
	if booger3 <> [] then
		print(`The following cases have Mathieu issues!`):
		prprint(booger3):
	fi:

else
	printf("All cases okay!\n"):
fi:

printf(" Number of problems= %d  Number with [zero, one, at least 2] solutions = [%d, %d, %d].\n",
         numProbs[1]+numProbs[2]+numProbs[3],numProbs[1],numProbs[2],numProbs[3]):

end:

#----------------------------------------------------------------------
# prprint(l)
#----------------------------------------------------------------------
# prettyprints output
# (invoked in biggalois)
prprint := proc(l)
local i,li, lold, lnew, j, isnew, n:

lold := []:lnew := []: 

n := nops(l[1][2][1]):

for i from 1 to nops(l) do
	isnew := true:
	for j from 1 to nops(l[i][1]) do
		if l[i][2][j][1] = 1 or l[i][2][j][n] = 0 then
			isnew := false:
		fi:
	od:
	if isnew= true then
		lnew := [op(lnew), l[i]]:
	else
		lold := [op(lold), l[i]]:

	fi:
od:

print(`Old:`):
for i from 1 to nops( lold) do
	print(`---`):	
	for j from 1 to nops(lold[i][1]) do
		print(lold[i][1][j],`times`,seqtosub(lold[i][2][j])):
	od:
od:

print(`New:`):
for i from 1 to nops( lnew) do
	print(`---`):	
	for j from 1 to nops(lnew[i][1]) do
		print(lnew[i][1][j],`times`,seqtosub(lnew[i][2][j])):
	od:
od:

end:

#----------------------------------------------------------------------
# prprint3(l)
#----------------------------------------------------------------------
# prettyprints output
# (invoked in biggalois)
prprint3 := proc(l)
local i,li, lold, lnew, j, isnew, n:

lold := []:lnew := []: 

n := nops(l[1][2][1]):

for i from 1 to nops(l) do
	isnew := true:
	for j from 1 to nops(l[i][1]) do
		if l[i][2][j][1] = 1 or l[i][2][j][n] = 0 then
			isnew := false:
		fi:
	od:
	if isnew= true then
		lnew := [op(lnew), l[i]]:
	else
		lold := [op(lold), l[i]]:
	fi:
od:

#print(`Old:`):
#for i from 1 to nops( lold) do
#	print(`---`):	
#	for j from 1 to nops(lold[i][1]) do
#		print(lold[i][1][j],`times`,seqtosub(lold[i][2][j])):
#	od:
#	print(lold[i][3],`issue(s)`):
#od:

#print(`New:`):
for i from 1 to nops( lnew) do
	printf(" %5d = ",lnew[i][4]):
	for j from 1 to nops(lnew[i][1]) do
	  printf("%A^%d ",seqtosub(lnew[i][2][j]), lnew[i][1][j]):
	od:
	printf("    \n"):
#	print(lnew[i][3],`issue(s)`):
od:




end:

#----------------------------------------------------------------------
# addtotlist(a,b,g,d,e,number,answer)
#----------------------------------------------------------------------
# adds abgde to tlist
addtotlist := proc(a,b,g,d,e,number,answer)
global tlist1, tlist2, tlist3:
local place:
place := findinlist( [a,b,g,d,e], tlist1):

if place[1] = false then  # place[1] _should_ be false!

	tlist1 := [ op(tlist1[1..(place[2]-1)]), [a,b,g,d,e], 
		op(tlist1[place[2]..nops(tlist1)])]:
	tlist2 := [ op(tlist2[1..(place[2]-1)]), number, 
		op(tlist2[place[2]..nops(tlist2)])]:
	tlist3 := [ op(tlist3[1..(place[2]-1)]), answer,
		op(tlist3[place[2]..nops(tlist3)])]:

fi:
end:



#----------------------------------------------------------------------
# addtolist2(alpha, beta, c)
#----------------------------------------------------------------------
# adds c to list
# returns [newalpha,newbeta]
addtolist2 := proc(alpha,beta,c)
local newalpha, newbeta, place:

newalpha := copy(alpha):
newbeta := copy(beta):
place := findinlist( c, beta):

if place[1] = true then
	newalpha[ place[2] ] := newalpha[ place[2] ] +1:
else
	newalpha := [ op(newalpha[1..(place[2]-1)]), 1, 
		op(newalpha[place[2]..nops(newalpha)])]:
	newbeta := [ op(newbeta[1..(place[2]-1)]), c,
		op(newbeta[place[2]..nops(newbeta)])]:
fi:
return [newalpha,newbeta]:
end:

#----------------------------------------------------------------------
# addtolist3(alpha, beta, gamma, c, d)
# adds c and d to list
# returns [newalpha,newbeta, newgamma]
# (usually used for gamma, delta, epsilon)
# Before calling, make sure to type 
addtolist3 := proc(alpha,beta,gamma,c,d)
local newalpha, newbeta, newgamma, place:

newalpha := copy(alpha):
newbeta := copy(beta):
newgamma := copy(gamma):
place := findinlist(c,beta):
if place[1] = true then
	newalpha[ place[2] ] := newalpha[ place[2] ] +1:
else
	newalpha := [ op(newalpha[1..(place[2]-1)]), 1, 
		op(newalpha[place[2]..nops(newalpha)])]:
	newbeta := [ op(newbeta[1..(place[2]-1)]), c, 
		op(newbeta[place[2]..nops(newbeta)])]:
	newgamma := [ op(newgamma[1..(place[2]-1)]), d,
		op(newgamma[place[2]..nops(newgamma)])]:
fi:
return [newalpha,newbeta,newgamma]:
end:



#----------------------------------------------------------------------
# checkzerocycle(l)
#----------------------------------------------------------------------
# checks is a list is a bunch of 0's followed by a bunch of 1's
# if so, returns 1
# if not, returns 0
checkzerocycle := proc(l)
local i:
if nops(l)<2 then return 1: fi:
for i from 1 to nops(l)-1 do
	if l[i] = 0 and l[i+1]=1  then 
		return 0: 
		print(`Non-point class found!`):
	fi:
od:
return 1:
end:

#----------------------------------------------------------------------
# btournamentstart(al,be)
#----------------------------------------------------------------------
# This is the right way to start the btournament
#   when we don't worry about sorting in advance
btournamentstart := proc(al,be)
local alpha, beta,i:
i := sortorder(al,be):
alpha := i[1]: beta := i[2]:
return btournament(alpha,beta):
end:

#----------------------------------------------------------------------
# btournament(alpha,beta)  "better tournament"
#----------------------------------------------------------------------
# beta[1], ..., beta[n] are the sequences
# alpha[1], ..., alpha[n] are the multiplicities
btournament := proc(alpha,beta)
global slist1, slist2:
local i, at, bt, s1, s2, lranswer, a1, b1, place, answer, loc:

if nops(alpha)=1 and alpha[1]=1 then return checkzerocycle(beta[1]): fi:

i := findinlist([alpha,beta],slist1):
if i[1] = true then
	return slist2[ i[2] ]:
fi:

at := copy(alpha): bt := copy(beta):
if nops(alpha)>1 then
	s1 := beta[1]: s2 := beta[2]:
	at[1] := at[1]-1:  at[2] := at[2] - 1:

#!print(`We are intersecting the first two.  The remaining alpha, beta =`,at,bt):

else

	s1 := beta[1]:  s2 := beta[1]:
	at[1] := at[1]-2:

#!print(`We are intersecting two of the first.  The remaining alpha, beta =`,at,bt):

fi:
# at and bt form the truncated list
i := trim(at, bt): at := i[1]: bt := i[2]:

#!print(`After truncation, we have at, bt=`,at,bt):

lranswer := checkergame(s1,s2):

#!print(`checkergame result:`,lranswer):

if lranswer = false then return 0: fi:

answer := 0:

for i from 1 to nops(lranswer) do

	loc := addtolist2( at, bt, lranswer[i]):
	a1 := loc[1]: b1 := loc[2]:

#!print(`We are ready for the next btournament, with a1, b1=`,a1,b1):

	answer := answer + btournament(a1,b1):
od:


place := findinlist([alpha,beta],slist1)[2]:

slist1 := [ op(slist1[1..(place-1)]), [alpha,beta], 
		op(slist1[place..nops(slist1)])]:
slist2 := [ op(slist2[1..(place-1)]), answer, 
		op(slist2[place..nops(slist2)])]:

#!print(`slists:`,slist1,slist2):

return answer:

end:

#----------------------------------------------------------------------
# tournament(c)
#----------------------------------------------------------------------
# c is a list of sequences
tournament := proc(c)
local newc, answer, i, a:

if nops(c) = 1 then 
	return checkzerocycle(c[1]):
else

	a := checkergame(c[1],c[2]):

	if a = false then return 0: fi:
	newc := op(c[3..nops(c)]):
	answer := 0:
	for i from 1 to nops(a) do
		answer := answer + tournament([a[i],newc]):
	od:
	return answer:
fi:
end:

#----------------------------------------------------------------------
# findinlist(a,l)
#----------------------------------------------------------------------
# looks for a in sorted list l
# if it's in the list, returns the [true, #]
# otherwise, returns [false, # where it should be]
findinlist := proc(a,l)
return filsub(a,l,0,nops(l)+1):
end:

#----------------------------------------------------------------------
# filsub(a,l,q,r)
#----------------------------------------------------------------------
# looks for a in list between positions q and r, not inclusive
# if it's in the list, returns the #
# otherwise, returns integer + 1/2 to indicate where to place it
filsub := proc(a,l,q,r)
local s:

# print(`a,l,q,r=`,a,l,q,r):

if r = q+1 then
	return [false, r]:
fi:
s := floor( (q+r)/2):

# print(`s,l[s],a`,s,l[s],a):


if l[s] = a then
	return [true,s]:
elif compare(l[s],a)=1  then
	return filsub(a,l,q,s):
else
	return filsub(a,l,s,r):
fi:
end:

#----------------------------------------------------------------------
# compare(a,b)
#----------------------------------------------------------------------
# this returns -1,0,1 if a<b, a=b, a>b
# where a and b are nested lists of integers
compare := proc(a,b)
local i, j:

if a=b then return 0: fi:
if whattype(a)=integer and whattype(b)=integer then return sign(a-b): fi:
if whattype(a)=integer or whattype(b)=integer then
	print(`Trying to compare integers and lists!  Caution!`):
	if whattype(a)=integer then return -1: fi:
	if whattype(b)=integer then return 1: fi:
fi:
if nops(a) <> nops(b) then 
	return sign(nops(a)-nops(b)):
fi:
for i from 1 to nops(a) do
	j := compare(a[i],b[i]):
	if j<> 0 then return j: fi:
od:
end:

#----------------------------------------------------------------------
# sortorder(alpha,beta)
# sorts by beta, carrying alpha along for the ride
# output is [newalpha,newbeta]
# This will typically be used at the start of a calculation.
sortorder := proc(alpha,beta)
local i, j, newalpha,newbeta, temp:
newalpha := alpha: newbeta:=beta:
for i from 1 to nops(alpha)-1 do
	for j from i+1 to nops(alpha) do



		if compare(beta[i],beta[j])=1 then
			temp := newbeta[i]: newbeta[i] := newbeta[j]: newbeta[j]:=temp:
			temp := newalpha[i]: newalpha[i] := newalpha[j]:newalpha[j]:=temp:
		fi:
	od:
od:
return [newalpha, newbeta]:
end:

#----------------------------------------------------------------------
# trim(alpha,beta)
#----------------------------------------------------------------------
# removes all entries where alpha[i]=0
# returns [newalpha,newbeta]
trim := proc(alpha,beta)
local newalpha,newbeta,i, j:
newalpha := alpha: newbeta := beta:
j := nops(alpha):
for i from j by -1 to 1 do
	if alpha[i] = 0 then
		newalpha := [ op(newalpha[1..(i-1)]),op(newalpha[(i+1)..(nops(newalpha))])]:
		newbeta := [ op(newbeta[1..(i-1)]),op(newbeta[(i+1)..(nops(newbeta))])]:
	fi:
od:
return [newalpha,newbeta]:
end:

#----------------------------------------------------------------------
# trim3(alpha,beta, gamma)
#----------------------------------------------------------------------
# removes all entries where alpha[i]=0
# returns [newalpha,newbeta,newgamma]
trim3 := proc(alpha,beta, gamma)
local newalpha,newbeta,i, j,newgamma:
newalpha := alpha: newbeta := beta: newgamma := gamma:
j := nops(alpha):
for i from j by -1 to 1 do
	if alpha[i] = 0 then
		newalpha := [ op(newalpha[1..(i-1)]),op(newalpha[(i+1)..(nops(newalpha))])]:
		newbeta := [op(newbeta[1..(i-1)]),op(newbeta[(i+1)..(nops(newbeta))])]:
		newgamma := [ op(newgamma[1..(i-1)]),op(newgamma[(i+1)..(nops(newgamma))])]:
	fi:
od:
return [newalpha,newbeta,newgamma]:
end:

#----------------------------------------------------------------------
# checktlist()
#----------------------------------------------------------------------
# This program looks for problems (features, not bugs!) in tlist
# (gal and biggalois)
# usage:  checktlist():
checktlist := proc()
local i:
global tlist1, tlist2, tlist3:
for i from 1 to nops(tlist3) do
	if tlist3[i] = 2 then
		print(`2-trans. needed:`,tlist1[i]):
	elif tlist3[i] = 3 then
		print(`possible Mathieu issue:`,tlist[i]):
	fi:
od:
end:


#----------------------------------------------------------------------
# lr(list=[ul,vl,w,vr,ur])
#----------------------------------------------------------------------
# produces a list [true/false, 1 or 2 things]
#  first entry:  are we there yet?
#  if false, next entries = list of [ul,vl,w,vr,ur]
#  if true, next entry is a subset
#
lr := proc(input)
local output, f1, fT, i, ul, vl, w, vr, ur, ul1,vl1,w1,vr1,ur1;

ul := input[1];
vl := input[2];
w := input[3];
vr := input[4];
ur := input[5];


if w <> [] and vr = [] and ur <>[] then


	if op(w) = 0 and ur[1] = 0 then
		vl := [op(vl),0];
		w := [];  # cut later
		ur := ur[2..nops(ur)]:


		output := wempty(ul,vl,vr,ur);
		return output;

	fi;
	if op(w) = 1 and ur[1] = 1 then

		vl := [op(vl),1];
		w := []; # cut later
		ur := ur[2..nops(ur)]:




		output := wempty(ul,vl,vr,ur);
		return output;
	fi;
	if op(w) = 1 and ur[1] = 0 then
		vl := [op(vl),-1];
		w := [];  # cut later
		ur := ur[2..nops(ur)]:


		output := wempty(ul,vl,vr,ur);
		return output;
	fi;
fi;

f1 := infinity;
fT := infinity;
for i from nops(vr) by -1 to 1  do
	if vr[i] = -1 then fT := i; fi;
	if vr[i] = 1 then f1 := i; fi;
od;

if op(w)=0 then


	if vr[1] = 0 then  # case 1
		vl := [op(vl),0];
		w := [0];
		vr := vr[2..nops(vr)];
		output := [false, [ul,vl,w,vr,ur]];
		return output;
	fi:
	if vr[1] = -1 then  # case 2
		vl := [op(vl),0];
		w := [1];
		vr := vr[2..nops(vr)];
		output := [false, [ul,vl,w,vr,ur]];
		return output;
	fi:
	if vr[1] = 1 then  # case 3
		vl := [op(vl),1];
		w := [0];
		vr := vr[2..nops(vr)];
		output := [false, [ul,vl,w,vr,ur]];
		return output;
	fi:
fi:

if op(w)=1 then
	if fT = infinity and ur[1] = 1 then  # case 4
		vl := [op(vl),1,op(vr)];
		w := [];  # cut later
		vr := [];
		ur := ur[2..nops(ur)];


		output := wempty(ul,vl,vr,ur);
		return output;
	fi:
	if fT=infinity and f1 = infinity and ur[1]=0 then # case 5
		vl := [op(vl), -1,op(vr)];
		w := [];  # cut later
		vr := [];
		ur := ur[2..nops(ur)];

		output := wempty(ul,vl,vr,ur);
		return output;
	fi:

	if fT<f1 then # case 6 (implicitly fT<infinity)
		vl := [op(vl),-1,op(vr[1..(fT-1)])];
#		w := [1]; 
		vr := vr[(fT+1)..nops(vr)];
#		ur := ur[2..nops(ur)];
		output := [false,[ul,vl,w,vr,ur]];
		return output;
	fi:

	if f1=1 and fT<infinity then # case 7a
		vl := [op(vl),1];
#		w := [1]; 
		vr := vr[2..nops(vr)];
#		ur := ur[2..nops(ur)];
		output := [false,[ul,vl,w,vr,ur]];
		return output;
	fi:
	if f1=1 and fT=infinity and ur[1]=0 then # case 7b
		vl := [op(vl),1];
#		w := [1]; 
		vr := vr[2..nops(vr)];
#		ur := ur[2..nops(ur)];
		output := [false,[ul,vl,w,vr,ur]];
		return output;
	fi:



	if f1=1 and fT=infinity and ur[1]=0 then # case 7b
		vl := [op(vl),1];
#		w := [1]; 
		vr := vr[2..nops(vr)];
#		ur := ur[2..nops(ur)];
		output := [false,[ul,vl,w,vr,ur]];
		return output;
	fi:

	if 1<f1 and f1<fT and fT<infinity then # case 8a
		ul1 := ul;
		vl1 := [op(vl),-1];
		w1 := [0];
		vr1 := vr[2..nops(vr)];
		ur1 := ur;
		
		vl := [op(vl),1,0,op(vr[2..(f1-1)])];
		w := [1];
		vr := vr[(f1+1)..nops(vr)];
		output := [false,[ul1,vl1,w1,vr1,ur1],[ul,vl,w,vr,ur]];
		return output;
	fi:

	if 1<f1 and f1<fT and fT=infinity and ur[1]=0 then # case 8b
			# Note that this case is identical
			# to the previous one!
		ul1 := ul;
		vl1 := [op(vl),-1];
		w1 := [0];
		vr1 := vr[2..nops(vr)];
		ur1 := ur;
		
		vl := [op(vl),1,0,op(vr[2..(f1-1)])];
		w := [1];
		vr := vr[(f1+1)..nops(vr)];
		output := [false,[ul1,vl1,w1,vr1,ur1],[ul,vl,w,vr,ur]];
		return output;
	fi:

	print(`PROBLEM!  THIS CASE FELL THROUGH THE CRACKS!`);



fi:


end:
