from visual import *
from numpy import *

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)))

M=5000        # up to this total ring size

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

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

        l1=lcm(n,l)/n
        l2=lcm(l-n+1,l)/(l-n+1)

        a = choose(l,n-1)
        b = choose(l,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("Total number of sets at size:",l,"is:",sump)
##    print("Total number of new sets at size:",l,"is:",sump-sumo)
##    print("Ratio of new contained sets to old contained sets:",float(new_set_contains)/float(previous_set_contains))        
    
##    Fi = F
##    F = 2*F+1+n*(n+1)/2
##    print("There are:",F,"sets of size:",n)
##    print("The ratio of previous to present is:")
    
##    G = G*int(choose(N,n))
##    print("There are:",G,"set chains at size:",n)
##    print("Their ratio is:",F/G,"sets at size:",n)
