Kemeny constant for parallel paths between cliques

1784 days ago by butler

def B(k,A,B,C): G = Graph() for i in range(A): new_verts=range(i*k,i*k+k) #G.add_vertices(new_verts) if (i == 0) or (i == A-1): for p in range(k): for q in range(p+1,k): G.add_edge([new_verts[p],new_verts[q]]) if i: for p in new_verts: G.add_edge([p,p-k]) new_verts = range(k*A,k*A+B) for p in range(B): for q in range(k): G.add_edge([q,new_verts[p]]) for q in range(p+1,B): G.add_edge([new_verts[p],new_verts[q]]) new_verts = range(k*A+B,k*A+B+C) for p in range(C): for q in range(k): G.add_edge([(A-1)*k+q,new_verts[p]]) for q in range(p+1,C): G.add_edge([new_verts[p],new_verts[q]]) return G def D(a,b,c,d): G = Graph() G.add_edge([0,1]) G.add_edge([0,a+b]) G.add_edge([a+b+c+d+1,a]) G.add_edge([a+b+c+d+1,a+1]) for i in range(1,a): G.add_edge([i,i+1]) for i in range(a+1,a+b): G.add_edge([i,i+1]) for i in range(a+b+1,a+b+c+1): G.add_edge([i,0]) for j in range(i+1,a+b+c+1): G.add_edge([i,j]) for i in range(a+b+c+1,a+b+c+d+1): G.add_edge([i,a+b+c+d+1]) for j in range(i+1,a+b+c+d+1): G.add_edge([i,j]) return G def Kemeny(G): L=Matrix(QQ,G.order(),G.order()) for i in range(G.order()): for j in G.neighbors(i): L[i,j] = -1/G.degree(i) L[i,i] = 1 f=L.characteristic_polynomial(x) return -1*f[2]/f[1] def educated_guess(k,a,b,c): return -(-a^3*b^4*c^2*k^2 - a^3*b^2*c^4*k^2 - 2*a^3*b^4*c*k^3 - 2*a^4*b^2*c^2*k^3 - 4*a^3*b^3*c^2*k^3 - 4*a^3*b^2*c^3*k^3 - 2*a^3*b*c^4*k^3 - a^3*b^4*k^4 - 4*a^4*b^2*c*k^4 - 8*a^3*b^3*c*k^4 - 4*a^4*b*c^2*k^4 - 12*a^3*b^2*c^2*k^4 - 8*a^3*b*c^3*k^4 - a^3*c^4*k^4 - 2*a^4*b^2*k^5 - 4*a^3*b^3*k^5 - 8*a^4*b*c*k^5 - 16*a^3*b^2*c*k^5 - 2*a^4*c^2*k^5 - 16*a^3*b*c^2*k^5 - 4*a^3*c^3*k^5 - 4*a^4*b*k^6 - 7*a^3*b^2*k^6 - 4*a^4*c*k^6 - 16*a^3*b*c*k^6 - 7*a^3*c^2*k^6 - 2*a^4*k^7 - 6*a^3*b*k^7 - 6*a^3*c*k^7 - 2*a^3*k^8 - 3*a^2*b^4*c^4 - 2*a^3*b^4*c^2*k - 12*a^2*b^4*c^3*k - 2*a^3*b^2*c^4*k - 12*a^2*b^3*c^4*k - 4*a^3*b^4*c*k^2 - 7*a^3*b^3*c^2*k^2 - 15*a^2*b^4*c^2*k^2 - 7*a^3*b^2*c^3*k^2 - 48*a^2*b^3*c^3*k^2 - 4*a^3*b*c^4*k^2 - 15*a^2*b^2*c^4*k^2 - 2*a^3*b^4*k^3 - 14*a^3*b^3*c*k^3 - 6*a^2*b^4*c*k^3 - 10*a^3*b^2*c^2*k^3 - 60*a^2*b^3*c^2*k^3 - 14*a^3*b*c^3*k^3 - 60*a^2*b^2*c^3*k^3 - 2*a^3*c^4*k^3 - 6*a^2*b*c^4*k^3 - 7*a^3*b^3*k^4 - 7*a^3*b^2*c*k^4 - 24*a^2*b^3*c*k^4 - 7*a^3*b*c^2*k^4 - 72*a^2*b^2*c^2*k^4 - 7*a^3*c^3*k^4 - 24*a^2*b*c^3*k^4 - 2*a^3*b^2*k^5 + 12*a^3*b*c*k^5 - 24*a^2*b^2*c*k^5 - 2*a^3*c^2*k^5 - 24*a^2*b*c^2*k^5 + 9*a^3*b*k^6 + 3*a^2*b^2*k^6 + 9*a^3*c*k^6 + 3*a^2*c^2*k^6 + 6*a^3*k^7 + 6*a^2*b*k^7 + 6*a^2*c*k^7 + 3*a^2*k^8 + 3*a^2*b^4*c^3 + 3*a^2*b^3*c^4 + 6*a*b^4*c^4 + 2*a^3*b^3*c^2*k + 15*a^2*b^4*c^2*k - 3*a*b^5*c^2*k + 2*a^3*b^2*c^3*k + 24*a^2*b^3*c^3*k + 21*a*b^4*c^3*k + 15*a^2*b^2*c^4*k + 21*a*b^3*c^4*k - 3*a*b^2*c^5*k + 4*a^3*b^3*c*k^2 + 18*a^2*b^4*c*k^2 - 6*a*b^5*c*k^2 + 12*a^3*b^2*c^2*k^2 + 66*a^2*b^3*c^2*k^2 + 10*a*b^4*c^2*k^2 + 4*a^3*b*c^3*k^2 + 66*a^2*b^2*c^3*k^2 + 72*a*b^3*c^3*k^2 + 18*a^2*b*c^4*k^2 + 10*a*b^2*c^4*k^2 - 6*a*b*c^5*k^2 + 2*a^3*b^3*k^3 + 6*a^2*b^4*k^3 - 3*a*b^5*k^3 + 12*a^3*b^2*c*k^3 + 60*a^2*b^3*c*k^3 - 19*a*b^4*c*k^3 + 12*a^3*b*c^2*k^3 + 98*a^2*b^2*c^2*k^3 + 52*a*b^3*c^2*k^3 + 2*a^3*c^3*k^3 + 60*a^2*b*c^3*k^3 + 52*a*b^2*c^3*k^3 + 6*a^2*c^4*k^3 - 19*a*b*c^4*k^3 - 3*a*c^5*k^3 + 2*a^3*b^2*k^4 + 15*a^2*b^3*k^4 - 14*a*b^4*k^4 + 40*a^2*b^2*c*k^4 - 28*a*b^3*c*k^4 + 2*a^3*c^2*k^4 + 40*a^2*b*c^2*k^4 + 24*a*b^2*c^2*k^4 + 15*a^2*c^3*k^4 - 28*a*b*c^3*k^4 - 14*a*c^4*k^4 - 8*a^3*b*k^5 - 7*a^2*b^2*k^5 - 29*a*b^3*k^5 - 8*a^3*c*k^5 - 40*a^2*b*c*k^5 - 53*a*b^2*c*k^5 - 7*a^2*c^2*k^5 - 53*a*b*c^2*k^5 - 29*a*c^3*k^5 - 8*a^3*k^6 - 38*a^2*b*k^6 - 38*a*b^2*k^6 - 38*a^2*c*k^6 - 68*a*b*c*k^6 - 38*a*c^2*k^6 - 22*a^2*k^7 - 30*a*b*k^7 - 30*a*c*k^7 - 10*a*k^8 - 3*a^2*b^3*c^3 - 12*a*b^4*c^3 - 12*a*b^3*c^4 - 3*b^4*c^4 - 21*a^2*b^3*c^2*k - 25*a*b^4*c^2*k + 3*b^5*c^2*k - 21*a^2*b^2*c^3*k - 90*a*b^3*c^3*k - 9*b^4*c^3*k - 25*a*b^2*c^4*k - 9*b^3*c^4*k + 3*b^2*c^5*k + 4*a^3*b^2*c*k^2 - 30*a^2*b^3*c*k^2 - 8*a*b^4*c*k^2 + 6*b^5*c*k^2 + 4*a^3*b*c^2*k^2 - 60*a^2*b^2*c^2*k^2 - 134*a*b^3*c^2*k^2 + 6*b^4*c^2*k^2 - 30*a^2*b*c^3*k^2 - 134*a*b^2*c^3*k^2 - 24*b^3*c^3*k^2 - 8*a*b*c^4*k^2 + 6*b^2*c^4*k^2 + 6*b*c^5*k^2 + 4*a^3*b^2*k^3 - 12*a^2*b^3*k^3 + 5*a*b^4*k^3 + 3*b^5*k^3 + 16*a^3*b*c*k^3 - 30*a^2*b^2*c*k^3 - 22*a*b^3*c*k^3 + 27*b^4*c*k^3 + 4*a^3*c^2*k^3 - 30*a^2*b*c^2*k^3 - 122*a*b^2*c^2*k^3 + 12*b^3*c^2*k^3 - 12*a^2*c^3*k^3 - 22*a*b*c^3*k^3 + 12*b^2*c^3*k^3 + 5*a*c^4*k^3 + 27*b*c^4*k^3 + 3*c^5*k^3 + 12*a^3*b*k^4 + 9*a^2*b^2*k^4 + 34*a*b^3*k^4 + 15*b^4*k^4 + 12*a^3*c*k^4 + 75*a^2*b*c*k^4 + 94*a*b^2*c*k^4 + 60*b^3*c*k^4 + 9*a^2*c^2*k^4 + 94*a*b*c^2*k^4 + 60*b^2*c^2*k^4 + 34*a*c^3*k^4 + 60*b*c^3*k^4 + 15*c^4*k^4 + 8*a^3*k^5 + 75*a^2*b*k^5 + 107*a*b^2*k^5 + 33*b^3*k^5 + 75*a^2*c*k^5 + 246*a*b*c*k^5 + 93*b^2*c*k^5 + 107*a*c^2*k^5 + 93*b*c^2*k^5 + 33*c^3*k^5 + 54*a^2*k^6 + 138*a*b*k^6 + 42*b^2*k^6 + 138*a*c*k^6 + 84*b*c*k^6 + 42*c^2*k^6 + 60*a*k^7 + 30*b*k^7 + 30*c*k^7 + 9*k^8 + 6*a*b^4*c^2 + 18*a*b^3*c^3 + 9*b^4*c^3 + 6*a*b^2*c^4 + 9*b^3*c^4 + 6*a*b^4*c*k - 3*b^5*c*k + 12*a^2*b^2*c^2*k + 79*a*b^3*c^2*k + 6*b^4*c^2*k + 79*a*b^2*c^3*k + 60*b^3*c^3*k + 6*a*b*c^4*k + 6*b^2*c^4*k - 3*b*c^5*k - 3*b^5*k^2 + 62*a*b^3*c*k^2 - 30*b^4*c*k^2 + 150*a*b^2*c^2*k^2 + 36*b^3*c^2*k^2 + 62*a*b*c^3*k^2 + 36*b^2*c^3*k^2 - 30*b*c^4*k^2 - 3*c^5*k^2 - 12*a^2*b^2*k^3 + a*b^3*k^3 - 27*b^4*k^3 - 60*a^2*b*c*k^3 - 15*a*b^2*c*k^3 - 102*b^3*c*k^3 - 12*a^2*c^2*k^3 - 15*a*b*c^2*k^3 - 72*b^2*c^2*k^3 + a*c^3*k^3 - 102*b*c^3*k^3 - 27*c^4*k^3 - 60*a^2*b*k^4 - 92*a*b^2*k^4 - 87*b^3*k^4 - 60*a^2*c*k^4 - 294*a*b*c*k^4 - 258*b^2*c*k^4 - 92*a*c^2*k^4 - 258*b*c^2*k^4 - 87*c^3*k^4 - 48*a^2*k^5 - 223*a*b*k^5 - 156*b^2*k^5 - 223*a*c*k^5 - 330*b*c*k^5 - 156*c^2*k^5 - 130*a*k^6 - 147*b*k^6 - 147*c*k^6 - 54*k^7 - 6*a*b^3*c^2 - 9*b^4*c^2 - 6*a*b^2*c^3 - 21*b^3*c^3 - 9*b^2*c^4 - 6*a*b^3*c*k + 3*b^4*c*k - 66*a*b^2*c^2*k - 69*b^3*c^2*k - 6*a*b*c^3*k - 69*b^2*c^3*k + 3*b*c^4*k + 9*b^4*k^2 + 12*a^2*b*c*k^2 - 25*a*b^2*c*k^2 + 21*b^3*c*k^2 - 25*a*b*c^2*k^2 - 78*b^2*c^2*k^2 + 21*b*c^3*k^2 + 9*c^4*k^2 + 12*a^2*b*k^3 + 29*a*b^2*k^3 + 57*b^3*k^3 + 12*a^2*c*k^3 + 152*a*b*c*k^3 + 177*b^2*c*k^3 + 29*a*c^2*k^3 + 177*b*c^2*k^3 + 57*c^3*k^3 + 12*a^2*k^4 + 159*a*b*k^4 + 177*b^2*k^4 + 159*a*c*k^4 + 423*b*c*k^4 + 177*c^2*k^4 + 124*a*k^5 + 252*b*k^5 + 252*c*k^5 + 120*k^6 + 3*b^4*c + 15*b^3*c^2 + 15*b^2*c^3 + 3*b*c^4 + 6*a*b^2*c*k + 15*b^3*c*k + 6*a*b*c^2*k + 84*b^2*c^2*k + 15*b*c^3*k - 9*b^3*k^2 - 24*a*b*c*k^2 + 12*b^2*c*k^2 + 12*b*c^2*k^2 - 9*c^3*k^2 - 42*a*b*k^3 - 60*b^2*k^3 - 42*a*c*k^3 - 186*b*c*k^3 - 60*c^2*k^3 - 48*a*k^4 - 177*b*k^4 - 177*c*k^4 - 120*k^5 - 3*b^3*c - 6*b^2*c^2 - 3*b*c^3 - 18*b^2*c*k - 18*b*c^2*k + 3*b^2*k^2 + 6*b*c*k^2 + 3*c^2*k^2 + 42*b*k^3 + 42*c*k^3 + 48*k^4)/(3*a*b^4*c^2*k + 3*a*b^2*c^4*k + 6*a*b^4*c*k^2 + 6*a^2*b^2*c^2*k^2 + 12*a*b^3*c^2*k^2 + 12*a*b^2*c^3*k^2 + 6*a*b*c^4*k^2 + 3*a*b^4*k^3 + 12*a^2*b^2*c*k^3 + 24*a*b^3*c*k^3 + 12*a^2*b*c^2*k^3 + 36*a*b^2*c^2*k^3 + 24*a*b*c^3*k^3 + 3*a*c^4*k^3 + 6*a^2*b^2*k^4 + 12*a*b^3*k^4 + 24*a^2*b*c*k^4 + 48*a*b^2*c*k^4 + 6*a^2*c^2*k^4 + 48*a*b*c^2*k^4 + 12*a*c^3*k^4 + 12*a^2*b*k^5 + 21*a*b^2*k^5 + 12*a^2*c*k^5 + 48*a*b*c*k^5 + 21*a*c^2*k^5 + 6*a^2*k^6 + 18*a*b*k^6 + 18*a*c*k^6 + 6*a*k^7 - 3*a*b^3*c^2*k - 3*b^4*c^2*k - 3*a*b^2*c^3*k - 3*b^2*c^4*k - 6*a*b^3*c*k^2 - 6*b^4*c*k^2 - 30*a*b^2*c^2*k^2 - 12*b^3*c^2*k^2 - 6*a*b*c^3*k^2 - 12*b^2*c^3*k^2 - 6*b*c^4*k^2 - 3*a*b^3*k^3 - 3*b^4*k^3 - 51*a*b^2*c*k^3 - 24*b^3*c*k^3 - 51*a*b*c^2*k^3 - 36*b^2*c^2*k^3 - 3*a*c^3*k^3 - 24*b*c^3*k^3 - 3*c^4*k^3 - 24*a*b^2*k^4 - 12*b^3*k^4 - 84*a*b*c*k^4 - 48*b^2*c*k^4 - 24*a*c^2*k^4 - 48*b*c^2*k^4 - 12*c^3*k^4 - 39*a*b*k^5 - 21*b^2*k^5 - 39*a*c*k^5 - 48*b*c*k^5 - 21*c^2*k^5 - 18*a*k^6 - 18*b*k^6 - 18*c*k^6 - 6*k^7 + 3*b^4*c*k + 6*b^3*c^2*k + 6*b^2*c^3*k + 3*b*c^4*k + 3*b^4*k^2 + 6*a*b^2*c*k^2 + 24*b^3*c*k^2 + 6*a*b*c^2*k^2 + 42*b^2*c^2*k^2 + 24*b*c^3*k^2 + 3*c^4*k^2 + 6*a*b^2*k^3 + 18*b^3*k^3 + 24*a*b*c*k^3 + 84*b^2*c*k^3 + 6*a*c^2*k^3 + 84*b*c^2*k^3 + 18*c^3*k^3 + 18*a*b*k^4 + 48*b^2*k^4 + 18*a*c*k^4 + 120*b*c*k^4 + 48*c^2*k^4 + 12*a*k^5 + 57*b*k^5 + 57*c*k^5 + 24*k^6 - 3*b^3*c*k - 6*b^2*c^2*k - 3*b*c^3*k - 3*b^3*k^2 - 27*b^2*c*k^2 - 27*b*c^2*k^2 - 3*c^3*k^2 - 21*b^2*k^3 - 66*b*c*k^3 - 21*c^2*k^3 - 42*b*k^4 - 42*c*k^4 - 24*k^5) 
       
