with(combinat):

#  FIXED POINTS


# check for fixed points exhaustively
fixed_points := proc(n):
	local cs,bb,ss,c,s,b,foo,ok,bb2,ss2,b2,s2:
	cs := cases(n):
	
	bb := []:
	ss := []:
	bb2 := {}:
	ss2 := {}:

	
	for c in cs do
		if isBBfp(n,c) then
			foo := BB(n,c):
			ok := true:
			for b2 in bb2 do
				if foo = b2 then
					ok := false:
					break:
				fi:
			od:
			if ok then
				bb := [op(bb),foo]:
				bb2 := bb2 union {MultiSet(op(foo))}:
			fi:
		fi:
		if isSSfp(n,c) then
			foo := SS(n,c):
			ok := true:
			for s2 in ss2 do
				if foo = s2 then
					ok := false:
					break:
				fi:
			od:
			if ok then
				ss := [op(ss),foo]:
				ss2 := ss2 union {MultiSet(op(foo))}:
			fi:
		fi:
	
	od:
	printf("Shapely-Shubik fixed points:"):
	for s in ss do
		print(s):
	od:
	printf("Banzhaf fixed points:"):
	for b in bb do
		print(b):
	od:

end:


f := proc(s, q)
	local qi:
	for qi in q do
		if qi minus s = {} then
			return (1):
		fi:
	od:
	return(0):
end:


SS := proc(n, q)
	local PS:
	PS:=powerset(n) minus {{}}:
	[seq(add((nops(S)-1)!*(n-nops(S))!*(f(S,q)-f(S minus {i},q)),S in PS)/n!,i=1..n)]:

end:


BB := proc(n,q)
	local scores,c,PS,S,i:
	scores := [0$n]:
	c:=0:
	PS:=powerset(n) minus {{}}:
	for S in PS do
		if f(S,q) = 0 then
			next:
		fi:
		for i from 1 to n do
			if f(S minus {i},q) = 0 then
				scores[i] := scores[i]+1:
				c:=c+1:
			fi:
			
		od:
	od:
	[seq(scores[i]/c,i=1..n)]:
end:


isSSfp :=proc(n,q) 
	local ss,PS,S,sum,s:
	ss := SS(n,q):
	PS:=powerset(n) minus {{}}:
	
	for S in PS do
		sum := 0:
		for s in S do
			sum := sum + ss[s]:
		od:
		if (f(S,q) = 0 and sum > 1/2) or (f(S,q) = 1 and sum <= 1/2) then
			return(false):
		fi:
		
	od:
	true:
end:

isBBfp :=proc(n,q) 
	local bb,PS,S,sum,s:
	bb := BB(n,q):
	PS:=powerset(n) minus {{}}:
	
	for S in PS do
		sum := 0:
		for s in S do
			sum := sum + bb[s]:
		od:
		if (f(S,q) = 0 and sum > 1/2) or (f(S,q) = 1 and sum <= 1/2) then
			return(false):
		fi:
		
	od:
	true:
end:

apply_perm := proc(q, pi)
	local ans,s,news,i:
	ans := []:
	for s in q do
		news := {}:
		for i in s do
			news := news union {pi[i]}:
		od:
		ans := [op(ans),news]:
	od:

end:

cleanse_symmetry := proc(L,n)

	local ans,bad,ps,q,pi,newq:

	ans := []:
	bad := {}:
	ps := permute(n):
	for q in L do
		if q in bad then
			next:
		fi:
		ans := [op(ans),q]:
		for pi in ps do
			newq := apply_perm(q,pi):
			
			bad := bad union {newq}:
			
		od:
	
	od:
	ans:
end:

