12.py

'''
Author:    Michael Sherif Naguib
Date:      November 3, 2018
@:         University of Tulsa

Question #12:  What is the value of the first triangle number to have over five hundred divisors?
Example:     
    The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:

    1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

    Let us list the factors of the first seven triangle numbers:

    1: 1                    
    3: 1,3                  
    6: 1,2, 3,6             
    10: 1,2, 5,10
    15: 1,3, 5,15
    21: 1,3, 7,21
    28: 1,2 ,4,7 ,14,28

    We can see that 28 is the first triangle number to have over five divisors.
'''

'''
Ideas: 

generate triangle numbers and test
predict by looking at how many factors each number has.... ?
find how many factors (not necessarily what they are) and generate a base number test up.....
    observation: n^k (where n>1 k>0 n,k are integers)  has 1+k factors 
    what about n_1^k_1*n_i^k_i....

Note!
repeted factors of numbers are not counted however products of those factors are
ex. 2^2*3 ===> 1,2*2*3, 2*3,3,2*2  

n^k*m^i ==> (1+k)(1+i) factors.... minimize product, maximize factors ==> good starting point

2^2*3^2 = 36 ... 1,2,3,4,6,9,12,18,36 ==>9 factors  (1+k)(1+i)
3^3*5^1 = 135... 1,3,5,9,15,27,45,135 ==>8 factors (1+k)(1+i)
2^1*5^1*3^2 = 90                      ==>12 factors (predicted...)

the number of factors a composite number N is the product of 1 + each coefficient
in the prime factors of a number ... i.e    n^k*m^i*l^j  has (k+1)(i+1)(j+1) factors...
... all numbers k,i,j affect the divisors N has at an equal rate.... however
if n<m<l then j influences product size the most.... making it larger the fastest...
and n makeing it larger the slower...

 
init [i,j,k]  all=0
init [n,m,l] as increasing primes

iterate 
    if n^k+1>m^i
        check n^k*m^i*l^j
        k++
    else:
        if m^(i+1) <l^j
            i++
            check n^k*m^i*l^j
        else:
            if ...... etc ......
    k++
        
i.e ensure n^k<m^i<l^j holds before adding to i j k and add such as to maintain that condition
    increment and meet that condition... 



OR


upper and lower bounds for divisor function... calculate those for every triangle number
if the target number of divisiors is in the range... factor N and check divisors if divisors checks out
return N...except... the lower bound is 1 i.e the number is prime and the upper bound ? 

OR

 generate a list of numbers by knowing their prime factorization and calculate their divisors... 
grab all that have divisors >=500 ... order this subset.... check if it is a triangle num


'''
import math
ten = __import__("10")
eleven  = __import__("11")
isPrime = __import__("3").isPrime
#read a list of numbers from a file into an array
def numArrayFromFile(filename):
    #gotta love one-liners: read file, get list of lines, concat lines, replace "," with " " , split on those, convert items in that list to int...finish
    return list(map(lambda x: int(x),reduce(lambda x,y: y+x, open(filename,"r").readlines()).replace(","," ").split(" ")))

#returns the Nth triangle Number
def nthTriangleNum(N):
    return N*(N+1)/2# (1/2) N^2 +(1/2)N   ==> the triangle numbers increase at a rate of 2N

#tests if a number N where N>1, is a triangle number
def isTriangleNum(N):
    #note f(f+1)/2 ==> is gaurenteed to be a triangle num if f is a positive integer
    # N = f(f+1)/2 ==> 2N = f(f+1) ==> f^2 +f -2N = 0
    #if either or both of the roots are integers then the num is a triangle num because root(root+1)/2 = N
    des = 1+8*N # The Descriminant b^2 -4*a*c      
    root_one= -1 + math.sqrt(des)#    2*a = 2*0.5 = 1
    print
    #root_two ... -0.5 -math.sqrt(des)   ===> is always negative because Range sqrt [0,inf)
    return float(int(root_one)) == root_one

#creates a number given a list of primes and their coefficient powers
def make_num(primes,powers):
    num=1
    for i in range(0,len(primes)):
        p=primes[i]
        n = powers[i]
        if(n != 0):#skip lengthy calculation for p^n n=0 ==> p^n=1
            num = num*math.pow(p,n)
    return num

