Sieve of Atkin

An implementation of Sieve of Atkin algorithm for prime number generation up to 100,000,000.

This is a parallel Java version of Sieve of Atkin algorithm optimized for multi-core platforms. Fork-Join framework (see java.util.concurrent) is applied to make full usage of CPU resources. Also, in order to accelerate disk IO, FileChannel and MappedByteBuffer are used to replace old IO classes. Even with BufferedOutputStream, MappedByteBuffer is still 25% faster.

The results are saved to a compact binary file filled with 32-bits integers. The file contains the first 5,761,455 prime numbers. The last one is 99,999,989.

The program totally ran in 1.8 sec including 1.1 sec for prime generation and 0.69 sec for file operations.

Source code's as follows:

import java.util.concurrent.*;
import java.io.*;
import java.nio.*;
import java.nio.channels.*;
import java.util.*;
import java.text.*;

public class SieveAtkin extends RecursiveAction {
    private int from, to;
    private int len;
    private boolean[] isPrime;
    private int fromGlobal;

    public SieveAtkin (int from, int to, boolean[] isPrime, int fromGlobal) {
        this.from = from;
        this.to = to;
        len = to - from + 1;
        this.isPrime = isPrime;
        this.fromGlobal = fromGlobal;
    }

    private void genPrimes () {
        for (int yy=1, i=1, n=0; yy<=to-4; yy+=(i+=2)) {
            int f = from - yy;
            int low = 4;
            int rlow = 1;
            if (f > 4) {
                int r = (int) Math.sqrt(f/4.0);
                if (r * r * 4 == f) {
                    low = f;
                } else {
                    r++;
                    low = r * r * 4;
                }
                rlow = r;
            }
            for (int xx4=low, j=8*rlow-4; (n=xx4+yy)<=to; xx4+=(j+=8)){
                int m = n % 12;
                if (m == 1 || m == 5) {
                    isPrime[n-fromGlobal] = !isPrime[n-fromGlobal];
                }
            }
        }
        for (int yy=1, i=1, n=0; yy<=to-3; yy+=(i+=2)) {
            int f = from - yy;
            int low = 3;
            int rlow = 1;
            if (f > 3) {
                int r = (int) Math.sqrt(f/3.0);
                if (r * r * 3 == f) {
                    low = f;
                } else {
                    r++;
                    low = r * r * 3;
                }
                rlow = r;
            }
            for (int xx3=low, j=6*rlow-3; (n=xx3+yy)<=to; xx3+=(j+=6)){
                int m = n % 12;
                if (m == 7) {
                    isPrime[n-fromGlobal] = !isPrime[n-fromGlobal];
                }
            }
        }
        for (int y=1, yy=1, i=1, t=0, n=0; (t=2*yy+6*y+3)<=to; y++, yy+=(i+=2)) {
            int f = from + yy;
            int low = (yy + i + 2) * 3;
            int rlow = y + 1;
            if (f > low) {
                int r = (int) Math.sqrt(f/3.0);
                if (r * r * 3 == f) {
                    low = f;
                } else {
                    r++;
                    low = r * r * 3;
                }
                rlow = r;
            }
            for (int xx3=low, j=6*rlow-3; (n=xx3-yy)<=to; xx3+=(j+=6)){
                int m = n % 12;
                if (m == 11) {
                    isPrime[n-fromGlobal] = !isPrime[n-fromGlobal];
                }
            }
        }
    }

    protected void compute () {
        if (len <= 10000000) {
            genPrimes();
            return;
        }
        invokeAll(
            new SieveAtkin(from, from + len / 2, isPrime, fromGlobal),
            new SieveAtkin(from + len / 2 + 1, to, isPrime, fromGlobal)
        );
    }

    public static void main (String args[]) throws IOException {
        int to = 100000000;
        int from = 2;
        int len = to - from + 1;

        boolean[] isPrime = new boolean[len];
        isPrime[0] = true;
        isPrime[1] = true;

        new ForkJoinPool().invoke(new SieveAtkin(from, to, isPrime, from));

        for (int i=5-from; i<len; i++) {
            if (isPrime[i]) {
                int p = i + from;
                int pp = p * p;
                if (pp > to) {
                    break;
                }
                for (int q=pp; q<=to; q+=pp) {
                    isPrime[q - from] = false;
                }
            }
        }

        int size = 0;
        for (int i=0; i<isPrime.length; i++) {
            if (isPrime[i]) {
                size += 8;
            }
        }
        RandomAccessFile raf = new RandomAccessFile("prime", "rw");
        FileChannel ch = raf.getChannel();
        MappedByteBuffer buff = ch.map(FileChannel.MapMode.READ_WRITE, 0, size);
        for (int i=0; i<isPrime.length; i++) {
            if (isPrime[i]) {
                buff.putInt(from+i);
            }
        }
        ch.close();
    }
}