Tuesday, December 11, 2012

Extensible sieve of Eratosthenes


I was creating a Carmichael triplets generator and needed a function to tell me if a number was prime. It wasn't immediately clear what the largest prime for checking would be so I decided to look into creating an extensible prime generator - If given a number to check for prime-ness that was outside its range, it would attempt to extend its range.

I stopped modifying the class as soon as it was working as the Carmichael generator task then ran in acceptable time.


class Isprime():
    '''
    Extensible sieve of Erastosthenes
    
    >>> isprime.check(11)
    True
    >>> isprime.multiples
    {2, 4, 6, 8, 9, 10}
    >>> isprime.primes
    [2, 3, 5, 7, 11]
    >>> isprime(13)
    True
    >>> isprime.multiples
    {2, 4, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 18, 20, 21, 22}
    >>> isprime.primes
    [2, 3, 5, 7, 11, 13, 17, 19]
    >>> isprime.nmax
    22
    >>> 
    '''
    multiples = {2}
    primes = [2]
    nmax = 2

    def __init__(self, nmax):
        if nmax > self.nmax:
            self.check(nmax)

    def check(self, n):
        if type(n) == float:
            if not n.is_integer(): return False
            n = int(n)
        multiples = self.multiples
        if n <= self.nmax:
            return n not in multiples
        else:
            # Extend the sieve
            primes, nmax = self.primes, self.nmax
            newmax = max(nmax*2, n)
            for p in primes:
                multiples.update(range(p*((nmax + p + 1) // p), newmax+1, p))
            for i in range(nmax+1, newmax+1):
                if i not in multiples:
                    primes.append(i)
                    multiples.update(range(i*2, newmax+1, i))
            self.nmax = newmax
            return n not in multiples

    __call__ = check



Now I wouldn't be surprised if some of the loop indices are not optimum and are maybe recomputing a few values that don't need to be, and this code is only efficient enough for my purposes; (a couple of Rosetta Code entries in other languages state that they are based on the Python but make changes so that is an indication).

1 comment:

  1. Some problems...

    In [72]: Isprime(2).check(41)
    Out[72]: True

    In [73]: isprime = Isprime(2)

    In [74]: isprime.
    isprime.__call__ isprime.__class__ isprime.__doc__ isprime.__init__ isprime.__module__ isprime.check isprime.multiples isprime.nmax isprime.primes

    In [74]: isprime.multiples
    Out[74]:
    set([2,
    4,
    6,
    8,
    9,
    10,
    12,
    14,
    15,
    16,
    18,
    20,
    21,
    22,
    24,
    25,
    26,
    27,
    28,
    30,
    32,
    33,
    34,
    35,
    36,
    38,
    39,
    40])

    In [75]: isprime.nmax
    Out[75]: 2

    In [76]: isprime.primes
    Out[76]: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]

    In [77]: isprime(43)
    Out[77]: True

    In [78]: isprime.nmax
    Out[78]: 43

    In [79]: isprime.primes
    Out[79]: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 3, 43]

    In [80]: isprime.multiples
    Out[80]:
    set([2,
    4,
    5,
    6,
    7,
    8,
    9,
    10,
    11,
    12,
    13,
    14,
    15,
    16,
    17,
    18,
    19,
    20,
    21,
    22,
    23,
    24,
    25,
    26,
    27,
    28,
    29,
    30,
    31,
    32,
    33,
    34,
    35,
    36,
    37,
    38,
    39,
    40,
    41,
    42])

    In [81]: isprime(41)
    Out[81]: False

    In [82]: isprime(29)
    Out[82]: False

    In [83]: isprime(23)
    Out[83]: False

    In [84]: isprime.nmax
    Out[84]: 43

    ReplyDelete