Friday, February 21, 2014

Semiprimes in Python

Someone added a Semiprime task to Rosetta Code:
    > Semiprime numbers are natural numbers that are products of exactly two (possibly equal) prime numbers. 
    > Example: 1679 = 23 × 73 (This particular number was chosen as the length of the Arecibo message). 
    > Write a function determining whether a given number is semiprime.
I first looked for a prime generator and found one in part of an existing Prime decomposition tasks example code, where decompose is a generator for the prime factors of its argument.
In [ ]:
from __future__ import print_function

import sys
 
def is_prime(n):
    return zip((True, False), decompose(n))[-1][0]
 
class IsPrimeCached(dict):
    def __missing__(self, n):
        r = is_prime(n)
        self[n] = r
        return r
 
is_prime_cached = IsPrimeCached()
 
def primes():
    yield 2
    n = 3
    while n < sys.maxint - 2:
        yield n
        n += 2
        while n < sys.maxint - 2 and not is_prime_cached[n]:
            n += 2
 
def decompose(n):
    for p in primes():
        if p*p > n: break
        while n % p == 0:
            yield p
            n /=p
    if n > 1:
        yield n
 
if __name__ == '__main__':
    # Example: calculate factors of Mersenne numbers to M59 #
 
    import time
 
    for m in primes():
        p = 2 ** m - 1
        print( "2**{0:d}-1 = {0:d}, with factors:".format(m, p) )
        start = time.time()
        for factor in decompose(p):
            print (factor, end=' ')
            sys.stdout.flush()
 
        print( "=> {0:.2f}s".format( time.time()-start ) )
        if m >= 59:
            break

Semiprime

The code for semiprime kind'a wrote itself from the description. I have a generator of prime factors. Lets check if the first two prime factors multiply to give n for a semiprime. Any errors and it isn't:
In [8]
from prime_decomposition import decompose

def semiprime(n):
    d = decompose(n)
    try:
        return d.next() * d.next() == n
    except:
        return False
You can then play around in the interpreter with it:
In [9]:
semiprime(1679)
Out[9]:
True
In [10]:
[n for n in range(1,101) if semiprime(n)]
Out[10]:
[4,
 6,
 9,
 10,
 14,
 15,
 21,
 22,
 25,
 26,
 33,
 34,
 35,
 38,
 39,
 46,
 49,
 51,
 55,
 57,
 58,
 62,
 65,
 69,
 74,
 77,
 82,
 85,
 86,
 87,
 91,
 93,
 94,
 95]

No comments:

Post a Comment