#Ok to post homework
#Tifany Tong, November 22nd, 2020, HW #21

# Question 1.

# L:=WtEdConGclever(50,a):
# coeff(L[50],a,70) = 557570125262936299552095051452832858951958306952392379552719516748979534894703561052813706583449695070191616000000

# -----------------------------------------------------------------
# Question 2.

helper:=proc(N,a,k) local f,i,x:
option remember:
f:=(log(add((1+a)^binomial(i,2)*x^i/i!,i=0..N+1)))^k/k!:
f:=taylor(f,x=0,N+1):
[seq(expand(coeff(f,x,i)*i!),i=1..N)]:
end:
           
# L:=helper(50,a,2):
# coeff(L[50],a,70) =
# 2267762282072205271778791270680596849921516338574571823883456519281730035101965414722739023408444177334403072000000

# L := helper(50, a, 3);
# coeff(L[50], a, 70) =
# 4084614615513849564059280183743190588234377301848961646892154560205839868330317696952259164363869748482983526400000

# L := helper(50, a, 4);
# coeff(L[50], a, 70) =
# 4336520950646656327224993268880264842114984745013581637572959168612209444462829083123209142308423023328834795520000

# L := helper(50, a, 5);
# coeff(L[50], a, 70) = 
# 3046643850794698885807347271650949921058341236747610399052555568087884858386903803880100124793197196674954096998400
# -----------------------------------------------------------------

# Question 3.

# My EstimateCutOff(30, 1000) got me 79 while FindCutOff(30, 1.05) = 84

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

# Question 4.

# evalf(FindCutOff(20, 1.05)/20) = 2.450000000

# evalf(FindCutOff(30, 1.05)/30) = 2.800000000

# evalf(FindCutOff(40, 1.05)/40) = 3.

# evalf(FindCutOff(50, 1.05)/50) = 3.160000000

# evalf(FindCutOff(60, 1.05)/60) = 3.300000000

# The general trend seems to be that the ratio between the cutoff for the number of edges when there is
# an average of 1.05 connected components and the number of vertices increases as the number of vertices
# goes up. 

# ------------------------------------------------------------------
# Question 5.

EstProbKcomps := proc(n, m, k, N) local i, co, G:
co := 0:
for i to N do 
    G := RandGr(n, m):
    if nops(CCs(G)) = k then 
        co := co + 1:
    fi:
od:
evalf(co/N):
end:


# ------------------------------------------------------------------
# Question 6.

WtEdConGcleverHelper := proc(n, a, k) local f, i, x:
option remember:
f := log(sum((a + 1)^(1/2*p*(p - 1))*x^p/p!, p = 0 .. infinity))^k/k!:
f := taylor(f, x = 0, n + 1):
[seq(expand(coeff(f, x, i)*i!), i = 1 .. n)]:
end:


ProbKcomponents := proc(n, m, k) local a:
evalf(coeff(WtEdConGcleverHelper(n, a, k)[n], a, m)/binomial(1/2*n*(n - 1), m)):
end:

# ProbKcomponents(40,50,2) and EstProbKcomps(40,50,2,300) are relatively close. ProbKcomponents(40,50,2) = 0.1033333333333333333333333333333333333333
# and ProbKcomponents(40,50,2) = 0.09801230117869971836432456621233340685540