Posts Tagged ‘prime’

Problem 35 – Circular Primes

November 1, 2009

In problem 35 we are asked to find the number of “circular primes” below a bound. A prime is circular if all of the numbers obtained by cyclic permutation of the original number are also prime.

In my first solution, I’ll use my PrimeOracle from a while ago:

from PrimeOracle import PrimeOracle

def iscircularprime(num, isprime):
    if not isprime[num]:
        return False
    numstr = str(num)
    cyclics = [ int(numstr[i:] + numstr[:i]) for i in range(0,len(numstr)) ]
    return reduce(lambda a,b:a and b, map(lambda n: isprime[n], cyclics))

def solve(bound):
    """ Find how many circular primes there are below bound """
    isprime = PrimeOracle(bound)
    return len(set([n for n in range(2,bound) if iscircularprime(n, isprime)]))

This is pretty straightforward. I like how we can make all of the cyclic permutations using that one line (line 7), using the slice notation. And testing if all elements in a list are “True”, using reduce, is nice. I had hoped that maybe instead of the lambda expression I wrote, I could just pass in “and” as the first argument to that reduce. No dice though. Anybody know a nice way? I’ve also wondered about passing in the identity function in that first spot (in a filter, not a reduce), but don’t know how (without a lambda function or so).

Anyway, I think it’s good to have a solution like this that is easy to write and easy to read (hopefully) and, therefore, easy to check for correctness. Then I like to try to make it faster. The code above checks each cyclic permutation of each number many times, which is pointless.

Here’s a slightly faster solution, which doesn’t invoke the PrimeOracle:

def solve(bound):
    """ Find how many circular primes there are below bound """
    # isprime[n] if 2n+3 is prime, eventually
    # so to test if m is prime, check isprime[(m-3)/2]
    isprime = [True for n in range(0, bound//2)]
    for val in xrange(3, bound, 2):
        if isprime[(val-3)//2]:
            for mul in xrange(val*val, bound, 2*val):
                isprime[(mul-3)//2] = False

    ret = 1 # 2 isn't included in isprime, but it is circular
            # we've assumed bound > 2

    # let isprime[n] mean 2n+3 is a prime we haven't checked for circular primality

    And = lambda a,b: a and b # since we'll use it twice

    for val in xrange(3, bound, 2):
        if isprime[(val-3)//2]:
            valstr = str(val)
            circs = [ int(valstr[i:] + valstr[:i])
                      for i in range(0,len(valstr))]
            if reduce(And, map(lambda d:d%2, map(int, valstr))): # no even digits
                if reduce(And, map(lambda n:isprime[(n-3)//2], circs)): # all circs are prime
                    ret += len(set(circs))
            for n in circs:
                isprime[(n-3)//2] = False

    return ret

In line 25, where we increment ret, my original thought was to add len(valstr), since that’s how many cyclic permutations there are. However, these cyclic permutations need not be distinct, so I switched to making a set and getting its length.

What do you think of lines 23 and 24, an if within an if, and not much else? Would it be better (more readable?) to have simply “if () and ()”, even though that’d look best on two lines anyway? I guess python wouldn’t have to make one more nested context for that inner if, if we used 1 instead of 2…

It has, after writing all of this, occured to me that there’s some issues if “bound” isn’t a power of 10. Should we consider only primes such that all cyclic permutations are also primes less than the bound? We’d also have to think about how big to make our isprime array. So it would probably be better to re-write the above code accepting the bound as a number of digits.

Problem 27 – Prime Generating Quadratics

September 27, 2009

With my freshly minted PrimeOracle, I’m ready to attack problem 27, in which we are asked to find polynomials that generate long strings of consecutive primes starting when you plug in 0.

Basically my idea for this problem was to iterate through all of the choices for the coefficients a and b in n^2+an+b, and then for each polynomial plug in values until we hit a composite. Record maximums as appropriate.

There are a few easy optimizations to make. First, b itself has to be prime, because we test n=0. Next we test n=1, so we need 1+a+b\geq 2, that is a\geq 1-b. Those are the only optimizations I made, though I’m sure there are some more advanced ones out there.

So here’s my code:

# assume is in the same directory as this file
from PrimeOracle import PrimeOracle

def solve(n):
    # seed our oracle to 2*n+2, since n=1 gives 1+a+b
    isprime = PrimeOracle(2*n + 2)
    maxp, maxn = 41, 40
    for b in [m for m in xrange(2, n) if isprime[m] ]:
        for a in xrange(1-b, n):
            idx = 1 # we already tesetd 0 with isprime[b]
                idx += 1
            if(idx - 1 > maxn):
                maxp, maxn = a*b, idx - 1
    return maxp

I had my oracle also tell me what the biggest number it was asked about was, and the answer was up around 12000. So we could change the starting seed in the oracle to be big enough to never have to extend, and probably save some time.

Prime Oracle

September 27, 2009

I sat down to work on problem 27, and realized it might be handy to have a class that’ll tell me when numbers are prime. Something that’ll keep some memory of tested values, and sieve when it needs to know about bigger numbers than have been tested. Here’s what I threw together:

""" A class that'll help you decide if things are prime

Make an object of this class, I'll call it oracle, and
oracle[n] is True or False depending on if n is prime or not

Assumes negatives are not prime

class PrimeOracle:
    # Takes a starting maximum to sieve to
    def __init__(self, startmax = 2000):
        # self.isprime[n] = True if 2n+1 is prime
        self.isprime = [True for n in xrange(0,startmax//2 + 1)]
        self.isprime[0] = False
        # run the sieve
        for val in xrange(3, startmax + 1, 2):
            if(self.isprime[val // 2]):
                for mul in xrange(val*val, startmax + 1, 2*val):
                    self.isprime[mul // 2] = False

    # is n prime?
    def __getitem__(self, n):
        if n == 2:
            return True
        if n<2 or n % 2 == 0:
            return False
        nidx = n//2
        while(nidx >= len(self.isprime)):
            # extend the sieve
            # currently, we know the primality of 2*(len-1)+1=2*len-1
            # we can extend that to its square, 4*len^2 - 4len + 1
            # which would be index (4*len^2-4*len+1) // 2 = 4len(len-1)+1 // 2
            # so the new length should be one more than this
            # we already have length len, so we need to add (prev line) - len
            curlen = len(self.isprime)
            newlen = (4*curlen*(curlen-1) + 1) // 2
            self.isprime += [True for n in xrange(0, newlen+1-curlen)]
            maxknown = 2*curlen - 1
            maxtoknow = 2*newlen - 1
            for idx in xrange(1, curlen):
                    val = 2*idx + 1
                    # start crossing off multiples larger than maxknown
                    sfact = 1 + (maxknown // val)
                    if(sfact % 2 == 0):
                        sfact += 1 # otherwise our indexing is wrong
                    startfact = max(val, sfact)
                    for mul in xrange(startfact * val, maxtoknow + 1, 2*val):
                        self.isprime[mul // 2] = False
        return self.isprime[nidx]

    # the largest index we know about
    def __len__(self):
        return 2*len(self.isprime) - 1

    def __str__(self):
        return "\n".join(["%d: %s" % (2*n+1, self.isprime[n])
                          for n in xrange(0,len(self.isprime))])

It’s basically just the usual sieve, and then it extends itself whenever necessary.

There’s room for improvement. We could remember the highest value we’ve sieved to, and then take our array and just tease out the primes from it. Then we only store the list of primes. When we go to extend, we make a new array to hold whatever the new values would be, sieve appropriately, and then just append to our list of primes. Here’s my modified __getitem__ method:

def __getitem__(self, n):
    if n == 2:
        return True
    if n<2 or n % 2 == 0:
        return False
    while(n > self.sievedto):
        # we have sieved to sst = self.sievedto (assume odd)
        # we could sieve to sst^2 (will also be odd)
        # that'd be [sst+1 .. sst^2]
        # but really we might as well do [sst+2 .. sst^2]
        # sst=2s+1 => [2s+3 .. 4s^2+4s+1] //2 is [s+1 .. 2s(s+1)]
        # which has length 2s^2+s
        # and index i corresponds to 2(s+1+i)+1
        # so that an odd integer m is at m//2 - (s+1)
        s = self.sievedto // 2
        newsst = self.sievedto * self.sievedto
        isprime = [True for idx in xrange(0,2*s*s+s+1)]
        for pr in self.primes.keys():
            sfact = 1 + (self.sievedto // pr)
            if(sfact % 2 == 0):
                sfact += 1
            for mul in xrange(max(pr, sfact) * pr, newsst + 1, 2 * pr):
                isprime[mul//2 - (s+1)] = False
        for m in [2*(s+1+i)+1 for i in xrange(0,2*s*s+s+1) if isprime[i]]:
            self.primes[m] = True
        self.sievedto = newsst
    return self.primes.has_key(n)

My first guess was that the second version would be better, because it would be storing less. However, it seems to be a little slower. I guess filtering out the primes at each extension has its cost. Or my implementation could use some tweaking? I wish I knew a good way to compare the size requirements of these scripts, in addition to how long they take to run.

By the way, I would be a little shocked to find that I’m not off by one somewhere (many somewheres) in the above code. Little tests seem to work out. But use are your own peril.