# the following is used to generate the formula for B(k,a,b,c) a,b,c,k,x=var('a,b,c,k,x') def p(n): return (-1/2)^(n-2-a)*(n+3)*(n+2)*(n+1)*n*(n-1)/120*x^2+(-1/2)^(n-1-a)*(n+2)*(n+1)*n*x/6+(-1/2)^(n-a)*(n+1) d = 0 d += p(a-2)*(x-(b+k+1)/(b+k))*(x-(c+k+1)/(c+k)) d += p(a-4)*(1/(2*(b+k)))*(1/(2*(c+k))) d -= p(a-3)*(x-(b+k+1)/(b+k))*(1/(2*(c+k))) d -= p(a-3)*(x-(c+k+1)/(c+k))*(1/(2*(b+k))) d = d*-1*(b*c + (b + c)*k + k^2)*3/4 d0 = d.coefficients(x)[0][0].canonicalize_radical().expand() d1 = d.coefficients(x)[1][0].canonicalize_radical().expand() e = 0 e += p(a-2)*(x-b*k/(b^2-b+b*k))*(x-(b+1)/(b+k))*(x-(c+1)/(c+k))*(x-c*k/(c^2-c+c*k)) e -= p(a-2)*(b^2*k/((b+k)*(b^2-b+b*k)))*(x-(c+1)/(c+k))*(x-c*k/(c^2-c+c*k)) e -= p(a-2)*(x-b*k/(b^2-b+b*k))*(x-(b+1)/(b+k))*(c^2*k/((c+k)*(c^2-c+c*k))) e += p(a-2)*(b^2*k/((b+k)*(b^2-b+b*k)))*(c^2*k/((c+k)*(c^2-c+c*k))) e -= p(a-3)*(x-b*k/(b^2-b+b*k))*(1/(2*(b+k)))*(x-(c+1)/(c+k))*(x-c*k/(c^2-c+c*k)) e += p(a-3)*(x-b*k/(b^2-b+b*k))*(1/(2*(b+k)))*(c^2*k/((c+k)*(c^2-c+c*k))) e -= p(a-3)*(x-b*k/(b^2-b+b*k))*(x-(b+1)/(b+k))*(1/(2*(c+k)))*(x-c*k/(c^2-c+c*k)) e += p(a-3)*(b^2*k/((b+k)*(b^2-b+b*k)))*(1/(2*(c+k)))*(x-c*k/(c^2-c+c*k)) e += p(a-4)*(x-b*k/(b^2-b+b*k))*(1/(2*(b+k)))*(1/(2*(c+k)))*(x-c*k/(c^2-c+c*k)) e *= (3/4)*(b^2*c^2 + 2*b^2*c*k + 2*b*c^2*k + b^2*k^2 + 4*b*c*k^2 + c^2*k^2 + 2*b*k^3 + 2*c*k^3 + k^4 - b^2*c - b*c^2 - b^2*k - 4*b*c*k - c^2*k - 3*b*k^2 - 3*c*k^2 - 2*k^3 + b*c + b*k + c*k + k^2) e1 = e.coefficients(x)[0][0].canonicalize_radical().expand() e2 = e.coefficients(x)[1][0].canonicalize_radical().expand() f = (-(b-1)*(b+k-1)*x+(b+k))*(-(c-1)*(c+k-1)*x+(c+k))*((k-1)*d1*x+d0)*(e2*x+e1) f *= (1/3) f1 = f.coefficients(x)[0][0].canonicalize_radical().expand() f2 = f.coefficients(x)[1][0].canonicalize_radical().expand() 
       
