Zeta File Creator | Python Fiddle

## Our goal is to answer the following question: Given any positive integer x, how many finite abelian groups exist with x elements (unique up to isomorphism)?
## We start out be saying that each x can be expressed as its unique prime factorization:
## x = p(0)^k(0) * ... * p(m)^k(m)
## where each p(j) is a unique prime, and each k(j) is a positive integer
## so we need to find the prime factorization for each x. in order to do that, we need to get prime numbers up to at least the square root of x
## From the fundamental theorem of finitely generated abelian groups we can easily show the following:
##
## Every finite abelian group G is isomorphic to a direct product of cyclic groups in the form Z(p(1)^r(1)) * ... * Z(p(n)^r(n))
## Where the p(j) are primes, not necessarily distinct, and the r(j) are positive integers.
## Note that each Z(p(j)^r(j)) is just modular arithmetic on the integers modulo (p(j)^r(j))
## That is, modular arithmetic where the modulus is some power of a prime number
##
## So, given our unique factorization of x from earlier, it can easily be shown that if we can create a function:
## psi(k), which, given a positive integer k, returns the number of finite abelian groups with p^k elements for any prime p,
## then the number of finite abelian groups with x elements, call this number zeta(x), will be equal to the product of the psi(k(j)) for all of the k(j) in the prime factorization of x.
##
## We want to create a list of the zeta(x) for all positive integers x up to 1000000 (or some other "large" number)
## We will then study the behavior of zeta(x)
##
## in order to do this, we obviously need the prime factorization of each x
## and to do that, we will need a function that generates primes.

def primegenerator(ceil):
    t=6
    primes=[2,3,5]
    while t<=ceil:
        z=0
        a=0
        k=primes[a]
        y=(t**(0.5)+1)
        h=len(primes)-1
        while k < y:
            if t%primes[a] == 0:
                z+=1
                k=y+1
            else:
                a+=1
                k=primes[a]
        if z==0:
            primes.append(t)
        t+=1
    return primes

## ok, now in order to create our psi(n) function from earlier, we are going to have to make another function called gamma
## gamma will be defined recursively, so in order to keep the runtime low, we're going to generate a file which contains lists of all the gamma values we will need
## we will then refer to this file when we need our gamma values later to calculate our psi values, rather than repeatedly calculate the gammas.
## bear with me. I will provide an exlanation of the gamma function soon.

def gfilemaker(n):
    glist = [0]
    gglist = str(glist)
    gfile = open('gammafile.txt','w')
    gfile.write(gglist)
    def gammaloader(s,k):
        a = len(glist)
        if s<a:
            return glist[s][k-1]
        elif k==1:
            return 1
        elif k==s:
            return 1
        else:
            y = 0
            for x in range(1,min(s-k+1,k+1)):
                y+=gammaloader(s-k,x)
            return y
        ## gammaloader is our gamma function.
        ## It calculates, recursively, the number of finite abelian groups of order p^s for any prime p, FOR WHICH THE HIGHEST ORDER SUBGROUP is of order p^k
        ## technically, the first "if" statement is unnecessary, because it can always be calculated recursively, but referring to the list limits the recursion
        ## this keeps the runtime down, otherwise it can take a LONG time
        ## given any positive integer s, we define psi(s)=sum([gamma(s,k) for k in range(1,s+1)])
        ## psi(s) is equal to the total number of finite abelian groups with p^s elements for ANY prime integer p.
    for a in range(1,n+1):
        alist = []
        for k in range(1,a+1):
            alist.append(gammaloader(a,k))
        glist.append(alist)
        aglist = str(alist)
        gfile.write('\n')
        gfile.write(aglist)
    gfile.close()

# gfilemaker(500)
## uncomment the above line and run it to create a file with all the gammaloader(s,k) for all positive integers s up to 500 (or whatever argument you want).

## This is more gammas than we would ever need when creating our list of zeta values, since we would never calculate zeta for all integers up to 2^500.
## Still, it will be nice to have these gamma values (and psi values) for the sake of studying the behavior of these functions. Who knows what insight these lists might provide.
##
## Take note of that definition of our psi function above, we will use it in the below function
## The following function will create a file containing a list of our zeta(x) values for all x up to n

def psiglistmaker(n):
    gammafile = open('gammafile.txt','r')
    ## importing our gamma file to use so that we don't have to go through insane levels of recursion
    lines = gammafile.readlines()
    def gline(a):
        gline=lines[a].replace('[','').replace(']','').replace(',','')
        gline = gline.split()
        for x in range(len(gline)):
            gline[x]=int(gline[x])
        return gline
    ## the gline function gets the line we need from the gamma file and returns it as a list
    def quickpsi(num):
        guesses=primegenerator(num**0.5)
        ## guesses is the list of primes that we will need. we only need to generate this list once
        def quickprimer(trial):
            primedict = {}
            def factor(p, s):
                if not(p in primedict):
                    primedict[p]=1
                else:
                    primedict[p]+=1
            for p in guesses:
                while trial % p == 0:
                    factor(p, trial)
                    trial = int(trial/p)
                if trial == 1:
                    return primedict
            if trial != 1:
                primedict[trial] = 1
            return primedict
        ## quickprimer returns a dictionary containing the prime factorization of the argument, where the keys are p(j) and values are k(j)
        ## remember, we need the prime factoriztion so we can calculate the psi values and ultimately the zeta value.
        psifactlist=[0]
        for z in range(1,num+1):
            b=1
            for a in quickprimer(z):
                c = sum(gline(quickprimer(z)[a]))
                ## HERE is our psi value: psi(x) the sum of the gamma(x,k)for all positive integers k up to x.
                b*=c
                ## And HERE is our zeta value.
                ## zeta(x) is equal to the product of the psi(k(j)) values for all the k(j) in our prime factorization of x
                ## So, this loop will ultimately set our variable b equal to the zeta(x)
                ## we will add zeta(x) to a running list of zetas
            psifactlist.append(b) 
        return psifactlist
    ## quickpsi returns a list of our zeta(x) for all positive integer x up to the argument.
    ## now all that is left to do is to store the list in a file.
    zetafile = open('zetafile.txt','w')
    zetalist = quickpsi(n)
    gammafile.close()
    zetatext = str(zetalist)
    gammafile.close()
    zetafile.write(zetatext)
    zetafile.close()
    return zetalist

#zetas = psiglistmaker(1000)
## uncomment the above line and run to create a list and a file containing a list of the zeta values up to 1000000 (or whatever argument you want)

## beware: when the argument was 1000000, this took a couple minutes on my laptop (which is a Lenovo ideapad U430. Not exactly a super powerful machine, but still).
## we can access this file later to study the behavior of zeta
## Furthermore, let sigma(n) = sum(zeta(x) for all positive integers x up to n)
## I am very interested to see the behavior of the function sigma(n)/n
## I want to know if this function has a limit as n approaches infinity
## I will use the file of zeta values that I created to study this sigma(n)/n function

Python Fiddle

Python Cloud IDE