Sieb von Atkin – Nachlese


Hier eine etwas performantere Version des Siebs von Atkin:

object atkin {

    class IsPrime(limit:Int) extends
                       scala.collection.Traversable[Boolean] {
        import scala.annotation.switch

        private val array1 = Array.fill(1 + limit/6)(false)
        private val array5 = Array.fill(1 + limit/6)(false)

        def apply(index:Int) =
        if (index == 2 || index == 3 || index == 5) true
        else (index % 6: @switch) match {
            case 1 => array1(index/6)
            case 5 => array5(index/6)
            case _ => false
        }

        def update(index:Int, value:Boolean) {
            (index % 6: @switch) match {
                case 1 => array1(index/6) = value
                case 5 => array5(index/6) = value
                case _ => if (value) error("Illegal Index")
            }
        }

        def foreach[U](f: Boolean => U) {
            (0 until limit).foreach(index => f(apply(index)))
        }
        def zipWithIndex =
           (0 until limit).map(index => (apply(index), index))
        def primes = (0 until limit).filter(apply(_)).force
    }
    
    def sieve(limit:Int) = {
        val isPrime = new IsPrime(limit)
        for{x <- 1 to (Math.sqrt(limit/4).toInt)+1
            y <- 1 to Math.sqrt(limit - 4*x*x).toInt by 2
            if (x % 3 + y % 3 != 0)} {
            val n = 4*x*x+y*y
            if (n < limit) isPrime(n) = ! isPrime(n)
        }
        for{x <- 1 to Math.sqrt(limit/3).toInt+1 by 2
            y <- 1 to Math.sqrt(limit - 3*x*x).toInt
            if y % 3 != 0 && y % 2 == 0} {
            val n = 3*x*x+y*y
            if (n < limit) isPrime(n) = ! isPrime(n)
        }
        for{x <- 1 to Math.sqrt(limit/2).toInt+1
            y <- 1 until x if y % 3 != 0 && x % 2 + y % 2 == 1} {
            val n = 3*x*x-y*y
            if (n < limit) isPrime(n) = ! isPrime(n)
        }
        for{n <- 5 until Math.sqrt(limit).toInt+1 by 2 if isPrime(n)
            s <- n*n until limit by n*n}
               isPrime(s) = false
        isPrime.primes
    }
}

Nicht, dass ich damit irgendwelche Geschwindigkeitsrekorde brechen wollte – ich habe auch keine Performance-Messungen vorgenommen. Trotzdem ist es interessant, diese Version mit meiner naiven Variante zu vergleichen – dabei kann man sich vielleicht den einen oder anderen Trick abschauen.

Zuerst einmal war das alte Boolean-Array eine enorme Speicherplatzverschwendung, denn das Sieb arbeitet ja nur mit Zahlen 1 und 5 (mod 6), d.h. 2/3 des Arrays werden gar nicht genutzt – es steht von vornherein fest, dass diese Zahlen keine Primzahlen sind. Also habe ich mir eine Wrapper-Klasse geschrieben, die intern zwei Boolean-Arrays für die Zahlen 1 und 5 (mod 6) enthält und sich nach außen „in etwa“ wie ein normales Array verhält. Die „magische“ apply-Methode habe ich inzwischen oft genug erwähnt, sie hat aber auch ein „schreibendes“ Gegenstück namens update. In Java müssen wir bei einer Liste z.B. list.set(10, „X“) schrieben, anstatt ähnlich wie bei Arrays einfach list(10) = „X“. Die update-Methode erlaubt es in Scala, diese Schreibweise zu verwenden.

Auffällig ist die Annotation @switch (in Scala 2.8). Ähnlich wie @tailrec einen Fehler ausgibt, wenn der Compiler wider Erwarten keine Endrekursions-Optimierung durchführen kann, sorgt @switch dafür, dass man gewarnt wird, wenn der Compiler die match-Anweisung nicht in ein performantes switch-Statement, sondern in eine „teurere“ if-Kaskade übersetzt. Natürlich kann das nur klappen, wenn man nur einfache Ausdrücke in den Fall-Klauseln hat – mit Fallklassen oder Guards ist diese Optimierung nicht möglich.

Auch wenn diese Funktionalität in diesem Programm gar nicht verwendet wird, habe ich das Trait Traversable (in Scala 2.8) implementiert. Es verlangt nur eine foreach-Methode, aber dafür bekommt man eine Menge anderer Methoden „umsonst“. Man kann beispielsweise ein IsPrime zu anderen Collections wie Set mit ++ hinzufügen u.s.w.

Die Sieb-Methode selbst ist nicht viel länger geworden. Die erste Schleife mit den drei Test ist in drei einzelne Schleifen aufgeteilt worden, wobei Variable jeweils nur soweit gezählt werden, dass die (n < limit)-Bedinung nicht verletzt wird. Mit Guards und dem "by" für Schrittweiten wird erreicht, dass nur x-y-Paare erzeugt werden, die auch die Test-Bedingungen erfüllen (womit auch die alte invert-Methode überflüssig wird). Um das richtig hinzubekommen, habe ich eine Weile mit Open Office Calc (dem Excel-Pendant) herumgespielt, um zu sehen, welche Kombinationen überhaupt möglich sind – dafür reicht mein Zahlentheorie-Wissen gerade noch.

So, viel mehr fällt mir zum Thema Optimierung nicht ein. Für Verbesserungsvorschläge bin ich natürlich immer dankbar.

Advertisements

Kommentar verfassen

Trage deine Daten unten ein oder klicke ein Icon um dich einzuloggen:

WordPress.com-Logo

Du kommentierst mit Deinem WordPress.com-Konto. Abmelden / Ändern )

Twitter-Bild

Du kommentierst mit Deinem Twitter-Konto. Abmelden / Ändern )

Facebook-Foto

Du kommentierst mit Deinem Facebook-Konto. Abmelden / Ändern )

Google+ Foto

Du kommentierst mit Deinem Google+-Konto. Abmelden / Ändern )

Verbinde mit %s