from visual import *
from numpy import *

k=0
limit = 8
Cs=list()
Bs=list()
As=list()
superlist=list()
while limit > k:
    k+=1
    As.append(list())
    Cs.append(list())
    
k=0
As[1].append(1)
As[2].append(1)
As[2].append(2)

def gcd(m,n):
    while n:
        m,n=n,m%n
    return m

def lcm(a, b):
    return a*b/gcd(a,b)

def fac(N):
    m=1
    for x in xrange(1,int(N)):
        m=m*(int(N)-x+1)
    return m

def choose(a,b):
    return int(fac(a)/(fac(b)*fac(a-b)))

def intersection(list1,list2):
    outlist = list()
    for x in xrange(len(list2)):
        if list1.count(list2[x]) == 0:
            pass
        else:
            if outlist.count(list2[x]) == 0:
                outlist.append(list2[x])
            else:
                pass
    return outlist

def subtraction(list1,list2):               ## list 2 from list 1
    outlist = list()
    for x in xrange(len(list1)):
        outlist.append(list1[x])
    for x in xrange(len(list2)):
        if outlist.count(list2[x]) == 0:
            pass
        else:
            outlist.pop(outlist.index(list2[x]))
    return outlist



M=1234570        # up to this total ring size
limit = M
sumo = 1     # new total number of sets with property of commensurability
sump = 1     # new total number of sets with or without property of commensurability

set_ratio = 1
set_with = 0                # average number of sets each new set contains
net_without = 1             # average number of sets each previous set contained

agreement_ratio = 1
property_ratio = 1
kappa = 1
beta = 1

                            # Accumulated average of ratios of combinatorial sets with the property of commensurability to total sets within sets.
                            # This is equivalent to finding, those sets that agree (or ones contained in each other), in combinatorial number.
                            # This ratio is the number of sets without this property to those with and without this property
                            # These are the contained and containing because it is the number with the agreement condition to those with or without it.
                            # Each set takes part in both processes of containment and being contained.  The total is the number of sets.
                            # This is the proportion of sets that are neither contained nor the container.
                            # This is with the agreement condition.

totalwithproperty = 1
total = 1
k = 0                       # for the average



while limit-1 > k:
    k+=1
    if As[k].count(k) == 0:
        As[k].append(k)
    else:
        pass
    for x in xrange(len(As[k-1])):
        for y in xrange(len(As[k-1])):
            if As[k-1][x] != 1 and As[k-1][y]:
                b=As[k-1][x]*As[k-1][y]
                if As[k].count(b) == 0:
                    As[k].append(b)
                    if superlist.count(b) == 0:
                        superlist.append(b)
                else:
                    pass
                As[k].sort()
    As[k].sort()
    
k=0
while limit-2 > k:
    k+=1
##    print(k," and ",k+1,"lists")
##    print(As[k])
##    print(As[k+1])
##    print("intersection")
##    print(intersection(As[k],As[k+1]))
##    print("subtraction")
##    print(subtraction(As[k+1],intersection(As[k+1],As[k])))
    r = subtraction(As[k+1],intersection(As[k+1],As[k]))
    if Bs.count(r[0]) == 0 and superlist.count(r[0]) == 0:
       Bs.append(r[0])


for l in xrange(1,M+1):
    for n in xrange(1,int((l+1)/2)+1):
        k += 1
        
        l2=lcm(n,l)/n
        l1=lcm(l-n+1,l)/(l-n+1)

        a = choose(l,l-n-1)
        b = choose(l,n-1)

        if n == (l+1)/2:        
            total += 1
        else:
            total += 2
            
        sump += (a*b)
        
        if l1==l2 and n != (l+1)/2:
            totalwithproperty += 2
            sumo += 2*(a*b)
        if l1==l2 and n == (l+1)/2:
            totalwithproperty += 1
            sumo += (a*b)
##        if l1 == l2:
##            print("LCM",l1,l2,"equal for:",n,"&",l-n+1,"of",l)

        kappa = float(sump-sumo)/float(sump)
        beta = float(totalwithproperty)/float(total)
        
        set_ratio = set_ratio + kappa
        property_ratio = property_ratio + beta

        print("Property Ratio (pi):",1/(property_ratio/k),"Set Ratio (pi*?):",1/(set_ratio/k),"Their Ratio (?):",set_ratio/property_ratio)
        print(Bs[k])


print(Bs)
print(len(Bs))
