Herausforderung angenommen
Posted by Jesco Freund at March 27, 2009 12:01 p.m.
André hat es provoziert – was er mit Perl hinkriegt, muss ich mit Python doch auch hinbekommen 
. Natürlich habe ich mich nicht damit begnügt, einen Abklatsch seiner Perl-Implementierung nach Python zu portieren – ein bisschen Optimierung musste schon sein. In Zeiten moderner Multicore-CPUs, so dachte ich, wäre es ideal, das Problem so zu formulieren, dass es parallelisiert abgearbeitet werden kann. Bei einem iterativen Algorithmus, dessen Iterationen auf die Ergebnisse der jeweiligen Vorstufen angewiesen sind, allerdings kein ganz triviales Unterfangen…
Zur Theorie
Dass es dennoch nicht unmöglich ist, liegt an einer Eigenschaft der natürlichen Zahlen bezüglich Multiplikation: Der kleinste Primfaktor einer natürlichen Zahl ist kleiner oder gleich der Quadratwurzel dieser Zahl (sofern sie überhaupt faktorisierbar ist). Nun ist das Ermitteln des ganzzahligen Anteils der Quadratwurzel einer sehr großen natürlichen Zahl (sagen wir mal, mit 600 Dezimalstellen oder so…) alles andere als erquicklich oder gar schnell erledigt. Dreht man den Spieß aber um, geht es nicht nur fixer, sondern man bekommt auch gleich ganze Zahlenreihen, die parallel abgearbeitet werden können.
Ich beginne daher in meinem Algorithmus mit zwei aufeinanderfolgenden Primzahlen p und q. Dabei setze ich voraus, dass alle Primzahlen ≤ p bekannt sind. Trivial wäre etwa ein Startpärchen (2, 3) – beides sind Primzahlen, es gibt keine Primzahl < 2 und zwischen 2 und 3 liegt definitiv keine weitere Primzahl. Betrachte ich nun das Intervall [p2; q2 – 1] (im Beispiel also [4; 8]), müssen alle Nicht-Primzahlen in diesem Intervall einen Primfaktor ≤ p aufweisen.
Da alle Primfaktoren ≤ p bereits bekannt sind, können alle Kandidaten im Intervall [p2; q2 – 1] parallel auf Primalität getestet werden. Bei kleinen p und q mag dies vielleicht kaum Geschwindigkeitsvorteile gegenüber einem linearen Test bringen – die Kandidatenliste wächst aber recht rasch, so dass der Parallelisierungsvorteil umso deutlicher ins Gewicht fällt, je größer p und q werden. Nach Abarbeitung des Intervalls [p2; q2 – 1] müssen p und q neu vorbelegt werden. pneu wird dabei zu qalt und qneu wird die kleinste gefundene Primzahl > pneu zugewiesen.
Zur Implementierung
Meine Implementierung dieses Algorithmus ist etwas umfassender geraten – ein paar Erläuterungen sollten also nicht schaden. Eine zentrale Rolle spielt dabei die Klasse PTStatus. Diese Klasse ist als Singleton ausgelegt und bildet den Zustand des Primzahl-Automaten ab – also die Belegung von p und q, die Liste aller Primzahlen ≤ p und die Liste aller gefundenen Primzahlen ≥ q . Aus diesen Informationen heraus generiert die Methode candidate_range eine iterierbare Folge von möglichen Kandidaten des Intervalls [p2; q2 – 1] , die auf Primalität getestet werden sollen. Bei der Überprüfung festgestellte Primzahlen können mit der Methode register der Liste der gefundenen Primzahlen hinzugefügt werden. Sind alle Kandidaten überprüft, kommt die Methode forward zum Einsatz: Sie verschiebt die Werte von p und q auf das nächste zu prüfende Intervall und stockt die Liste der prüfungsrelevanten Primzahlen entsprechend auf.
Da Singleton-Objekte, die mithilfe innerer Klassen implementiert wurden, von Pythons pickle-Modul nicht serialisiert werden können, besitzt die Klasse darüber hinaus eine Methode save, die den aktuellen Zustand des Objekts in eine Datei schreibt. Bei der Initialisierung des Status-Objekts wird versucht, den gespeicherten Zustand aus einer ggf. vorhandenen Sicherungsdatei auszulesen und wieder herzustellen, so dass der Automat die Arbeit an der Stelle wieder aufnehmen kann, bei der zuletzt abgebrochen wurde.
Die eigentliche Arbeit wird allerdings von den Worker Threads übernommen, also Instanzen der Klasse PrimeWorker. Diesen wird bei der Erstellung eine Referenz auf ein Queue-Objekt mitgegeben. Die Worker-Threads lesen nun aus dieser Queue mögliche Primzahl-Kandidaten und testen sie auf Primalität. Durch das Queue-Modell blieb es mir erspart, für jeden Thread ein festes Arbeitspaket zu schnüren. Der Vorteil darin: Jeder Thread holt sich dann die nächste Aufgabe, wenn er die vorhergehende abgearbeitet hat. So werden zu große Idle-Zeiten vermieden. Selbst bei ungünstigem Scheduling ist die Auslastung noch nahezu optimal: Bekommt ein Thread überoproportional viel Rechenzeit zugeteilt, erledigt er seine Aufgaben auch schneller und kann insgesamt mehr Kandidaten prüfen.
Alle übrigen Klassen und Funktionen sind mehr oder weniger Hilfskonstrukte und schmückendes Beiwerk zur besseren Handhabung. Ich möchte daher nur kurz auf zwei Kleinigkeiten eingehen. Da wie bereits oben angedeutet der Status des Automaten in Form einer Pseudo-Objektserialisierung gespeichert wird, gibt es keine Klartext-Ergebnisliste mit den gefundenen Primzahlen. Um diese dennoch erzeugen zu können, habe ich die Funktion dump eingebaut. Sie gibt alle im Statusobjekt enthaltenen Primzahlen als Auflistung auf der Konsole aus.
Außerdem möchte ich noch eine Bemerkung zur Obergrenze machen: diese kann zwar über eine Kommandozeilenoption gesetzt werden; das Programm indes wird in der Regel jedoch nicht exakt an dieser angegebenen Obergrenze abbrechen. Der Grund hierfür liegt in der Parallelisierung: Das Programm muss ein Quadrat-Intervall immer vollständig abarbeiten. Zwar verhindert das Status-Objekt das Hinzufügen von bereits bekannten Primzahlen, aber niemand kann garantieren, dass ein Thread noch an der Prüfung eines Kandidaten kaut, der kleiner als die vorgegebene Höchstgrenze ist, während ein anderer Thread einen Kandidaten zu fassen kriegt, der über der Höchstgrenze liegt und per Event alle anderen Threads beendet. Damit wäre die Vollständigkeit der Primzahl-Listen im Statusobjekt nicht mehr gewährleistet; diese ist jedoch Voraussetzung dafür, dass der Algorithmus mit den Quadratintervallen überhaupt funktionieren kann.
Quelltext
Doch nun genug der Worte; hier das komplette Skript (getestet mit Python 2.5 und 2.6 unter FreeBSD und Windows). Zu Lehrzwecken darf es gerne kopiert und weitergegeben werden, solange der Urheber dabei erkennbar bleibt.
#!/usr/bin/env python
# -*- coding: utf-8 -*-
'''
primaton - the prime number automaton
Copyright (c) 2009 Jesco Freund. All rights reserved.
This program generates prime numbers. It uses a modified "Sieve of Eratosthenes"
algorithm meaning it will find all prime numbers within a given range of
natural numbers, and it is deterministic (meaning all indicated primes are
definitively primes, not only on a probability base).
The implementation comes with some deviations of the original sieve. So it can
be run without indicating an upper limit for prime research, meaning it is only
limited by memory and processing time. To benefit from the power of modern
SMP systems, this implementation makes use of "square chunks" which can be
investigated concurrently without risk of mistake.
For a better understandig of these square chunks, please look at the following
explanation: It is a mathematical law that the biggest prime factor of a natural
number cannot be larger than its square root (clear thing - a factor larger than
a number`s square root would have a counterpart that is smaller than the sqrt).
So if p and q are subsequent prime numbers (e. g. 2 and 3), then one can
conclude that all naturals within [p**2, (q**2)-1] that cannot be devided
neither by p nor by any prime smaller than p are primes.
This means that for testing all naturals in [p**2, (q**2)-1], we only need to
know all primes smaller than q (which includes p). Therefore, investigation of
such a range can be parallelized without the risk of missing a non-prime due to
the fact that not all primes smaller than the candidate natural are already
known.
'''
# Stuff we need for multithreading
import threading
import Queue
# Other stuff we need...
import logging
import os
import pickle
import sys
from optparse import OptionParser
# Global stuff for multithreading
SuperLock = threading.Lock()
class PTConf(object):
'''
Configuration storage class providing some kind of singleton
Python dictionary with a simplified interface (only get() and set()
are implemented).
'''
class __impl(object):
'''
Inner class implementing the singleton behaviour
It simply provides a protected dictionary storing configuration
properties as key => value assignments. It can be accessed by
respective get() and set() methods
'''
def __init__(self):
'''
Class constructor creating the (empty) Python dictionary
'''
self.__conf = {}
def get(self, property):
'''
Method for reading a single value from the dictionary that
corresponds to the given key
'''
if property in self.__conf:
return self.__conf[property]
else:
return ''
def set(self, property, value):
'''
Method for storing a key => value pair within the dictionary.
If a pair with the same key already exists, the previously stored
value is overwritten by the new one.
'''
self.__conf[property] = value
# Storage for the singleton instance reference
__instance = None
def __init__(self):
'''
Constructor for the `outer` class
If no instance of the inner class exists, create a new one.
If there is already one, just do - nothing...
'''
if PTConf.__instance is None:
PTConf.__instance = PTConf.__impl()
# Store reference to implementation class instance
self.__dict__['_PTConf__instance'] = PTConf.__instance
def reload(self):
'''
This method resets the inner class and reinitialises it.
This behaviour can be used for implementing reactions
related to SIGHUP signals or similar.
'''
PTConf.__instance = None
self.__dict__['_PTConf__instance'] = None
self.__init__()
# We need delegating setattr and getattr methods, so we overwrite them
def __getattr__(self, attr):
return getattr(self.__instance, attr)
def __setattr__(self, attr, value):
return setattr(self.__instance, attr, value)
class PTLogger(object):
'''
Frontend class for the Python logging module
granting some singleton behaviour to serialize logging from
several threads into a single file without resource deadlocks
'''
class __impl(object):
'''
Inner class implementing the singleton behaviour
'''
def __init__(self):
# define loglevel dictionary
self.__level = {}
self.__level['DEBUG'] = logging.DEBUG
self.__level['INFO'] = logging.INFO
self.__level['WARNING'] = logging.WARNING
self.__level['ERROR'] = logging.ERROR
self.__level['CRITICAL'] = logging.CRITICAL
# configure Logger
conf = PTConf()
logging.basicConfig(
filename = conf.get('logfile'),
level = self.__level[conf.get('loglevel')],
format = conf.get('logformat'),
datefmt = conf.get('logdatefmt')
)
self.__logger = logging.getLogger('Primaton')
def log(self, level, message):
self.__logger.log(self.__level[level], message)
# Storage for the singleton instance reference
__instance = None
def __init__(self):
'''
Constructor for the `outer` class
If no instance of the inner class exists, create a new one.
If there is already one, just do - nothing...
'''
if PTLogger.__instance is None:
PTLogger.__instance = PTLogger.__impl()
# Store reference to implementation class instance
self.__dict__['_PTLogger__instance'] = PTLogger.__instance
def reload(self):
'''
This method resets the inner class and reinitialises it.
This behaviour can be used for implementing reactions
related to SIGHUP signals or similar.
'''
PTLogger.__instance = None
self.__dict__['_PTLogger__instance'] = None
self.__init__()
# We need delegating setattr and getattr methods, so we overwrite them
def __getattr__(self, attr):
return getattr(self.__instance, attr)
def __setattr__(self, attr, value):
return setattr(self.__instance, attr, value)
class PTStatus(object):
'''
Singleton object used to reflect the current operation status of
primaton. Using this object makes it possible to pause and later resume
prime calculation.
'''
class __impl(object):
'''
Inner class implementing the singleton behaviour
'''
def __init__(self):
'''
Class constructor defining starting conditions for prime research
'''
conf = PTConf()
lg = PTLogger()
try:
'''
Let`s have a go and try to load the stuff from file
'''
fobj = open(conf.get('primedb'), 'r')
readbuffer = ''
run = 'p'
# iterate over all lines in file
for line in fobj:
# in the first place, we expect P= and q= assignments
# then the relevant and last the irrelevant prime list pickle
if len(line) >= 2:
if run == 'p' and line[0:2] == 'p=':
self.__p = long(line.partition('=')[-1])
run = 'q'
continue
elif run == 'q' and line[0:2] == 'q=':
self.__q = long(line.partition('=')[-1])
run = 'rel'
continue
elif run == 'rel' and line[0:2] == 'a.':
run = 'irr'
readbuffer += line
self.__relevant_primes = pickle.loads(readbuffer)
lg.log('DEBUG', "Relevant Primes: %s" % str(self.__relevant_primes))
readbuffer = ''
continue
elif run == 'irr' and line[0:2] == 'a.':
readbuffer += line
self.__irrelevant_primes = pickle.loads(readbuffer)
lg.log('DEBUG', "Irrelevant Primes: %s" % str(self.__irrelevant_primes))
break
# anything else must belong to the prime list pickles
readbuffer += line
fobj.close()
lg.log('INFO', "Successfully loaded data from %s" % conf.get('primedb'))
except IOError, e:
lg.log('WARNING', "Could not open primes database. Starting at primordial ooze...")
self.__p = 2
self.__q = 3
self.__relevant_primes = [2,]
self.__irrelevant_primes = [3,]
def forward(self):
'''
This method is used to go one step ahead and turn
towards the next square chunk
'''
# move first irrelevant prime into relevant primes
self.__irrelevant_primes.sort()
self.__relevant_primes.append(self.__irrelevant_primes.pop(0))
self.__p = self.__q
self.__q = self.__irrelevant_primes[0]
def get_pq(self):
'''
This method returns a textual representation of p and q
'''
return "[%d, %d]" % (self.__p, self.__q)
def register(self, prime):
'''
register a newly found prime
'''
if prime < self.__p and prime not in self.__relevant_primes:
self.__relevant_primes.append(prime)
elif prime > self.__q and prime not in self.__irrelevant_primes:
self.__irrelevant_primes.append(prime)
def candidate_range(self):
'''
Generator yielding all candidates from the (p, q) square chunk
'''
i = self.__p ** 2
stop = self.__q ** 2 # only calculate this once and store result
if i % 2 == 0:
i += 1
while i < stop:
yield i
i += 2
def relevant_primerange(self):
'''
Generator yielding all relevant primes
'''
i = 0
while i < len(self.__relevant_primes):
yield self.__relevant_primes[i]
i += 1
def complete_primerange(self):
'''
Generator yielding all primes (relevant + irrelevant ones)
'''
i = 0
while i < len(self.__relevant_primes):
yield self.__relevant_primes[i]
i += 1
i = 0
while i < len(self.__irrelevant_primes):
yield self.__irrelevant_primes[i]
i += 1
def latest_prime(self):
'''
This method returns the largest prime that has been found so far
'''
return self.__irrelevant_primes[-1]
def save(self):
'''
This method dumps the current status into the database file
'''
conf = PTConf()
lg = PTLogger()
try:
fobj = open(conf.get('primedb'), 'w')
fobj.write("p=%d\n" % self.__p)
fobj.write("q=%d\n" % self.__q)
fobj.write("%s\n" % pickle.dumps(self.__relevant_primes))
fobj.write(pickle.dumps(self.__irrelevant_primes))
fobj.close()
except IOError, e:
lg.log('ERROR', "Could not open primes database. %d: %s" % (e.args[0], e.args[1]))
# Storage for the singleton instance reference
__instance = None
def __init__(self):
'''
Constructor for the `outer` class
If no instance of the inner class exists, create a new one.
If there is already one, just do - nothing...
'''
if PTStatus.__instance is None:
PTStatus.__instance = PTStatus.__impl()
# Store reference to implementation class instance
self.__dict__['_PTStatus__instance'] = PTStatus.__instance
def reload(self):
'''
This method resets the inner class and reinitialises it.
This behaviour can be used for implementing reactions
related to SIGHUP signals or similar.
'''
PTStatus.__instance = None
self.__dict__['_PTStatus__instance'] = None
self.__init__()
# We need delegating setattr and getattr methods, so we overwrite them
def __getattr__(self, attr):
return getattr(self.__instance, attr)
def __setattr__(self, attr, value):
return setattr(self.__instance, attr, value)
class PrimeWorker(threading.Thread):
'''
Worker thread class for prime evaluation
'''
def __init__(self, jobqueue, status):
'''
Class constructor. Call static constructor of parent class and set
up the job queue and the result set pointer
'''
threading.Thread.__init__(self)
self.jobqueue = jobqueue
self.status = status
def run(self):
'''
This is the method that is actually executed when running the thread
'''
# main loop of worker thread: read from job queue and process
# whatever came from there
while True:
# wait until there is something in the job queue for this thread
candidate = self.jobqueue.get()
# a negative value in the job queue acts like a SIGKILL
if candidate < 0:
break
if self.is_prime(candidate):
SuperLock.acquire()
self.status.register(candidate)
SuperLock.release()
self.jobqueue.task_done()
def is_prime(self, candidate):
'''
This function does the actual testing for primality.
It returns a boolean test result (True or False).
'''
for i in self.status.relevant_primerange():
if candidate % i == 0:
return False
return True
def die(message=''):
'''
Implementation of the poplar die() function shamelessly stolen from PHP
This function prints a message (if one has been given) and exits with
a status != 0
'''
if message != '':
sys.stderr.write("%s\n" % message)
sys.exit(1)
def daemonize():
'''
This function creates a daemon process (works only on Linux/Unix operating
systems). This makes it possible to run this program independently
of any terminal.
'''
conf = PTConf()
# first fork
try:
pid = os.fork()
if pid > 0:
sys.exit(0) # exit first parent
except OSError, e:
die("fork #1 failed: (%d) %s\n" % (e.errno, e.strerror))
# Set process environment
os.chdir("/")
os.umask(0)
os.setsid()
try:
pid = os.fork()
if pid > 0:
sys.exit(0) # Exit second parent.
except OSError, e:
die("fork #2 failed: (%d) %s\n" % (e.errno, e.strerror))
# Redirect standard file descriptors.
si = open('/dev/null', 'r')
so = open(conf.get("dumpfile"), 'a+')
se = open(conf.get("dumpfile"), 'a+', 0)
os.dup2(si.fileno(), sys.stdin.fileno())
os.dup2(so.fileno(), sys.stdout.fileno())
os.dup2(se.fileno(), sys.stderr.fileno())
# Write PID to specified file
try:
fobj = open(conf.get("pidfile"), "r")
die("Fatal error: %s already exists" % conf.get("pidfile"))
except IOError:
try:
fobj = open(conf.get("pidfile"), "w")
fobj.write("%d\n" % os.getpid())
fobj.close
except IOError:
die("Fatal error: cannot open %s for writing" % conf.get("pidfile"))
def dump():
'''
This function dumps the database content to stdout
'''
s = PTStatus()
for i in s.complete_primerange():
sys.stdout.write("%d\n" % i)
def main():
'''
The main function managing the whole stuff, including thread creation
and termination, queue handling etc.
'''
# CLI argument handling
parser = OptionParser("%s [options]" % os.path.basename(sys.argv[0]))
if sys.platform != 'win32':
parser.add_option("-d", "--daemonize", action="store_true", dest="daemonize", default=False, help="run in daemon mode")
parser.add_option("-e", "--loglevel", dest="loglevel", default="WARNING", help="set the log level")
parser.add_option("-f", "--file", dest="primedb", default="primaton.db", help="file to store detected primes in")
parser.add_option("-l", "--log", dest="logfile", default="primaton.log", help="log file for primaton activity log")
parser.add_option("-o", "--out", action="store_true", dest="output", default=False, help="dump all detected primes from database file to stdout")
parser.add_option("-p", "--pid", dest="pidfile", default="/var/run/primaton.pid", help="PID file when running in daemon mode")
parser.add_option("-u", "--ulimit", dest="ulimit", default=0, help="upper limit to check for prime numbers")
parser.add_option("-w", "--worker", dest="worker", default=1, help="number of worker threads to create")
(options, args) = parser.parse_args()
conf = PTConf()
conf.set("primedb", os.path.normcase(os.path.abspath(options.primedb)))
conf.set("pidfile", os.path.normcase(os.path.abspath(options.pidfile)))
conf.set("logfile", os.path.normcase(os.path.abspath(options.logfile)))
conf.set("loglevel", options.loglevel)
conf.set("logformat", "%(asctime)s [%(levelname)-8s] %(message)s")
conf.set("logdatefmt", "%Y-%m-%d %H:%M:%S")
# This is only for the case a programmer does something stupid
# like using print() in daemon mode...
conf.set("dumpfile", os.path.normcase(os.path.abspath("primaton.out")))
# initialise primaton state
PrimatonState = PTStatus()
# is it only required to dump the prime status to standard output?
if options.output:
dump()
sys.exit(0)
# switch into daemon mode if desired
if sys.platform != 'win32':
if options.daemonize:
daemonize()
FeedQueue = Queue.Queue()
# create worker threads (daemon mode)
worker_threads = [PrimeWorker(FeedQueue, PrimatonState) for i in range(int(options.worker))]
for thread in worker_threads:
thread.setDaemon(True)
thread.start()
# start processing of primes
lg = PTLogger()
lg.log('INFO', "Starting prime research at %d" % PrimatonState.latest_prime())
while long(options.ulimit) == 0 or PrimatonState.latest_prime() < long(options.ulimit):
# log what we're doing
lg.log('INFO', "Processing interval from [p, q] = %s" % PrimatonState.get_pq())
# put current square chunk into the processing queue
for i in PrimatonState.candidate_range():
FeedQueue.put(i)
# wait until all queued items have been processed
FeedQueue.join()
# Forward to next square chunk and save state
SuperLock.acquire()
PrimatonState.forward()
PrimatonState.save()
SuperLock.release()
if __name__ == '__main__': main()
No comments | Defined tags for this entry: algorithms, programming, Project Euler, python
Comments
No comments
Content is subject to the conditions of the