Number Theory Inhibitor Calculator

Question

I've been working on generating primes and prime products quickly to aid me in my research on prime numbers, their density, etc. The answer to my Large Number Limit Extravaganza question proved to be extremely useful in writing my programs. I have now expanded on this program to look at the difference in the amount of primes up to a certain point according to a lower bound, and the amount of "straight inhibitors" up to a certain point. A straight inhibitor is a figment of my imagination that I have used in my math and is extremely important for it.

import math
grandtotal = 1000
def prime_sieve(limit):
    a = [True] * limit
    a[0] = a[1] = False

    for i, isprime in enumerate(a):
        if isprime:
            yield i
            for n in xrange(i * i, limit, i):
                a[n] = False

def getDotArrayLog():
    e=55
    data = []
    while e<grandtotal:
        e+=1
        data.append((2 * e / (math.log(2.0 * e, 10) + 2)) - e / (math.log(e, 10) - 4))
    return data

def getDotArrayProduct():
    e=55
    data = []
    product = 1.0
    while True:
        if e==grandtotal: break
        e+=1
        for p in prime_sieve(e):
            product *= (1.0 - 1.0 / p)
        data.append(product)
        product = 1.0
    return data


if __name__ == "__main__":
    logs = getDotArrayLog()
    producta = getDotArrayProduct()
    print "Tested all the way up to", grandtotal
    counter = 0
    for e in logs:
        a = (1-producta[counter])*(counter+55)
        print "Amount of primes:", e
        print "Amount of inhibitors:", a
        print e-a
        counter+=1

I'm looking for some optimizations, condensations, and probably a great lecture on proper conventions.

Answer 1

2 small enhancements that I can see:

First, you could avoid to compute log(e,10) twice in each iteration by noticing that log(a*b) = log(a)+log(b). Since one term is fixed, you save 1 computation:

log2 = math.log(2.0,10)

def getDotArrayLog():
    e=55
    data = []

    while e<grandtotal:
        e+=1
        loge = math.log(e, 10)  # compute once only
        data.append((2 * e / (loge + log2 + 2)) - e / (loge - 4))
    return data

then, in your final loop, you could zip both arrays instead of using a counter and index access, and use enumerate to do the index counting for you, since you still need it for your computations:

for counter,(e,p) in enumerate(zip(logs,producta)):
    a = (1-p)*(counter+55)
    print "Amount of primes:", e
    print "Amount of inhibitors:", a
    print e-a

Answer 2

Now that you are using enumerate and zip, thanks to @Jean-FrancoisFabre, you can make your functions generators instead of building a list and returning it at the end:

log2 = math.log(2.0,10)

def getDotArrayLog(start, end
    for e in xrange(start, end):
        loge = math.log(e, 10) # compute once only
        yield (2 * e / (loge + log2 + 2)) - e / (loge - 4)

I also made the while loop into a for loop over a range.

Similarly for the product:

def getDotArrayProduct(start, end):
    for e in xrange(start, end):
        product = 1.
        for p in prime_sieve(e):
            product *= (1.0 - 1.0 / p)
        yield product

Note that I also made grandtotal and 55 + 1 a parameter of the functions.

This last function has one efficiency problem, though: It recomputes all primes (and their product) every loop, even though it is just the last product multiplied by the primes larger than e - 1. For this we can use generator.send and make our generator a co-routine. This makes it more complicated again, though:

def primes_less_than(limit):
    primes = prime_sieve(limit)
    p = next(primes)
    while True:
        n = yield
        if p < n:
            n = yield p
            p = next(primes)


def getDotArrayProduct(start, end):
    primes = primes_less_than(end)
    primes.send(None)
    product = 1.
    for e in xrange(start, end):
        while True:
            p = primes.send(e)
            next(primes)
            if p is not None:
                product *= (1.0 - 1.0 / p)
            else:
                break
        yield product


if __name__ == "__main__":
    start = 55
    end = 10000
    logs = getDotArrayLog(start + 1, end)
    producta = getDotArrayProduct(start + 1, end)
    print "Tested all the way up to", end
    for counter, (e, p) in enumerate(zip(logs, producta), start=start):
        a = (1 - p) * counter
        print "Amount of primes:", e
        print "Amount of inhibitors:", a
        print e - a

This has a generator called primes_less_than(total). The first thing it does is use the prime sieve to generate all the primes (and keep them in memory). Next it yields all primes less than the value being send to it.

The magic is in the n = yield, which receives the values sent to the generator which it then uses to check whether to yield the next prime or not.

Here an example:

>>> primes = primes_less_than(100)
>>> primes.send(None)
>>> primes.send(10)
2
>>> primes.send(10)
>>> primes.send(10)
3
>>> primes.send(3)
>>> primes.send(3)
None
>>> primes.send(10)
>>> primes.send(10)
5
>>> primes.send(10)
>>> primes.send(10)
7
>>> primes.send(10)
>>> primes.send(10)
None
>>> primes.send(20)
11
...

This way we can loop over all primes less than e and keep a running product. So we save iterating over all previous primes. Note that we have to either use send twice or put a next after every send.

Additionally, I made the magic number 55 a parameter of the function now and also used it to start the enumeration with the right value.