print "d0: ", d0, "\n" print "d1: ", d1, "\n" print "e1: ", e1, "\n" print "e2: ", e2, "\n" print "f1: ", f1, "\n" print "f2: ", f2, "\n" 
       
d0:  -3*a*b*c - 3*a*b*k - 3*a*c*k - 3*a*k^2 + 3*b*c + 3*b*k + 3*c*k +
3*k^2 - 3*b - 3*c - 6*k 

d1:  a^3*b*c + a^3*b*k + a^3*c*k + a^3*k^2 - 3*a^2*b*c - 3*a^2*b*k -
3*a^2*c*k - 3*a^2*k^2 + 3*a^2*b + 3*a^2*c + 8*a*b*c + 6*a^2*k + 8*a*b*k
+ 8*a*c*k + 8*a*k^2 - 9*a*b - 9*a*c - 6*b*c - 18*a*k - 6*b*k - 6*c*k -
6*k^2 + 6*a + 9*b + 9*c + 18*k - 12 

e1:  -3*b^2*k - 3*c^2*k - 6*a*k^2 - 6*b*k^2 - 6*c*k^2 - 6*k^3 + 3*b*k +
3*c*k + 12*k^2 

e2:  3*a*b^2*c^2 + 3*a^2*b^2*k + 6*a*b^2*c*k + 3*a^2*c^2*k + 6*a*b*c^2*k
+ 2*a^3*k^2 + 6*a^2*b*k^2 + 3*a*b^2*k^2 + 6*a^2*c*k^2 + 12*a*b*c*k^2 +
3*a*c^2*k^2 + 6*a^2*k^3 + 6*a*b*k^3 + 6*a*c*k^3 + 3*a*k^4 - 3*a*b^2*c -
3*a*b*c^2 - 3*b^2*c^2 - 3*a^2*b*k - 9*a*b^2*k - 3*a^2*c*k - 12*a*b*c*k -
6*b^2*c*k - 9*a*c^2*k - 6*b*c^2*k - 12*a^2*k^2 - 21*a*b*k^2 - 3*b^2*k^2
- 21*a*c*k^2 - 12*b*c*k^2 - 3*c^2*k^2 - 18*a*k^3 - 6*b*k^3 - 6*c*k^3 -
3*k^4 + 3*a*b*c + 6*b^2*c + 6*b*c^2 + 15*a*b*k + 12*b^2*k + 15*a*c*k +
24*b*c*k + 12*c^2*k + 34*a*k^2 + 30*b*k^2 + 30*c*k^2 + 24*k^3 - 3*b^2 -
9*b*c - 3*c^2 - 12*a*k - 30*b*k - 30*c*k - 48*k^2 + 3*b + 3*c + 24*k 

f1:  3*a*b^4*c^2*k + 3*a*b^2*c^4*k + 6*a*b^4*c*k^2 + 6*a^2*b^2*c^2*k^2 +
12*a*b^3*c^2*k^2 + 12*a*b^2*c^3*k^2 + 6*a*b*c^4*k^2 + 3*a*b^4*k^3 +
12*a^2*b^2*c*k^3 + 24*a*b^3*c*k^3 + 12*a^2*b*c^2*k^3 + 36*a*b^2*c^2*k^3
+ 24*a*b*c^3*k^3 + 3*a*c^4*k^3 + 6*a^2*b^2*k^4 + 12*a*b^3*k^4 +
24*a^2*b*c*k^4 + 48*a*b^2*c*k^4 + 6*a^2*c^2*k^4 + 48*a*b*c^2*k^4 +
12*a*c^3*k^4 + 12*a^2*b*k^5 + 21*a*b^2*k^5 + 12*a^2*c*k^5 + 48*a*b*c*k^5
+ 21*a*c^2*k^5 + 6*a^2*k^6 + 18*a*b*k^6 + 18*a*c*k^6 + 6*a*k^7 -
3*a*b^3*c^2*k - 3*b^4*c^2*k - 3*a*b^2*c^3*k - 3*b^2*c^4*k -
6*a*b^3*c*k^2 - 6*b^4*c*k^2 - 30*a*b^2*c^2*k^2 - 12*b^3*c^2*k^2 -
6*a*b*c^3*k^2 - 12*b^2*c^3*k^2 - 6*b*c^4*k^2 - 3*a*b^3*k^3 - 3*b^4*k^3 -
51*a*b^2*c*k^3 - 24*b^3*c*k^3 - 51*a*b*c^2*k^3 - 36*b^2*c^2*k^3 -
3*a*c^3*k^3 - 24*b*c^3*k^3 - 3*c^4*k^3 - 24*a*b^2*k^4 - 12*b^3*k^4 -
84*a*b*c*k^4 - 48*b^2*c*k^4 - 24*a*c^2*k^4 - 48*b*c^2*k^4 - 12*c^3*k^4 -
39*a*b*k^5 - 21*b^2*k^5 - 39*a*c*k^5 - 48*b*c*k^5 - 21*c^2*k^5 -
18*a*k^6 - 18*b*k^6 - 18*c*k^6 - 6*k^7 + 3*b^4*c*k + 6*b^3*c^2*k +
6*b^2*c^3*k + 3*b*c^4*k + 3*b^4*k^2 + 6*a*b^2*c*k^2 + 24*b^3*c*k^2 +
6*a*b*c^2*k^2 + 42*b^2*c^2*k^2 + 24*b*c^3*k^2 + 3*c^4*k^2 + 6*a*b^2*k^3
+ 18*b^3*k^3 + 24*a*b*c*k^3 + 84*b^2*c*k^3 + 6*a*c^2*k^3 + 84*b*c^2*k^3
+ 18*c^3*k^3 + 18*a*b*k^4 + 48*b^2*k^4 + 18*a*c*k^4 + 120*b*c*k^4 +
48*c^2*k^4 + 12*a*k^5 + 57*b*k^5 + 57*c*k^5 + 24*k^6 - 3*b^3*c*k -
6*b^2*c^2*k - 3*b*c^3*k - 3*b^3*k^2 - 27*b^2*c*k^2 - 27*b*c^2*k^2 -
3*c^3*k^2 - 21*b^2*k^3 - 66*b*c*k^3 - 21*c^2*k^3 - 42*b*k^4 - 42*c*k^4 -
24*k^5 