cases := proc(n)

	local q1,ans,q,ok1,ok2,i,s,q2,v,ok,c,j,k,qj:

	q1 := iterate({},[op(powerset(n) minus {{}})]):
	
	# n players, no two disjoint coalitions meets the threshold
	
	
	
	q2 := []:
	for q in q1 do
		
		ok1 := true:
		
		for i from 1 to n do
			ok2 := false:
			for s in q do
				if i in s then
					ok2 := true:
					break:
				fi:
			od:
			if not(ok2) then
				ok1 := false:
				break:
			fi:
			
		od:
		if not(ok1) then 
			next:
		fi:
		q2 := [op(q2),q]:
	od:
	
	# each player is relevant
	
	

	
	for i from nops(q2) to 1 by -1 do
	
		v := nops(q2[i][1]):
		if  v > (n/2) + 1 then
			q2 := subsop( i = NULL, q2):
		fi:
	od:
	

	# not too many people are required
	
	
	for i from nops(q2) to 1 by -1 do
	
		ok := true:
		c:=0:
		for j from 1 to n do
			qj := {seq(i,i=1..n)} minus {j}:
			if f(qj, q2[i] ) = 0 then
				if n > 2 and c = 1 then
					ok := false:
					break:
				else 
					c := 1:
				fi:
				for k in qj do
					if not( {j,k} in q2[i] ) then
						ok := false:
						break:
					fi:
				od:
			fi:
			if not(ok) then
				break:
			fi:
		od:
	
		if not(ok) then
			q2 := subsop( i = NULL, q2):
		fi:
	od:
	

	return(q2):
	#q2:=cleanse_symmetry(q2,n):

	# each case is now distinct up to permutations of n


end:




iterate := proc(wc, opts)

	#print(wc,opts):

	local L,i,j,new,newopts:

	L := [wc]:
	
	if nops(wc) = 0 then
		L := []:
	fi:
	
	if nops(opts) = 0 then
		return( L ):
	fi:

	for i from 1 to nops(opts) do
	
		new := opts[i]:
		
		newopts := []:
		for j from i+1 to nops(opts) do
			if new minus opts[j] = {} then
				next:
			fi:
			if new intersect opts[j] = {} then 
				next:
			fi:
			newopts := [op(newopts),opts[j]]:
		
		od:
		L := [op(L),op(iterate([op(wc),new],newopts))]:
	
	od:

	L:


end:






# DIVISOR GAMES






with(combinat):
with(NumberTheory):
EqualEntriesNew := proc(L, Q) local i, S; S := {}; for i to nops(L) do if L[i] <> Q[i] then S := S union {i}; end if; end do; if S = {} then true, S; else false, S; end if; end proc:
SSVc := proc(L, q) local n, z, x, f, J, V, i1, k; n := nops(L); V := []; for J to n do f := mul(1 + z*x^L[i1], i1 = 1 .. J - 1)*mul(1 + z*x^L[i1], i1 = J + 1 .. n); f := add(coeff(f, x, i1), i1 = q - L[J] .. q - 1); V := [op(V), add(coeff(f, z, k - 1)*(k - 1)!*(n - k)!, k = 1 .. n)]; end do; V/n!; end proc:
BBVc := proc(L, q) local n, z, x, f, J, V, i1, k, i; n := nops(L); V := []; for J to n do f := mul(1 + x^L[i1], i1 = 1 .. J - 1)*mul(1 + x^L[i1], i1 = J + 1 .. n); V := [op(V), add(coeff(f, x, i1), i1 = q - L[J] .. q - 1)]; end do; V/add(V[i], i = 1 .. nops(V)); end proc:
QueryInteger := proc(n) EqualEntriesNew(BBVc(Divisors(n), floor(1/2*SumOfDivisors(n) + 1/2)), SSVc(Divisors(n), floor(1/2*SumOfDivisors(n) + 1/2))), evalf(SumOfDivisors(n)/n); end proc:
f1:=proc(): #checks for integers n which SS and BB values differ, but for not al
local i;
for i from 4000 to 10000 do
    if SumOfDivisors(i)/i < 2.18 then if nops(QueryInteger(i)[2]) <> nops(Divisors(i)) and nops(QueryInteger(i)[2]) <> 0 then print(i, QueryInteger(i)); end if; end if;
end do; end proc:
f2:=proc(): # computes QueryInteger(n) for n = 2,719
local n;
for n from 2 to 719 do
    print(n, QueryInteger(n));
end do; end proc: