#!/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()