f2:  -a^3*b^4*c^2*k^2 - a^3*b^2*c^4*k^2 - 2*a^3*b^4*c*k^3 -
2*a^4*b^2*c^2*k^3 - 4*a^3*b^3*c^2*k^3 - 4*a^3*b^2*c^3*k^3 -
2*a^3*b*c^4*k^3 - a^3*b^4*k^4 - 4*a^4*b^2*c*k^4 - 8*a^3*b^3*c*k^4 -
4*a^4*b*c^2*k^4 - 12*a^3*b^2*c^2*k^4 - 8*a^3*b*c^3*k^4 - a^3*c^4*k^4 -
2*a^4*b^2*k^5 - 4*a^3*b^3*k^5 - 8*a^4*b*c*k^5 - 16*a^3*b^2*c*k^5 -
2*a^4*c^2*k^5 - 16*a^3*b*c^2*k^5 - 4*a^3*c^3*k^5 - 4*a^4*b*k^6 -
7*a^3*b^2*k^6 - 4*a^4*c*k^6 - 16*a^3*b*c*k^6 - 7*a^3*c^2*k^6 - 2*a^4*k^7
- 6*a^3*b*k^7 - 6*a^3*c*k^7 - 2*a^3*k^8 - 3*a^2*b^4*c^4 -
2*a^3*b^4*c^2*k - 12*a^2*b^4*c^3*k - 2*a^3*b^2*c^4*k - 12*a^2*b^3*c^4*k
- 4*a^3*b^4*c*k^2 - 7*a^3*b^3*c^2*k^2 - 15*a^2*b^4*c^2*k^2 -
7*a^3*b^2*c^3*k^2 - 48*a^2*b^3*c^3*k^2 - 4*a^3*b*c^4*k^2 -
15*a^2*b^2*c^4*k^2 - 2*a^3*b^4*k^3 - 14*a^3*b^3*c*k^3 - 6*a^2*b^4*c*k^3
- 10*a^3*b^2*c^2*k^3 - 60*a^2*b^3*c^2*k^3 - 14*a^3*b*c^3*k^3 -
60*a^2*b^2*c^3*k^3 - 2*a^3*c^4*k^3 - 6*a^2*b*c^4*k^3 - 7*a^3*b^3*k^4 -
7*a^3*b^2*c*k^4 - 24*a^2*b^3*c*k^4 - 7*a^3*b*c^2*k^4 -
72*a^2*b^2*c^2*k^4 - 7*a^3*c^3*k^4 - 24*a^2*b*c^3*k^4 - 2*a^3*b^2*k^5 +
12*a^3*b*c*k^5 - 24*a^2*b^2*c*k^5 - 2*a^3*c^2*k^5 - 24*a^2*b*c^2*k^5 +
9*a^3*b*k^6 + 3*a^2*b^2*k^6 + 9*a^3*c*k^6 + 3*a^2*c^2*k^6 + 6*a^3*k^7 +
6*a^2*b*k^7 + 6*a^2*c*k^7 + 3*a^2*k^8 + 3*a^2*b^4*c^3 + 3*a^2*b^3*c^4 +
6*a*b^4*c^4 + 2*a^3*b^3*c^2*k + 15*a^2*b^4*c^2*k - 3*a*b^5*c^2*k +
2*a^3*b^2*c^3*k + 24*a^2*b^3*c^3*k + 21*a*b^4*c^3*k + 15*a^2*b^2*c^4*k +
21*a*b^3*c^4*k - 3*a*b^2*c^5*k + 4*a^3*b^3*c*k^2 + 18*a^2*b^4*c*k^2 -
6*a*b^5*c*k^2 + 12*a^3*b^2*c^2*k^2 + 66*a^2*b^3*c^2*k^2 +
10*a*b^4*c^2*k^2 + 4*a^3*b*c^3*k^2 + 66*a^2*b^2*c^3*k^2 +
72*a*b^3*c^3*k^2 + 18*a^2*b*c^4*k^2 + 10*a*b^2*c^4*k^2 - 6*a*b*c^5*k^2 +
2*a^3*b^3*k^3 + 6*a^2*b^4*k^3 - 3*a*b^5*k^3 + 12*a^3*b^2*c*k^3 +
60*a^2*b^3*c*k^3 - 19*a*b^4*c*k^3 + 12*a^3*b*c^2*k^3 +
98*a^2*b^2*c^2*k^3 + 52*a*b^3*c^2*k^3 + 2*a^3*c^3*k^3 + 60*a^2*b*c^3*k^3
+ 52*a*b^2*c^3*k^3 + 6*a^2*c^4*k^3 - 19*a*b*c^4*k^3 - 3*a*c^5*k^3 +
2*a^3*b^2*k^4 + 15*a^2*b^3*k^4 - 14*a*b^4*k^4 + 40*a^2*b^2*c*k^4 -
28*a*b^3*c*k^4 + 2*a^3*c^2*k^4 + 40*a^2*b*c^2*k^4 + 24*a*b^2*c^2*k^4 +
15*a^2*c^3*k^4 - 28*a*b*c^3*k^4 - 14*a*c^4*k^4 - 8*a^3*b*k^5 -
7*a^2*b^2*k^5 - 29*a*b^3*k^5 - 8*a^3*c*k^5 - 40*a^2*b*c*k^5 -
53*a*b^2*c*k^5 - 7*a^2*c^2*k^5 - 53*a*b*c^2*k^5 - 29*a*c^3*k^5 -
8*a^3*k^6 - 38*a^2*b*k^6 - 38*a*b^2*k^6 - 38*a^2*c*k^6 - 68*a*b*c*k^6 -
38*a*c^2*k^6 - 22*a^2*k^7 - 30*a*b*k^7 - 30*a*c*k^7 - 10*a*k^8 -
3*a^2*b^3*c^3 - 12*a*b^4*c^3 - 12*a*b^3*c^4 - 3*b^4*c^4 -
21*a^2*b^3*c^2*k - 25*a*b^4*c^2*k + 3*b^5*c^2*k - 21*a^2*b^2*c^3*k -
90*a*b^3*c^3*k - 9*b^4*c^3*k - 25*a*b^2*c^4*k - 9*b^3*c^4*k +
3*b^2*c^5*k + 4*a^3*b^2*c*k^2 - 30*a^2*b^3*c*k^2 - 8*a*b^4*c*k^2 +
6*b^5*c*k^2 + 4*a^3*b*c^2*k^2 - 60*a^2*b^2*c^2*k^2 - 134*a*b^3*c^2*k^2 +
6*b^4*c^2*k^2 - 30*a^2*b*c^3*k^2 - 134*a*b^2*c^3*k^2 - 24*b^3*c^3*k^2 -
8*a*b*c^4*k^2 + 6*b^2*c^4*k^2 + 6*b*c^5*k^2 + 4*a^3*b^2*k^3 -
12*a^2*b^3*k^3 + 5*a*b^4*k^3 + 3*b^5*k^3 + 16*a^3*b*c*k^3 -
30*a^2*b^2*c*k^3 - 22*a*b^3*c*k^3 + 27*b^4*c*k^3 + 4*a^3*c^2*k^3 -
30*a^2*b*c^2*k^3 - 122*a*b^2*c^2*k^3 + 12*b^3*c^2*k^3 - 12*a^2*c^3*k^3 -
22*a*b*c^3*k^3 + 12*b^2*c^3*k^3 + 5*a*c^4*k^3 + 27*b*c^4*k^3 + 3*c^5*k^3
+ 12*a^3*b*k^4 + 9*a^2*b^2*k^4 + 34*a*b^3*k^4 + 15*b^4*k^4 +
12*a^3*c*k^4 + 75*a^2*b*c*k^4 + 94*a*b^2*c*k^4 + 60*b^3*c*k^4 +
9*a^2*c^2*k^4 + 94*a*b*c^2*k^4 + 60*b^2*c^2*k^4 + 34*a*c^3*k^4 +
60*b*c^3*k^4 + 15*c^4*k^4 + 8*a^3*k^5 + 75*a^2*b*k^5 + 107*a*b^2*k^5 +
33*b^3*k^5 + 75*a^2*c*k^5 + 246*a*b*c*k^5 + 93*b^2*c*k^5 + 107*a*c^2*k^5
+ 93*b*c^2*k^5 + 33*c^3*k^5 + 54*a^2*k^6 + 138*a*b*k^6 + 42*b^2*k^6 +
138*a*c*k^6 + 84*b*c*k^6 + 42*c^2*k^6 + 60*a*k^7 + 30*b*k^7 + 30*c*k^7 +
9*k^8 + 6*a*b^4*c^2 + 18*a*b^3*c^3 + 9*b^4*c^3 + 6*a*b^2*c^4 + 9*b^3*c^4
+ 6*a*b^4*c*k - 3*b^5*c*k + 12*a^2*b^2*c^2*k + 79*a*b^3*c^2*k +
6*b^4*c^2*k + 79*a*b^2*c^3*k + 60*b^3*c^3*k + 6*a*b*c^4*k + 6*b^2*c^4*k
- 3*b*c^5*k - 3*b^5*k^2 + 62*a*b^3*c*k^2 - 30*b^4*c*k^2 +
150*a*b^2*c^2*k^2 + 36*b^3*c^2*k^2 + 62*a*b*c^3*k^2 + 36*b^2*c^3*k^2 -
30*b*c^4*k^2 - 3*c^5*k^2 - 12*a^2*b^2*k^3 + a*b^3*k^3 - 27*b^4*k^3 -
60*a^2*b*c*k^3 - 15*a*b^2*c*k^3 - 102*b^3*c*k^3 - 12*a^2*c^2*k^3 -
15*a*b*c^2*k^3 - 72*b^2*c^2*k^3 + a*c^3*k^3 - 102*b*c^3*k^3 - 27*c^4*k^3
- 60*a^2*b*k^4 - 92*a*b^2*k^4 - 87*b^3*k^4 - 60*a^2*c*k^4 -
294*a*b*c*k^4 - 258*b^2*c*k^4 - 92*a*c^2*k^4 - 258*b*c^2*k^4 -
87*c^3*k^4 - 48*a^2*k^5 - 223*a*b*k^5 - 156*b^2*k^5 - 223*a*c*k^5 -
330*b*c*k^5 - 156*c^2*k^5 - 130*a*k^6 - 147*b*k^6 - 147*c*k^6 - 54*k^7 -
6*a*b^3*c^2 - 9*b^4*c^2 - 6*a*b^2*c^3 - 21*b^3*c^3 - 9*b^2*c^4 -
6*a*b^3*c*k + 3*b^4*c*k - 66*a*b^2*c^2*k - 69*b^3*c^2*k - 6*a*b*c^3*k -
69*b^2*c^3*k + 3*b*c^4*k + 9*b^4*k^2 + 12*a^2*b*c*k^2 - 25*a*b^2*c*k^2 +
21*b^3*c*k^2 - 25*a*b*c^2*k^2 - 78*b^2*c^2*k^2 + 21*b*c^3*k^2 +
9*c^4*k^2 + 12*a^2*b*k^3 + 29*a*b^2*k^3 + 57*b^3*k^3 + 12*a^2*c*k^3 +
152*a*b*c*k^3 + 177*b^2*c*k^3 + 29*a*c^2*k^3 + 177*b*c^2*k^3 +
57*c^3*k^3 + 12*a^2*k^4 + 159*a*b*k^4 + 177*b^2*k^4 + 159*a*c*k^4 +
423*b*c*k^4 + 177*c^2*k^4 + 124*a*k^5 + 252*b*k^5 + 252*c*k^5 + 120*k^6
+ 3*b^4*c + 15*b^3*c^2 + 15*b^2*c^3 + 3*b*c^4 + 6*a*b^2*c*k + 15*b^3*c*k
+ 6*a*b*c^2*k + 84*b^2*c^2*k + 15*b*c^3*k - 9*b^3*k^2 - 24*a*b*c*k^2 +
12*b^2*c*k^2 + 12*b*c^2*k^2 - 9*c^3*k^2 - 42*a*b*k^3 - 60*b^2*k^3 -
42*a*c*k^3 - 186*b*c*k^3 - 60*c^2*k^3 - 48*a*k^4 - 177*b*k^4 - 177*c*k^4
- 120*k^5 - 3*b^3*c - 6*b^2*c^2 - 3*b*c^3 - 18*b^2*c*k - 18*b*c^2*k +
3*b^2*k^2 + 6*b*c*k^2 + 3*c^2*k^2 + 42*b*k^3 + 42*c*k^3 + 48*k^4 
d0:  -3*a*b*c - 3*a*b*k - 3*a*c*k - 3*a*k^2 + 3*b*c + 3*b*k + 3*c*k + 3*k^2 - 3*b - 3*c - 6*k 

d1:  a^3*b*c + a^3*b*k + a^3*c*k + a^3*k^2 - 3*a^2*b*c - 3*a^2*b*k - 3*a^2*c*k - 3*a^2*k^2 + 3*a^2*b + 3*a^2*c + 8*a*b*c + 6*a^2*k + 8*a*b*k + 8*a*c*k + 8*a*k^2 - 9*a*b - 9*a*c - 6*b*c - 18*a*k - 6*b*k - 6*c*k - 6*k^2 + 6*a + 9*b + 9*c + 18*k - 12 

e1:  -3*b^2*k - 3*c^2*k - 6*a*k^2 - 6*b*k^2 - 6*c*k^2 - 6*k^3 + 3*b*k + 3*c*k + 12*k^2 

e2:  3*a*b^2*c^2 + 3*a^2*b^2*k + 6*a*b^2*c*k + 3*a^2*c^2*k + 6*a*b*c^2*k + 2*a^3*k^2 + 6*a^2*b*k^2 + 3*a*b^2*k^2 + 6*a^2*c*k^2 + 12*a*b*c*k^2 + 3*a*c^2*k^2 + 6*a^2*k^3 + 6*a*b*k^3 + 6*a*c*k^3 + 3*a*k^4 - 3*a*b^2*c - 3*a*b*c^2 - 3*b^2*c^2 - 3*a^2*b*k - 9*a*b^2*k - 3*a^2*c*k - 12*a*b*c*k - 6*b^2*c*k - 9*a*c^2*k - 6*b*c^2*k - 12*a^2*k^2 - 21*a*b*k^2 - 3*b^2*k^2 - 21*a*c*k^2 - 12*b*c*k^2 - 3*c^2*k^2 - 18*a*k^3 - 6*b*k^3 - 6*c*k^3 - 3*k^4 + 3*a*b*c + 6*b^2*c + 6*b*c^2 + 15*a*b*k + 12*b^2*k + 15*a*c*k + 24*b*c*k + 12*c^2*k + 34*a*k^2 + 30*b*k^2 + 30*c*k^2 + 24*k^3 - 3*b^2 - 9*b*c - 3*c^2 - 12*a*k - 30*b*k - 30*c*k - 48*k^2 + 3*b + 3*c + 24*k 

