# This code can be used to calculate the values of C_Fi for each graph in F_4, and then verify that equation (8) is true for all 4 <= r <= 1000

var('k')

# Each entry in this list corresponds to a different graph in F_4. The enumeration of these graphs matches Figure 3 at the beginning of the proof of Theorem 1.3.
C_F = [0]*11

# This list contains the number of P3 subgraphs in each graph in F_4. This corresponds to equation (4).
P3count = [0]*11
P3count[6] = 1
P3count[7] = 2
P3count[8] = 4
P3count[9] = 6
P3count[10] = 12


# These are the functions for P_0, obtained using Corollary 1.5 of Zykov's Theorem.
P_0 = [(k^3 - 6*k^2 + 11*k - 6)/k^3]*11
P_0[10] = (-6*k^2 + 11*k - 6)/k^3



# These are the coefficient values of P_1, P_2, and P_3.
P_1 = [0]*34
P_2 = [0]*34
P_3 = [0]*34

P_1[0] = 6*k^2 - 12*k + 6
P_1[1] = k^2 - 2*k + 1 
P_1[2] = 1 - k
P_1[3] = 3-3*k
P_1[8] = 2
P_1[9] = 1

P_2[3] = 3
P_2[7] = 1
P_2[6] = -1
P_2[8] = -4

P_3[3] = 3*k^2 - 12*k+12
P_3[7] = k^2 - 8*k + 12
P_3[6] = k^2 - 6*k + 12
P_3[8] = 4*k^2 - 16*k + 16
P_3[9] = 20-8*k
P_3[10] = 24


# Values of the polynomials p_0,p_1,p_2,p_3.
p_0 = 18*(k^2 - 2*k + 1)/(3*k^2 - 11*k + 9)
p_1 = (3*k^3 - 10*k^2 + 7*k)/( (3*k^2 - 11*k + 9)*k^3 )
p_2 = 1/4*(9*k^5 - 32*k^4 + 25*k^3)/( (3*k^2 - 11*k + 9)*k^3) 
p_3 = 1/4*(15*k^3 - 24*k^2 + 7*k)/( (3*k^2 - 11*k + 9)*k^3 ) 


#This combines the terms from equation (6) and calculates the values of C_Fi.
for i in [0..10]:
	C_F[i] = factor(expand(p_0*P_0[i] + p_1*P_1[i] + p_2*P_2[i] + p_3*P_3[i] + P3count[i]))
	print('the coefficient of','F_',i,'equals',C_F[i])

# Now we will verify that coefficients 0,3,8,9,10 are maximum. In order to do so, we will check each C_Fi against C_F0.

#This function will verify that an input function x, which will be C_F[0] - C_F[i], is nonnegative for all r = 4 up to 1000. 

def test_positive(x):
	for r in [4..1000]:
		if x.substitute(k=r) < 0:
			return('C_F_0 is not maximum, coeffcient','fails when','r = ',r)
	return('for all r from 4 to',1000,'opt is max')

this verifies equation (8) for all r from 4 to 1000.
for i in [1..10]:
	print('Testing C_F',i,':',test_positive(C_F[0] - C_F[i]))