#takes a list of power coefficients corresponding t the ordered set of prime factors for a number and calculates the number of
#divisors that number has
def divisors(k_list):
    #  i=n(exclusive) i=0 II (k_i+1)
    return eleven.prod_list(list(map(lambda x: x+1,k_list)))

#increase according to the rule: n^k<m^i<l^j
def incrementIndex(index_i,p,k,PRIMES_UNDER):
    if index_i == len(p)-1:
        #code here that adds another prime to the list and an another 0 coeefficient to k
        errMsg = "Not Enough Primes: used all primes under "+str(PRIMES_UNDER)
        print(errMsg)
        raise errMsg
    elif math.pow(p[index_i],k[index_i]+1) > math.pow(p[index_i+1], k[index_i+1] if k[index_i+1] != 0 else 1):
        print( "{}^({}) > {}^{} ".format(p[index_i],k[index_i]+1,p[index_i+1],k[index_i+1] if k[index_i+1] != 0 else 1))
        incrementIndex(index_i+1,p,k,PRIMES_UNDER)
    else:
        print("incrementing: {}^({}+1)".format(p[index_i],k[index_i]))
        if(index_i!=0):
            print("decrementing: {}^({}+1)".format(p[index_i-1],k[index_i-1]))
            k[index_i-1] = k[index_i-1] -1
        k[index_i] = k[index_i] + 1

#returns prime factors and their coefficients for a number N    NOTE N musbe be less than the greatest prime in the list of primes
def primeFactorsCoefficients(n,list_of_primes):
    #based on code in three.py
    if(n==1):
        return []
    for prime in list_of_primes:
        coeff=0
        if n//prime== n/prime:
            while(n//prime== n/prime):# reduce as much as possible
                coeff = coeff+1
                n/=prime
            return [(prime,coeff)] + primeFactorsCoefficients(int(n),list_of_primes)# add that number to the list then return the recursion
    return [(n,1)] #otherwise n is prime

if __name__=="__main__":
    #calculate a lot of primes
    print("Calculating primes ... could be read from a file")
    primes_list = ten.sieveOfEratosthenes(1000000)
    counter=1

    #loop
    l=len(primes_list)
    while counter < l:
        #get the next triangle num
        triangle_num = nthTriangleNum(counter)
        #find its factorication
        pandk = primeFactorsCoefficients(triangle_num,primes_list)# p1^k1*p2^k2 ... pn^kn
        #determine its divisors count
        divisors_count = divisors(list(map(lambda x: x[1], pandk)))
        #print("divisors {}".format(divisors_count))
        #does it meet the criteria?
        if divisors_count>=500:
            print(triangle_num)
            break
        else:
            #continue searching
            counter = counter + 1
    print("end")

    '''#OLD ALGORITHM THAT DID NOT WORK
    #print(list(map(lambda x: isTriangleNum(nthTriangleNum(x)),list(range(1,10)))))#a test : confirm identity 
    #settings
    TARGET_NUM_OF_DIVISORS=6
    PRIMES_UNDER = 30

    #primes and their coefficients... arbritrary number of primes... if any number in that list n^k > last number^1 then more primes needed...
    primes_list = ten.sieveOfEratosthenes(PRIMES_UNDER)
    len_primes=len(primes_list)
    power_coeff = [0]*len_primes
    result=1#smallest triangle num

    #Calculate the result
    while(True):
        #only check if there are enough divisors...
        doCheck = divisors(power_coeff) >= TARGET_NUM_OF_DIVISORS
        if doCheck:
            print("did check")
            #What is the number formed from those Divisors?
            N = make_num(primes_list,power_coeff)
            #is it a triangle num?
            if isTriangleNum(N):
                print(str(N))
                result=N
                break
        #otherwise increment to the next number for which we know the divisors...
        #this is a recursive function which checks at the selected start index then continues
        #to combe over the lists...
        incrementIndex(0,primes_list,power_coeff,PRIMES_UNDER)                
    '''