f1:  3*a*b^4*c^2*k + 3*a*b^2*c^4*k + 6*a*b^4*c*k^2 + 6*a^2*b^2*c^2*k^2 + 12*a*b^3*c^2*k^2 + 12*a*b^2*c^3*k^2 + 6*a*b*c^4*k^2 + 3*a*b^4*k^3 + 12*a^2*b^2*c*k^3 + 24*a*b^3*c*k^3 + 12*a^2*b*c^2*k^3 + 36*a*b^2*c^2*k^3 + 24*a*b*c^3*k^3 + 3*a*c^4*k^3 + 6*a^2*b^2*k^4 + 12*a*b^3*k^4 + 24*a^2*b*c*k^4 + 48*a*b^2*c*k^4 + 6*a^2*c^2*k^4 + 48*a*b*c^2*k^4 + 12*a*c^3*k^4 + 12*a^2*b*k^5 + 21*a*b^2*k^5 + 12*a^2*c*k^5 + 48*a*b*c*k^5 + 21*a*c^2*k^5 + 6*a^2*k^6 + 18*a*b*k^6 + 18*a*c*k^6 + 6*a*k^7 - 3*a*b^3*c^2*k - 3*b^4*c^2*k - 3*a*b^2*c^3*k - 3*b^2*c^4*k - 6*a*b^3*c*k^2 - 6*b^4*c*k^2 - 30*a*b^2*c^2*k^2 - 12*b^3*c^2*k^2 - 6*a*b*c^3*k^2 - 12*b^2*c^3*k^2 - 6*b*c^4*k^2 - 3*a*b^3*k^3 - 3*b^4*k^3 - 51*a*b^2*c*k^3 - 24*b^3*c*k^3 - 51*a*b*c^2*k^3 - 36*b^2*c^2*k^3 - 3*a*c^3*k^3 - 24*b*c^3*k^3 - 3*c^4*k^3 - 24*a*b^2*k^4 - 12*b^3*k^4 - 84*a*b*c*k^4 - 48*b^2*c*k^4 - 24*a*c^2*k^4 - 48*b*c^2*k^4 - 12*c^3*k^4 - 39*a*b*k^5 - 21*b^2*k^5 - 39*a*c*k^5 - 48*b*c*k^5 - 21*c^2*k^5 - 18*a*k^6 - 18*b*k^6 - 18*c*k^6 - 6*k^7 + 3*b^4*c*k + 6*b^3*c^2*k + 6*b^2*c^3*k + 3*b*c^4*k + 3*b^4*k^2 + 6*a*b^2*c*k^2 + 24*b^3*c*k^2 + 6*a*b*c^2*k^2 + 42*b^2*c^2*k^2 + 24*b*c^3*k^2 + 3*c^4*k^2 + 6*a*b^2*k^3 + 18*b^3*k^3 + 24*a*b*c*k^3 + 84*b^2*c*k^3 + 6*a*c^2*k^3 + 84*b*c^2*k^3 + 18*c^3*k^3 + 18*a*b*k^4 + 48*b^2*k^4 + 18*a*c*k^4 + 120*b*c*k^4 + 48*c^2*k^4 + 12*a*k^5 + 57*b*k^5 + 57*c*k^5 + 24*k^6 - 3*b^3*c*k - 6*b^2*c^2*k - 3*b*c^3*k - 3*b^3*k^2 - 27*b^2*c*k^2 - 27*b*c^2*k^2 - 3*c^3*k^2 - 21*b^2*k^3 - 66*b*c*k^3 - 21*c^2*k^3 - 42*b*k^4 - 42*c*k^4 - 24*k^5 

f2:  -a^3*b^4*c^2*k^2 - a^3*b^2*c^4*k^2 - 2*a^3*b^4*c*k^3 - 2*a^4*b^2*c^2*k^3 - 4*a^3*b^3*c^2*k^3 - 4*a^3*b^2*c^3*k^3 - 2*a^3*b*c^4*k^3 - a^3*b^4*k^4 - 4*a^4*b^2*c*k^4 - 8*a^3*b^3*c*k^4 - 4*a^4*b*c^2*k^4 - 12*a^3*b^2*c^2*k^4 - 8*a^3*b*c^3*k^4 - a^3*c^4*k^4 - 2*a^4*b^2*k^5 - 4*a^3*b^3*k^5 - 8*a^4*b*c*k^5 - 16*a^3*b^2*c*k^5 - 2*a^4*c^2*k^5 - 16*a^3*b*c^2*k^5 - 4*a^3*c^3*k^5 - 4*a^4*b*k^6 - 7*a^3*b^2*k^6 - 4*a^4*c*k^6 - 16*a^3*b*c*k^6 - 7*a^3*c^2*k^6 - 2*a^4*k^7 - 6*a^3*b*k^7 - 6*a^3*c*k^7 - 2*a^3*k^8 - 3*a^2*b^4*c^4 - 2*a^3*b^4*c^2*k - 12*a^2*b^4*c^3*k - 2*a^3*b^2*c^4*k - 12*a^2*b^3*c^4*k - 4*a^3*b^4*c*k^2 - 7*a^3*b^3*c^2*k^2 - 15*a^2*b^4*c^2*k^2 - 7*a^3*b^2*c^3*k^2 - 48*a^2*b^3*c^3*k^2 - 4*a^3*b*c^4*k^2 - 15*a^2*b^2*c^4*k^2 - 2*a^3*b^4*k^3 - 14*a^3*b^3*c*k^3 - 6*a^2*b^4*c*k^3 - 10*a^3*b^2*c^2*k^3 - 60*a^2*b^3*c^2*k^3 - 14*a^3*b*c^3*k^3 - 60*a^2*b^2*c^3*k^3 - 2*a^3*c^4*k^3 - 6*a^2*b*c^4*k^3 - 7*a^3*b^3*k^4 - 7*a^3*b^2*c*k^4 - 24*a^2*b^3*c*k^4 - 7*a^3*b*c^2*k^4 - 72*a^2*b^2*c^2*k^4 - 7*a^3*c^3*k^4 - 24*a^2*b*c^3*k^4 - 2*a^3*b^2*k^5 + 12*a^3*b*c*k^5 - 24*a^2*b^2*c*k^5 - 2*a^3*c^2*k^5 - 24*a^2*b*c^2*k^5 + 9*a^3*b*k^6 + 3*a^2*b^2*k^6 + 9*a^3*c*k^6 + 3*a^2*c^2*k^6 + 6*a^3*k^7 + 6*a^2*b*k^7 + 6*a^2*c*k^7 + 3*a^2*k^8 + 3*a^2*b^4*c^3 + 3*a^2*b^3*c^4 + 6*a*b^4*c^4 + 2*a^3*b^3*c^2*k + 15*a^2*b^4*c^2*k - 3*a*b^5*c^2*k + 2*a^3*b^2*c^3*k + 24*a^2*b^3*c^3*k + 21*a*b^4*c^3*k + 15*a^2*b^2*c^4*k + 21*a*b^3*c^4*k - 3*a*b^2*c^5*k + 4*a^3*b^3*c*k^2 + 18*a^2*b^4*c*k^2 - 6*a*b^5*c*k^2 + 12*a^3*b^2*c^2*k^2 + 66*a^2*b^3*c^2*k^2 + 10*a*b^4*c^2*k^2 + 4*a^3*b*c^3*k^2 + 66*a^2*b^2*c^3*k^2 + 72*a*b^3*c^3*k^2 + 18*a^2*b*c^4*k^2 + 10*a*b^2*c^4*k^2 - 6*a*b*c^5*k^2 + 2*a^3*b^3*k^3 + 6*a^2*b^4*k^3 - 3*a*b^5*k^3 + 12*a^3*b^2*c*k^3 + 60*a^2*b^3*c*k^3 - 19*a*b^4*c*k^3 + 12*a^3*b*c^2*k^3 + 98*a^2*b^2*c^2*k^3 + 52*a*b^3*c^2*k^3 + 2*a^3*c^3*k^3 + 60*a^2*b*c^3*k^3 + 52*a*b^2*c^3*k^3 + 6*a^2*c^4*k^3 - 19*a*b*c^4*k^3 - 3*a*c^5*k^3 + 2*a^3*b^2*k^4 + 15*a^2*b^3*k^4 - 14*a*b^4*k^4 + 40*a^2*b^2*c*k^4 - 28*a*b^3*c*k^4 + 2*a^3*c^2*k^4 + 40*a^2*b*c^2*k^4 + 24*a*b^2*c^2*k^4 + 15*a^2*c^3*k^4 - 28*a*b*c^3*k^4 - 14*a*c^4*k^4 - 8*a^3*b*k^5 - 7*a^2*b^2*k^5 - 29*a*b^3*k^5 - 8*a^3*c*k^5 - 40*a^2*b*c*k^5 - 53*a*b^2*c*k^5 - 7*a^2*c^2*k^5 - 53*a*b*c^2*k^5 - 29*a*c^3*k^5 - 8*a^3*k^6 - 38*a^2*b*k^6 - 38*a*b^2*k^6 - 38*a^2*c*k^6 - 68*a*b*c*k^6 - 38*a*c^2*k^6 - 22*a^2*k^7 - 30*a*b*k^7 - 30*a*c*k^7 - 10*a*k^8 - 3*a^2*b^3*c^3 - 12*a*b^4*c^3 - 12*a*b^3*c^4 - 3*b^4*c^4 - 21*a^2*b^3*c^2*k - 25*a*b^4*c^2*k + 3*b^5*c^2*k - 21*a^2*b^2*c^3*k - 90*a*b^3*c^3*k - 9*b^4*c^3*k - 25*a*b^2*c^4*k - 9*b^3*c^4*k + 3*b^2*c^5*k + 4*a^3*b^2*c*k^2 - 30*a^2*b^3*c*k^2 - 8*a*b^4*c*k^2 + 6*b^5*c*k^2 + 4*a^3*b*c^2*k^2 - 60*a^2*b^2*c^2*k^2 - 134*a*b^3*c^2*k^2 + 6*b^4*c^2*k^2 - 30*a^2*b*c^3*k^2 - 134*a*b^2*c^3*k^2 - 24*b^3*c^3*k^2 - 8*a*b*c^4*k^2 + 6*b^2*c^4*k^2 + 6*b*c^5*k^2 + 4*a^3*b^2*k^3 - 12*a^2*b^3*k^3 + 5*a*b^4*k^3 + 3*b^5*k^3 + 16*a^3*b*c*k^3 - 30*a^2*b^2*c*k^3 - 22*a*b^3*c*k^3 + 27*b^4*c*k^3 + 4*a^3*c^2*k^3 - 30*a^2*b*c^2*k^3 - 122*a*b^2*c^2*k^3 + 12*b^3*c^2*k^3 - 12*a^2*c^3*k^3 - 22*a*b*c^3*k^3 + 12*b^2*c^3*k^3 + 5*a*c^4*k^3 + 27*b*c^4*k^3 + 3*c^5*k^3 + 12*a^3*b*k^4 + 9*a^2*b^2*k^4 + 34*a*b^3*k^4 + 15*b^4*k^4 + 12*a^3*c*k^4 + 75*a^2*b*c*k^4 + 94*a*b^2*c*k^4 + 60*b^3*c*k^4 + 9*a^2*c^2*k^4 + 94*a*b*c^2*k^4 + 60*b^2*c^2*k^4 + 34*a*c^3*k^4 + 60*b*c^3*k^4 + 15*c^4*k^4 + 8*a^3*k^5 + 75*a^2*b*k^5 + 107*a*b^2*k^5 + 33*b^3*k^5 + 75*a^2*c*k^5 + 246*a*b*c*k^5 + 93*b^2*c*k^5 + 107*a*c^2*k^5 + 93*b*c^2*k^5 + 33*c^3*k^5 + 54*a^2*k^6 + 138*a*b*k^6 + 42*b^2*k^6 + 138*a*c*k^6 + 84*b*c*k^6 + 42*c^2*k^6 + 60*a*k^7 + 30*b*k^7 + 30*c*k^7 + 9*k^8 + 6*a*b^4*c^2 + 18*a*b^3*c^3 + 9*b^4*c^3 + 6*a*b^2*c^4 + 9*b^3*c^4 + 6*a*b^4*c*k - 3*b^5*c*k + 12*a^2*b^2*c^2*k + 79*a*b^3*c^2*k + 6*b^4*c^2*k + 79*a*b^2*c^3*k + 60*b^3*c^3*k + 6*a*b*c^4*k + 6*b^2*c^4*k - 3*b*c^5*k - 3*b^5*k^2 + 62*a*b^3*c*k^2 - 30*b^4*c*k^2 + 150*a*b^2*c^2*k^2 + 36*b^3*c^2*k^2 + 62*a*b*c^3*k^2 + 36*b^2*c^3*k^2 - 30*b*c^4*k^2 - 3*c^5*k^2 - 12*a^2*b^2*k^3 + a*b^3*k^3 - 27*b^4*k^3 - 60*a^2*b*c*k^3 - 15*a*b^2*c*k^3 - 102*b^3*c*k^3 - 12*a^2*c^2*k^3 - 15*a*b*c^2*k^3 - 72*b^2*c^2*k^3 + a*c^3*k^3 - 102*b*c^3*k^3 - 27*c^4*k^3 - 60*a^2*b*k^4 - 92*a*b^2*k^4 - 87*b^3*k^4 - 60*a^2*c*k^4 - 294*a*b*c*k^4 - 258*b^2*c*k^4 - 92*a*c^2*k^4 - 258*b*c^2*k^4 - 87*c^3*k^4 - 48*a^2*k^5 - 223*a*b*k^5 - 156*b^2*k^5 - 223*a*c*k^5 - 330*b*c*k^5 - 156*c^2*k^5 - 130*a*k^6 - 147*b*k^6 - 147*c*k^6 - 54*k^7 - 6*a*b^3*c^2 - 9*b^4*c^2 - 6*a*b^2*c^3 - 21*b^3*c^3 - 9*b^2*c^4 - 6*a*b^3*c*k + 3*b^4*c*k - 66*a*b^2*c^2*k - 69*b^3*c^2*k - 6*a*b*c^3*k - 69*b^2*c^3*k + 3*b*c^4*k + 9*b^4*k^2 + 12*a^2*b*c*k^2 - 25*a*b^2*c*k^2 + 21*b^3*c*k^2 - 25*a*b*c^2*k^2 - 78*b^2*c^2*k^2 + 21*b*c^3*k^2 + 9*c^4*k^2 + 12*a^2*b*k^3 + 29*a*b^2*k^3 + 57*b^3*k^3 + 12*a^2*c*k^3 + 152*a*b*c*k^3 + 177*b^2*c*k^3 + 29*a*c^2*k^3 + 177*b*c^2*k^3 + 57*c^3*k^3 + 12*a^2*k^4 + 159*a*b*k^4 + 177*b^2*k^4 + 159*a*c*k^4 + 423*b*c*k^4 + 177*c^2*k^4 + 124*a*k^5 + 252*b*k^5 + 252*c*k^5 + 120*k^6 + 3*b^4*c + 15*b^3*c^2 + 15*b^2*c^3 + 3*b*c^4 + 6*a*b^2*c*k + 15*b^3*c*k + 6*a*b*c^2*k + 84*b^2*c^2*k + 15*b*c^3*k - 9*b^3*k^2 - 24*a*b*c*k^2 + 12*b^2*c*k^2 + 12*b*c^2*k^2 - 9*c^3*k^2 - 42*a*b*k^3 - 60*b^2*k^3 - 42*a*c*k^3 - 186*b*c*k^3 - 60*c^2*k^3 - 48*a*k^4 - 177*b*k^4 - 177*c*k^4 - 120*k^5 - 3*b^3*c - 6*b^2*c^2 - 3*b*c^3 - 18*b^2*c*k - 18*b*c^2*k + 3*b^2*k^2 + 6*b*c*k^2 + 3*c^2*k^2 + 42*b*k^3 + 42*c*k^3 + 48*k^4 
for kk in [1,2,3,5,6]: for aa in [2,3,4,7,10]: for bb in [1,4,11,14,15]: for cc in [2,6,7,10,16]: from_graph = Kemeny(B(kk,aa,bb,cc)) from_theory = educated_guess(kk,aa,bb,cc) if from_graph != from_theory: print "OOPS!", kk,aa,bb,cc print "Qapla!" # Qapla is Klingon for success 
       
