#OK to post homework

#Zidong Zhang, 24/1/2021, Assignment 1



with(combinat):


#FP(pi) , that inputs a permutation pi of {1,...n}, (where n=nops(pi)) and outputs an integer #between 0 and n indicating the number of fixed points, i.e. the number of places i such that #pi[i]=i. Note that if pi is a derangement, then FP(pi)=0.

FP := proc(pi) local l,n,i:
n:=nops(pi):
l := []:
for i from 1 to n do
	if pi[i] = i then
		l := [op(l),pi[i]]:
	fi:
od:
nops(l):
end:


#RandAlmostDer(n,k) that inputs a positive integer n, another positive integer k, (between 0 #and n) and outputs a random permutation whose number of fixed points is <=k. Note that #RandAlmostder(n,0) does the same thing as RandDer(n)

RandAlmostDer := proc(n,k) local newpi, pi:

pi := randperm(n):
newpi:=pi[1..k]:


while FP(newpi) !=0 do 
	pi:=randperm(n):

od:
pi:	
end:

#EstimateAveFP(n,K) that inputs a positive integer n, and a large positive integer K, picks #random permutations of {1,...,n} (using randperm(n) in the combinat package of Maple) #and outputs the average number of fixed points. By taking K fairly large, say K=10000, and #various n, can you guess the exact value?

EstimateAveFP := proc(n,K) local  pi,m,i:

pi:= randperm(n):
m:=0:

for i from 1 to K do
	m:= m + FP(pi):
	pi:=randperm(n):
od:

m/K:

end:

#EstimateAveFP(5,100000)
#                             9979 
#                             -----
#                           10000
#EstimateAveFP(10,100000)
#                             19993
#                             -----
#                             20000

#EstimateAveFP(50,100000)
#                           99493 
#                           ------
#                           100000

#EstimateAveFP(100,100000)
#                           49921
#                           -----
#                           50000

#EstimateAveFP(500,100000)
#                            4979
#                            ----
#                            5000

# when n = 1000 , it took so long time and my computer's fan is "scearming", so I will stop #at n = 500. the average number is 1 based on my experiment.

##########################6(optional challenge)##############
#Since there is only 1/n of the probability for each pi[i]=i, and there is n of pi[i], so the #average number should goes to 1 I think?





#EstimateMomFP(n,r,K)£¬that inputs a positive integer n, another positive integer, r, and a #large positive integer K, picks randompermutations of {1,...,n} (using randperm(n) in the #combinat package of Maple) and outputs the average number of FP(pi)^r. Note that #EstimateMomFP(n,1,K) does the same thing as EstimateAve(n,K). By taking K fairly large, #say K=10000, r=2, and various n, can you guess the exact value of the average of FF#(pi)^2?


EstimateMomFP:=proc(n,r,K) local   pi,m,i:

pi:= randperm(n):
m:=0:

for i from 1 to K do
	m:= m + FP(pi)**r:
	pi:=randperm(n):
od:

m/K:

end:


#EstimateMomFP(5,1,100000)
#                            6249
#                            ----
#                            6250
### this is the condition when EstimateAveFP(5,100000)


#EstimateMomFP(5,2,100000)
#                            6249
#                            ----
#                            3125

#EstimateMomFP(10,2,100000)
#                           199009
#                           ------
#                           100000

#EstimateMomFP(50,2,100000)
#                           100267
#                           ------
#                           50000 

#EstimateMomFP(100,2,100000)
#                           99643
#                           -----
#                           50000

#as the experiment goes on, the numver should be 2.


##########################8[Bigger Optional challenge] ############

# The exact value of the averrage of random variable FP(pi)^2 should be 2 over all #permutations of length n.
# However, I don't know the answer