#!/usr/bin/env sage -python ## written by Eric Sommers in May 2016 from sage.all import * R = PolynomialRing(QQ, 'q,z') R.inject_variables() """ Testing the positivity part of Thm 1.6 in our paper with Vic Reiner. This program tries to determine a factorization that shows that the ratio of q-numbers of form \prod [m-a_i]/[b_i] is a polynomial with positive coefficients when m belongs to the congruence classes modulo lcm(b_1, b_2, ..., b_r) when m is prime to the Coxeter number h. It makes use of Corollary 6 in Krattenthaler-Mueller's paper 'CYCLIC SIEVING FOR GENERALISED NON-CROSSING PARTITIONS ASSOCIATED WITH COMPLEX REFLECTION GROUPS OF EXCEPTIONAL TYPE - the details'. Namely, ([c]_q [ab]_q )/ ( [a]_q [b]_q ) has positive coefficients when gcd(a,b)=1 and c \geq (a-1)(b-1) This program therefore is less general than what is done in the Krattenthaler-Mueller paper since it can only handle cases where all denominator values, the b_i's, can be assigned to divide a numerator value in such a way that the sequence of denominator values, say, b_1, b_2, .., b_k satisifes b_i and lcm(b_{i+1}, b_{i+1}... b_k) do not divide each other for each i. See the first function below. """ def divide_lcm_sequence(k, a, j): b = list(a) b.insert(j, k) if len(b) == 1: return False for i in range(len(b)-1): ''' key constraint on this program so that we only have to use Corollary 6 in [KM]. ''' if (b[i] % lcm(b[i+1:])) == 0 or (lcm(b[i+1:]) % b[i]) == 0: return True return False def assign_to_spot_robust(a, m, number_being_assigned, num_index, j): if ((m- num_shifts[num_index]) % number_being_assigned) != 0: return False if len(a) == 0 or divide_lcm_sequence(number_being_assigned, a, j) == False: return True return False def new_initial_assignment(k, m): if k == -1: yield [[] for i in range(rank)] else: for choice in new_initial_assignment(k-1, m): # going to insert the k-th denom value for num_index in range(rank): #if assign_to_spot(choice[num_index], m, denom[k], num_index): current_len = len(choice[num_index]) new_array = list(choice[num_index]) new_choice = list(choice) for j in range(current_len+1): if assign_to_spot_robust(choice[num_index], m, denom[k], num_index, j): temp_array = list(new_array) new_array.insert(j, denom[k]) new_choice[num_index] = list(new_array) new_array = list(temp_array) yield new_choice def array_levels(a): return [len(a[i]) for i in range(rank)] def adjust_array_levels(old_levels, new_array, old_array): new_levels = [old_levels[i] + len(new_array[i]) - len(old_array[i]) for i in range(rank)] for i in range(rank): if old_levels[i] >= 2: new_levels[i] -= 1 return new_levels def next_pass(a, k, array_levels, m): ## assign gcds in kth position if k == -1: new_choice = list(a) yield new_choice else: for choice in next_pass(a, k-1, array_levels, m): if array_levels[k] <= 1: new_choice = list(choice) yield new_choice else: level = array_levels[k]-1 Gcd = gcd((choice[k])[level-1],lcm((choice[k])[level:])) for num_index in range(rank): if assign_to_spot_robust(choice[num_index], m, Gcd, num_index, len(choice[num_index])): new_array = list(choice[num_index]) + [Gcd] new_choice = list(choice) new_choice[num_index] = new_array yield new_choice def assign_divisors(k, a,b): if k == 0: yield [] else: for choice in assign_divisors(k-1, a, b): for index, i in enumerate(b): if index not in choice: if (a[k-1] % i) == 0: yield choice + [index] def check_validity(a): for i in range(rank): if len(a[i]) == 0: a[i] = [1] levels = array_levels(a) gcd_list = [] for index, list in enumerate(a): if len(list) > 1: for level in range(len(list)-1): value = gcd((list)[level],lcm((list)[level+1:])) gcd_list.append(value) if levels.count(1) < len(gcd_list): return False if len(gcd_list) == 0: return ([], []) ones_list = [(a[i][0], i) for i in range(rank) if len(a[i]) == 1 ] for i in assign_divisors(len(gcd_list), gcd_list, [j[0] for j in ones_list]): if len(i) == len(gcd_list): return ones_list, i return False def cyclo(n): return q**n-1 def cyclo_list(l): polynomial = 1 for i in l: polynomial *= cyclo(i) return polynomial def qCat(t): return cyclo_list(map(lambda x: t-x, num_shifts))/cyclo_list(denom) def check_pos_coeffs(poly): if poly == 0: return 1; if poly.denominator() != 1: #print "oops", poly.denominator() return 0 else: poly = poly.numerator() diff = len(poly.coefficients()); for i in range(diff): if poly.coefficients()[i] < 0: return 0 return 1 def print_out(a, b, c, m): ## c is assignmnet list, b is ones list gcd_counter = 0 voided_indices = [] m_min_values = [] for index, list in enumerate(a): if index in voided_indices: continue if len(list) > 1: for level in range(len(list)-1): value = gcd((list)[level],lcm((list)[level+1:])) if level == 0: print "[ (m-", num_shifts[index], ")/", lcm(list), "]_", lcm(list) this_index = b[c[gcd_counter]][1] print "([ (m-", num_shifts[this_index], ")/", a[this_index][0], "]_", print a[this_index][0], "[", lcm(list[level:])/value, "]_", value,")/", ## provide the minimum value of m that guarantees positivity from Cor 6 in [KM] m_min_values.append((list[level]/value - 1)*(lcm((list)[level+1:])/value -1) * value + num_shifts[this_index]) print "([", list[level]/value, "]_", value, "[", lcm((list)[level+1:])/value, "]_", value, ")", voided_indices.append(this_index) gcd_counter += 1 print for index, list in enumerate(a): if index in voided_indices or len(list) > 1: continue if len(list) == 1: print "[ (m-", num_shifts[index], ")/", a[index][0], "]_", a[index][0] else: print "[m-", num_shifts[index], "]_", 1 print ## now check small values for positivity, if needed if len(m_min_values) == 0: #print "No Further" return m_min = max(m_min_values) #if m_min < m: #print "No Further", m_min # qCat(m) for m_mod in range(m, m_min, shift): #print "Checking", m_mod #print qCat(m) if check_pos_coeffs(qCat(m)) == 0: print "Issue with value", m_mod ''' only goes 4 levels down. increase if needed ''' def poss_spots(): for m in range(shift): #if gcd(m,cox_num) != 1: ''' Throw out those congruence classes for which every element is not coprime to coxeter number ''' if gcd(m,gcd(cox_num, shift)) != 1: continue flg = False for jj in new_initial_assignment(rank-1, m): if flg != False: break exist = adjust_array_levels([0 for i in range(rank)], jj, [[] for i in range(rank)]) for kk in next_pass(jj, rank-1, exist, m): if flg != False: break exist2 = adjust_array_levels(exist, kk, jj) #if exist2 == exist: #flg = check_validity(jj) #print "Should stop" for ll in next_pass(kk, rank-1, exist2, m): if flg != False: break exist3 = adjust_array_levels(exist2, ll, kk) for mm in next_pass(ll, rank-1, exist3, m): flg = check_validity(mm) if flg != False: print "Win for m equal to", m #, mm print_out(mm, flg[0], flg[1], m) break if flg == False: print "Failed value of m:", m print all_cases_F4 = \ [ 12, [1,5,7,11], [2,6,8,12], [1,5,7], [2,4,6], [1,5], [2,2], #from G2_Cat [1,5], [2,6], #G2_Cat [1,3], [2,4] #B2 Cat ] all_cases_E6 = \ [ 12, #[1,4,5,7,8,11], [2,5,6,8,9,12], ## Catalan case won't work in general with the method in this program. ##See additional lemmas in [KM]. Still, we know by other reasons this is positive. #[1,4,5,7,8], [2,3,4,5,6], #[1,4,5,7], [1,2,4,6], #[1,4,5], [2,2,3], #A2 [1,4,5,3], [2,3,4,6], #fails 5, 11 #[1,4,5,7], [2,3,4,6], #fails 5,11 #[1,4,5,2], [2,3,4,6], #fails, 1,5,7, 11 (all) #[1,4,5,8], [2,3,4,6], ] all_cases_E7 = \ [ #18, 6, ## same prime factors as Coxeter number #[1,5,7,9,11,13,17], [2,6,8,10,12, 14, 18] , #[1,5,7,9,11,13], [2,4,6,6,8,10] , #[1,5,7,9,11], [2,2,4,6,8] , #[1,5,7,11], [2,6,8,12] ,#f4 #[1,5,7,9], [2,2,4,6] , #A2 one or the other pair except for 23 #[1,5,7,9,3], [2,4,6,6,10] , #[1,5,7,9,13], [2,4,6,6,10] , #[1,5,7,9,5], [2,4,6,6,10] , #[1,5,7,9,11], [2,4,6,6,10] , #A2+A1 #[1,5,7,7], [2,2,4,6], #works #[1,5,7,9], [2,2,4,6], #works #[1,5,7,5], [2,2,4,6], #works #[1,5,7,11], [2,2,4,6], #works, F4 #A4 #[1,3,5], [2,2,6] #B3 #[1,5,5], [2,2,6] #G2 works #A4+A1: clear ] all_cases_E8 = \ [ 30, #[1,7,11,13,17,19,23,29], [2,8,12,14,18,20,24,30], ## need more [1,7,11,13,17,19,23], [2,6,8,10,12, 14, 18], #[1,7,11,13,17,19], [2,4,6,8,10,12], #[1,7,11,13,17], [2,2,6,8,12], #[1,7,11,13,17, 23], [2,6,8,10,12, 18], ## funny case #[1,7,11,13], [2,4,6,8], #[1,7,11,13], [2,4,6,2], # follows from previous #[1,7,9,11,13], [2,4,6,8,10], # follows from previous #[1,7,11], [2,2,6], #[1,7,9,11], [2,4,6,2], #one missing is exact F4 #[5,7,11,13], [2,4,6,12], #[7,11, 13,17], [2,6,8,12], #[1,7,9], [2,2,4], #[1,5,7], [2,4,6], #[5,7], [2,6], #[7,11], [2,6], #[1,7,11], [2,2,6], ##A4+A1 works #[1,5,7], [2,2,6], #[1,7], [2,4], #[5,5,7], [2,4,6], #[5,7,11], [2,4,6], #[3,5], [2,4] ## from B_2 ## cases which are the sum of two ##2A2 has an F4 related term and #[1,7,11,17], [2,4,6,12], ## which works! # A2+A1 ## one pair or the other works except for 29 and 59 #[1,7,11,13, 11], [2,4,6,6,10], # fails #[1,7,11,13, 17], [2,4,6,6,10], # fails, should try other factorization #[1,7,11,13, 9], [2,4,6,6,10], # fails #[1,7,11,13, 19], [2,4,6,6,10], # fails in general #A2 ## one pair or the other works, except for 29 + 30s # [1,7,11,13,17,5], [2,6,8,10,12,18], # [1,7,11,13,17,23], [2,6,8,10,12,18], # [1,7,11,13,17,9], [2,6,8,10,12,18], # [1,7,11,13,17,19], [2,6,8,10,12,18], #A4 #fails for 23 and 53 seems like A2 in E7 #[1,7,9,5], [2,4,6,10] #[1,7,9,11], [2,4,6,10] #[1,7,9,3], [2,4,6,10] ] array = all_cases_F4 cox_num = cases = array[0] print for index, cases in enumerate(array): if index % 2 == 1: num_shifts = cases denom = array[index+1] #shift = lcm(denom+[cox_num]) shift = lcm(denom) rank = len(num_shifts) print "NEW CASE:" print "numerator shifts:", num_shifts print "denominator values:", denom print "LCM of denominator", shift print poss_spots()