Qapla!
Qapla!
# Note that there will be a disagreement if you plug in a=1 # this is due to the structure of the graphs used to compute # the Kemeny constant as being different in that case. Other # cases are all fine. 
       
# the following is used to generate the formula for D(a,b,c,d) a,b,c,d,x=var('a,b,c,d,x') def pa(n): return (-1/2)^(n-2-a)*(n+3)*(n+2)*(n+1)*n*(n-1)/120*x^2+(-1/2)^(n-1-a)*(n+2)*(n+1)*n*x/6+(-1/2)^(n-a)*(n+1) def pb(n): return (-1/2)^(n-2-b)*(n+3)*(n+2)*(n+1)*n*(n-1)/120*x^2+(-1/2)^(n-1-b)*(n+2)*(n+1)*n*x/6+(-1/2)^(n-b)*(n+1) r = -2*(x-1/c)*(x-1/d)/((c+2)*(d+2)) # long cycle; we have pulled out (-1/2)^(a+b) which is inside the pa, pb functions #first row r += pa(a) *pb(b) *(x-1/c)*(x-1)*(x-1/d)*(x-1) r += pa(a) *pb(b) *(x-1/c)*(x-1)*(-1/(d+2)) r += pa(a) *pb(b-1)*(x-1/c)*(x-1)*(x-1/d)*(-1/(2*(d+2))) r += pa(a-1)*pb(b) *(x-1/c)*(x-1)*(x-1/d)*(-1/(2*(d+2))) #second row r += pa(a) *pb(b) *(-1/(c+2))*(x-1/d)*(x-1) r += pa(a) *pb(b) *(-1/(c+2))*(-1/(d+2)) r += pa(a) *pb(b-1)*(-1/(c+2))*(x-1/d)*(-1/(2*(d+2))) r += pa(a-1)*pb(b) *(-1/(c+2))*(x-1/d)*(-1/(2*(d+2))) #third row r += pa(a) *pb(b-1)*(x-1/c)*(-1/(2*(c+2)))*(x-1/d)*(x-1) r += pa(a) *pb(b-1)*(x-1/c)*(-1/(2*(c+2)))*(-1/(d+2)) r += pa(a) *pb(b-2)*(x-1/c)*(-1/(2*(c+2)))*(x-1/d)*(-1/(2*(d+2))) r += pa(a-1)*pb(b-1)*(x-1/c)*(-1/(2*(c+2)))*(x-1/d)*(-1/(2*(d+2))) #fourth row r += pa(a-1)*pb(b) *(x-1/c)*(-1/(2*(c+2)))*(x-1/d)*(x-1) r += pa(a-1)*pb(b) *(x-1/c)*(-1/(2*(c+2)))*(-1/(d+2)) r += pa(a-1)*pb(b-1)*(x-1/c)*(-1/(2*(c+2)))*(x-1/d)*(-1/(2*(d+2))) r += pa(a-2)*pb(b) *(x-1/c)*(-1/(2*(c+2)))*(x-1/d)*(-1/(2*(d+2))) r1 = r.coefficients(x)[1][0].canonicalize_radical() *(c+2)*c*(d+2)*d*3 r1 = r1.canonicalize_radical().expand() r2 = r.coefficients(x)[2][0].canonicalize_radical() *(c+2)*c*(d+2)*d*3 r2 = r2.canonicalize_radical().expand() s = (r2*x+r1)*(-(c-1)*c*x+(c+1))*(-(d-1)*d*x+(d+1)) s1 = s.coefficients(x)[0][0].expand() s2 = s.coefficients(x)[1][0].expand() 
       
print "r1: ", r1, "\n" print "r2: ", r2, "\n" print "s1: ", s1, "\n" print "s2: ", s2, "\n" 
       
r1:  -3*a*c^2 - 3*b*c^2 - 3*a*d^2 - 3*b*d^2 - 6*a^2 - 12*a*b - 6*b^2 -
3*a*c - 3*b*c - 6*c^2 - 3*a*d - 3*b*d - 6*d^2 - 24*a - 24*b - 6*c - 6*d
- 24 

r2:  3*a*b*c^2*d^2 + a^3*c^2 + 3*a^2*b*c^2 + 3*a*b^2*c^2 + b^3*c^2 +
3*a*b*c^2*d + a^3*d^2 + 3*a^2*b*d^2 + 3*a*b^2*d^2 + b^3*d^2 +
3*a*b*c*d^2 + 3*a*c^2*d^2 + 3*b*c^2*d^2 + a^4 + 4*a^3*b + 6*a^2*b^2 +
4*a*b^3 + b^4 + a^3*c + 3*a^2*b*c + 3*a*b^2*c + b^3*c + 6*a^2*c^2 +
12*a*b*c^2 + 6*b^2*c^2 + a^3*d + 3*a^2*b*d + 3*a*b^2*d + b^3*d +
3*a*b*c*d + 6*a*c^2*d + 6*b*c^2*d + 6*a^2*d^2 + 12*a*b*d^2 + 6*b^2*d^2 +
6*a*c*d^2 + 6*b*c*d^2 + 3*c^2*d^2 + 8*a^3 + 24*a^2*b + 24*a*b^2 + 8*b^3
+ 12*a^2*c + 24*a*b*c + 12*b^2*c + 14*a*c^2 + 14*b*c^2 + 12*a^2*d +
24*a*b*d + 12*b^2*d + 9*a*c*d + 9*b*c*d + 9*c^2*d + 14*a*d^2 + 14*b*d^2
+ 9*c*d^2 + 23*a^2 + 46*a*b + 23*b^2 + 35*a*c + 35*b*c + 12*c^2 + 35*a*d
+ 35*b*d + 15*c*d + 12*d^2 + 28*a + 28*b + 30*c + 30*d + 12 

s1:  -3*a*c^3*d - 3*b*c^3*d - 3*a*c*d^3 - 3*b*c*d^3 - 3*a*c^3 - 3*b*c^3
- 6*a^2*c*d - 12*a*b*c*d - 6*b^2*c*d - 6*a*c^2*d - 6*b*c^2*d - 6*c^3*d -
6*a*c*d^2 - 6*b*c*d^2 - 3*a*d^3 - 3*b*d^3 - 6*c*d^3 - 6*a^2*c - 12*a*b*c
- 6*b^2*c - 6*a*c^2 - 6*b*c^2 - 6*c^3 - 6*a^2*d - 12*a*b*d - 6*b^2*d -
30*a*c*d - 30*b*c*d - 12*c^2*d - 6*a*d^2 - 6*b*d^2 - 12*c*d^2 - 6*d^3 -
6*a^2 - 12*a*b - 6*b^2 - 27*a*c - 27*b*c - 12*c^2 - 27*a*d - 27*b*d -
36*c*d - 12*d^2 - 24*a - 24*b - 30*c - 30*d - 24 

s2:  3*a*b*c^3*d^3 + a^3*c^3*d + 3*a^2*b*c^3*d + 3*a*b^2*c^3*d +
b^3*c^3*d + 6*a*b*c^3*d^2 + a^3*c*d^3 + 3*a^2*b*c*d^3 + 3*a*b^2*c*d^3 +
b^3*c*d^3 + 6*a*b*c^2*d^3 + 3*a*c^3*d^3 + 3*b*c^3*d^3 + a^3*c^3 +
3*a^2*b*c^3 + 3*a*b^2*c^3 + b^3*c^3 + a^4*c*d + 4*a^3*b*c*d +
6*a^2*b^2*c*d + 4*a*b^3*c*d + b^4*c*d + 2*a^3*c^2*d + 6*a^2*b*c^2*d +
6*a*b^2*c^2*d + 2*b^3*c^2*d + 6*a^2*c^3*d + 15*a*b*c^3*d + 6*b^2*c^3*d +
3*a*c^4*d + 3*b*c^4*d + 2*a^3*c*d^2 + 6*a^2*b*c*d^2 + 6*a*b^2*c*d^2 +
2*b^3*c*d^2 + 12*a*b*c^2*d^2 + 12*a*c^3*d^2 + 12*b*c^3*d^2 + a^3*d^3 +
3*a^2*b*d^3 + 3*a*b^2*d^3 + b^3*d^3 + 6*a^2*c*d^3 + 15*a*b*c*d^3 +
6*b^2*c*d^3 + 12*a*c^2*d^3 + 12*b*c^2*d^3 + 3*c^3*d^3 + 3*a*c*d^4 +
3*b*c*d^4 + a^4*c + 4*a^3*b*c + 6*a^2*b^2*c + 4*a*b^3*c + b^4*c +
2*a^3*c^2 + 6*a^2*b*c^2 + 6*a*b^2*c^2 + 2*b^3*c^2 + 6*a^2*c^3 +
12*a*b*c^3 + 6*b^2*c^3 + 3*a*c^4 + 3*b*c^4 + a^4*d + 4*a^3*b*d +
6*a^2*b^2*d + 4*a*b^3*d + b^4*d + 10*a^3*c*d + 30*a^2*b*c*d +
30*a*b^2*c*d + 10*b^3*c*d + 24*a^2*c^2*d + 54*a*b*c^2*d + 24*b^2*c^2*d +
17*a*c^3*d + 17*b*c^3*d + 6*c^4*d + 2*a^3*d^2 + 6*a^2*b*d^2 +
6*a*b^2*d^2 + 2*b^3*d^2 + 24*a^2*c*d^2 + 54*a*b*c*d^2 + 24*b^2*c*d^2 +
36*a*c^2*d^2 + 36*b*c^2*d^2 + 18*c^3*d^2 + 6*a^2*d^3 + 12*a*b*d^3 +
6*b^2*d^3 + 17*a*c*d^3 + 17*b*c*d^3 + 18*c^2*d^3 + 3*a*d^4 + 3*b*d^4 +
6*c*d^4 + a^4 + 4*a^3*b + 6*a^2*b^2 + 4*a*b^3 + b^4 + 9*a^3*c +
27*a^2*b*c + 27*a*b^2*c + 9*b^3*c + 24*a^2*c^2 + 48*a*b*c^2 + 24*b^2*c^2
+ 14*a*c^3 + 14*b*c^3 + 6*c^4 + 9*a^3*d + 27*a^2*b*d + 27*a*b^2*d +
9*b^3*d + 35*a^2*c*d + 73*a*b*c*d + 35*b^2*c*d + 82*a*c^2*d + 82*b*c^2*d
+ 15*c^3*d + 24*a^2*d^2 + 48*a*b*d^2 + 24*b^2*d^2 + 82*a*c*d^2 +
82*b*c*d^2 + 60*c^2*d^2 + 14*a*d^3 + 14*b*d^3 + 15*c*d^3 + 6*d^4 + 8*a^3
+ 24*a^2*b + 24*a*b^2 + 8*b^3 + 29*a^2*c + 58*a*b*c + 29*b^2*c +
70*a*c^2 + 70*b*c^2 + 12*c^3 + 29*a^2*d + 58*a*b*d + 29*b^2*d + 53*a*c*d
+ 53*b*c*d + 78*c^2*d + 70*a*d^2 + 70*b*d^2 + 78*c*d^2 + 12*d^3 + 23*a^2
+ 46*a*b + 23*b^2 + 39*a*c + 39*b*c + 60*c^2 + 39*a*d + 39*b*d + 27*c*d
+ 60*d^2 + 28*a + 28*b + 18*c + 18*d + 12 
r1:  -3*a*c^2 - 3*b*c^2 - 3*a*d^2 - 3*b*d^2 - 6*a^2 - 12*a*b - 6*b^2 - 3*a*c - 3*b*c - 6*c^2 - 3*a*d - 3*b*d - 6*d^2 - 24*a - 24*b - 6*c - 6*d - 24 

r2:  3*a*b*c^2*d^2 + a^3*c^2 + 3*a^2*b*c^2 + 3*a*b^2*c^2 + b^3*c^2 + 3*a*b*c^2*d + a^3*d^2 + 3*a^2*b*d^2 + 3*a*b^2*d^2 + b^3*d^2 + 3*a*b*c*d^2 + 3*a*c^2*d^2 + 3*b*c^2*d^2 + a^4 + 4*a^3*b + 6*a^2*b^2 + 4*a*b^3 + b^4 + a^3*c + 3*a^2*b*c + 3*a*b^2*c + b^3*c + 6*a^2*c^2 + 12*a*b*c^2 + 6*b^2*c^2 + a^3*d + 3*a^2*b*d + 3*a*b^2*d + b^3*d + 3*a*b*c*d + 6*a*c^2*d + 6*b*c^2*d + 6*a^2*d^2 + 12*a*b*d^2 + 6*b^2*d^2 + 6*a*c*d^2 + 6*b*c*d^2 + 3*c^2*d^2 + 8*a^3 + 24*a^2*b + 24*a*b^2 + 8*b^3 + 12*a^2*c + 24*a*b*c + 12*b^2*c + 14*a*c^2 + 14*b*c^2 + 12*a^2*d + 24*a*b*d + 12*b^2*d + 9*a*c*d + 9*b*c*d + 9*c^2*d + 14*a*d^2 + 14*b*d^2 + 9*c*d^2 + 23*a^2 + 46*a*b + 23*b^2 + 35*a*c + 35*b*c + 12*c^2 + 35*a*d + 35*b*d + 15*c*d + 12*d^2 + 28*a + 28*b + 30*c + 30*d + 12 

s1:  -3*a*c^3*d - 3*b*c^3*d - 3*a*c*d^3 - 3*b*c*d^3 - 3*a*c^3 - 3*b*c^3 - 6*a^2*c*d - 12*a*b*c*d - 6*b^2*c*d - 6*a*c^2*d - 6*b*c^2*d - 6*c^3*d - 6*a*c*d^2 - 6*b*c*d^2 - 3*a*d^3 - 3*b*d^3 - 6*c*d^3 - 6*a^2*c - 12*a*b*c - 6*b^2*c - 6*a*c^2 - 6*b*c^2 - 6*c^3 - 6*a^2*d - 12*a*b*d - 6*b^2*d - 30*a*c*d - 30*b*c*d - 12*c^2*d - 6*a*d^2 - 6*b*d^2 - 12*c*d^2 - 6*d^3 - 6*a^2 - 12*a*b - 6*b^2 - 27*a*c - 27*b*c - 12*c^2 - 27*a*d - 27*b*d - 36*c*d - 12*d^2 - 24*a - 24*b - 30*c - 30*d - 24 

s2:  3*a*b*c^3*d^3 + a^3*c^3*d + 3*a^2*b*c^3*d + 3*a*b^2*c^3*d + b^3*c^3*d + 6*a*b*c^3*d^2 + a^3*c*d^3 + 3*a^2*b*c*d^3 + 3*a*b^2*c*d^3 + b^3*c*d^3 + 6*a*b*c^2*d^3 + 3*a*c^3*d^3 + 3*b*c^3*d^3 + a^3*c^3 + 3*a^2*b*c^3 + 3*a*b^2*c^3 + b^3*c^3 + a^4*c*d + 4*a^3*b*c*d + 6*a^2*b^2*c*d + 4*a*b^3*c*d + b^4*c*d + 2*a^3*c^2*d + 6*a^2*b*c^2*d + 6*a*b^2*c^2*d + 2*b^3*c^2*d + 6*a^2*c^3*d + 15*a*b*c^3*d + 6*b^2*c^3*d + 3*a*c^4*d + 3*b*c^4*d + 2*a^3*c*d^2 + 6*a^2*b*c*d^2 + 6*a*b^2*c*d^2 + 2*b^3*c*d^2 + 12*a*b*c^2*d^2 + 12*a*c^3*d^2 + 12*b*c^3*d^2 + a^3*d^3 + 3*a^2*b*d^3 + 3*a*b^2*d^3 + b^3*d^3 + 6*a^2*c*d^3 + 15*a*b*c*d^3 + 6*b^2*c*d^3 + 12*a*c^2*d^3 + 12*b*c^2*d^3 + 3*c^3*d^3 + 3*a*c*d^4 + 3*b*c*d^4 + a^4*c + 4*a^3*b*c + 6*a^2*b^2*c + 4*a*b^3*c + b^4*c + 2*a^3*c^2 + 6*a^2*b*c^2 + 6*a*b^2*c^2 + 2*b^3*c^2 + 6*a^2*c^3 + 12*a*b*c^3 + 6*b^2*c^3 + 3*a*c^4 + 3*b*c^4 + a^4*d + 4*a^3*b*d + 6*a^2*b^2*d + 4*a*b^3*d + b^4*d + 10*a^3*c*d + 30*a^2*b*c*d + 30*a*b^2*c*d + 10*b^3*c*d + 24*a^2*c^2*d + 54*a*b*c^2*d + 24*b^2*c^2*d + 17*a*c^3*d + 17*b*c^3*d + 6*c^4*d + 2*a^3*d^2 + 6*a^2*b*d^2 + 6*a*b^2*d^2 + 2*b^3*d^2 + 24*a^2*c*d^2 + 54*a*b*c*d^2 + 24*b^2*c*d^2 + 36*a*c^2*d^2 + 36*b*c^2*d^2 + 18*c^3*d^2 + 6*a^2*d^3 + 12*a*b*d^3 + 6*b^2*d^3 + 17*a*c*d^3 + 17*b*c*d^3 + 18*c^2*d^3 + 3*a*d^4 + 3*b*d^4 + 6*c*d^4 + a^4 + 4*a^3*b + 6*a^2*b^2 + 4*a*b^3 + b^4 + 9*a^3*c + 27*a^2*b*c + 27*a*b^2*c + 9*b^3*c + 24*a^2*c^2 + 48*a*b*c^2 + 24*b^2*c^2 + 14*a*c^3 + 14*b*c^3 + 6*c^4 + 9*a^3*d + 27*a^2*b*d + 27*a*b^2*d + 9*b^3*d + 35*a^2*c*d + 73*a*b*c*d + 35*b^2*c*d + 82*a*c^2*d + 82*b*c^2*d + 15*c^3*d + 24*a^2*d^2 + 48*a*b*d^2 + 24*b^2*d^2 + 82*a*c*d^2 + 82*b*c*d^2 + 60*c^2*d^2 + 14*a*d^3 + 14*b*d^3 + 15*c*d^3 + 6*d^4 + 8*a^3 + 24*a^2*b + 24*a*b^2 + 8*b^3 + 29*a^2*c + 58*a*b*c + 29*b^2*c + 70*a*c^2 + 70*b*c^2 + 12*c^3 + 29*a^2*d + 58*a*b*d + 29*b^2*d + 53*a*c*d + 53*b*c*d + 78*c^2*d + 70*a*d^2 + 70*b*d^2 + 78*c*d^2 + 12*d^3 + 23*a^2 + 46*a*b + 23*b^2 + 39*a*c + 39*b*c + 60*c^2 + 39*a*d + 39*b*d + 27*c*d + 60*d^2 + 28*a + 28*b + 18*c + 18*d + 12 
def educated_guess2(a,b,c,d): return 1/3*(3*a*b*c^3*d^3 + a^3*c^3*d + 3*a^2*b*c^3*d + 3*a*b^2*c^3*d + b^3*c^3*d + 6*a*b*c^3*d^2 + a^3*c*d^3 + 3*a^2*b*c*d^3 + 3*a*b^2*c*d^3 + b^3*c*d^3 + 6*a*b*c^2*d^3 + 3*a*c^3*d^3 + 3*b*c^3*d^3 + a^3*c^3 + 3*a^2*b*c^3 + 3*a*b^2*c^3 + b^3*c^3 + a^4*c*d + 4*a^3*b*c*d + 6*a^2*b^2*c*d + 4*a*b^3*c*d + b^4*c*d + 2*a^3*c^2*d + 6*a^2*b*c^2*d + 6*a*b^2*c^2*d + 2*b^3*c^2*d + 6*a^2*c^3*d + 15*a*b*c^3*d + 6*b^2*c^3*d + 3*a*c^4*d + 3*b*c^4*d + 2*a^3*c*d^2 + 6*a^2*b*c*d^2 + 6*a*b^2*c*d^2 + 2*b^3*c*d^2 + 12*a*b*c^2*d^2 + 12*a*c^3*d^2 + 12*b*c^3*d^2 + a^3*d^3 + 3*a^2*b*d^3 + 3*a*b^2*d^3 + b^3*d^3 + 6*a^2*c*d^3 + 15*a*b*c*d^3 + 6*b^2*c*d^3 + 12*a*c^2*d^3 + 12*b*c^2*d^3 + 3*c^3*d^3 + 3*a*c*d^4 + 3*b*c*d^4 + a^4*c + 4*a^3*b*c + 6*a^2*b^2*c + 4*a*b^3*c + b^4*c + 2*a^3*c^2 + 6*a^2*b*c^2 + 6*a*b^2*c^2 + 2*b^3*c^2 + 6*a^2*c^3 + 12*a*b*c^3 + 6*b^2*c^3 + 3*a*c^4 + 3*b*c^4 + a^4*d + 4*a^3*b*d + 6*a^2*b^2*d + 4*a*b^3*d + b^4*d + 10*a^3*c*d + 30*a^2*b*c*d + 30*a*b^2*c*d + 10*b^3*c*d + 24*a^2*c^2*d + 54*a*b*c^2*d + 24*b^2*c^2*d + 17*a*c^3*d + 17*b*c^3*d + 6*c^4*d + 2*a^3*d^2 + 6*a^2*b*d^2 + 6*a*b^2*d^2 + 2*b^3*d^2 + 24*a^2*c*d^2 + 54*a*b*c*d^2 + 24*b^2*c*d^2 + 36*a*c^2*d^2 + 36*b*c^2*d^2 + 18*c^3*d^2 + 6*a^2*d^3 + 12*a*b*d^3 + 6*b^2*d^3 + 17*a*c*d^3 + 17*b*c*d^3 + 18*c^2*d^3 + 3*a*d^4 + 3*b*d^4 + 6*c*d^4 + a^4 + 4*a^3*b + 6*a^2*b^2 + 4*a*b^3 + b^4 + 9*a^3*c + 27*a^2*b*c + 27*a*b^2*c + 9*b^3*c + 24*a^2*c^2 + 48*a*b*c^2 + 24*b^2*c^2 + 14*a*c^3 + 14*b*c^3 + 6*c^4 + 9*a^3*d + 27*a^2*b*d + 27*a*b^2*d + 9*b^3*d + 35*a^2*c*d + 73*a*b*c*d + 35*b^2*c*d + 82*a*c^2*d + 82*b*c^2*d + 15*c^3*d + 24*a^2*d^2 + 48*a*b*d^2 + 24*b^2*d^2 + 82*a*c*d^2 + 82*b*c*d^2 + 60*c^2*d^2 + 14*a*d^3 + 14*b*d^3 + 15*c*d^3 + 6*d^4 + 8*a^3 + 24*a^2*b + 24*a*b^2 + 8*b^3 + 29*a^2*c + 58*a*b*c + 29*b^2*c + 70*a*c^2 + 70*b*c^2 + 12*c^3 + 29*a^2*d + 58*a*b*d + 29*b^2*d + 53*a*c*d + 53*b*c*d + 78*c^2*d + 70*a*d^2 + 70*b*d^2 + 78*c*d^2 + 12*d^3 + 23*a^2 + 46*a*b + 23*b^2 + 39*a*c + 39*b*c + 60*c^2 + 39*a*d + 39*b*d + 27*c*d + 60*d^2 + 28*a + 28*b + 18*c + 18*d + 12)/(a*c^3*d + b*c^3*d + a*c*d^3 + b*c*d^3 + a*c^3 + b*c^3 + 2*a^2*c*d + 4*a*b*c*d + 2*b^2*c*d + 2*a*c^2*d + 2*b*c^2*d + 2*c^3*d + 2*a*c*d^2 + 2*b*c*d^2 + a*d^3 + b*d^3 + 2*c*d^3 + 2*a^2*c + 4*a*b*c + 2*b^2*c + 2*a*c^2 + 2*b*c^2 + 2*c^3 + 2*a^2*d + 4*a*b*d + 2*b^2*d + 10*a*c*d + 10*b*c*d + 4*c^2*d + 2*a*d^2 + 2*b*d^2 + 4*c*d^2 + 2*d^3 + 2*a^2 + 4*a*b + 2*b^2 + 9*a*c + 9*b*c + 4*c^2 + 9*a*d + 9*b*d + 12*c*d + 4*d^2 + 8*a + 8*b + 10*c + 10*d + 8) 
       
for aa in [1,2,3,5,6]: for bb in [2,3,4,7,10]: for cc in [1,4,11,14,15]: for dd in [2,6,7,10,16]: from_graph = Kemeny(D(aa,bb,cc,dd)) from_theory = educated_guess2(aa,bb,cc,dd) if from_graph != from_theory: print "OOPS!", aa,bb,cc,dd print "Qapla!" # Qapla is Klingon for success 
       
Qapla!
Qapla!