Create a gist now

Instantly share code, notes, and snippets.

@EricBurnett /array1.py
Created Mar 3, 2010

What would you like to do?
Download ZIP
Evolution of an implementation of the Sieve of Eratosthenes
Raw
array1.py
import math
from time import clock
# Calculates all the primes from 0 to stop_at using a sieve on an array.
def primes(stop_at):
if stop_at is None:
print("This algorithm doesn't support unbounded ranges")
return []
if stop_at <= 2: return []
# Candidates from 3 to stop_at-1, except only odd values.
length = (stop_at-1) // 2
s = [True] * length
root = int(math.sqrt(stop_at))
i = 0
p = 3
while p <= root:
if s[i]: # is p prime?
j = (p*p-3) // 2 # index of (p^2)
while j < length:
s[j]=False
j += p
i = i + 1
p = 2*i + 3
return [2]+[x*2+3 for x in range(length) if s[x]]
# Runs the function described by the string s, then prints out how long
# it took to run.
def time(s):
t = clock()
eval(s)
print(s + " took " + str(clock() - t))
Raw
array2.cs
using System;
using System.Collections;
using System.Collections.Generic;
using System.Diagnostics;
namespace PrimeGen
{
using NumType = Int64;
class PrimeGen {
static void Main(string[] args) {
NumType limit = 0;
if (args.Length < 1 || !NumType.TryParse(args[0], out limit)) {
Console.Error.WriteLine(
"Usage: " + System.AppDomain.CurrentDomain.FriendlyName +
" {limit}");
return;
}
Console.Out.WriteLine("Using limit of " + limit + ":");
Stopwatch sw = new Stopwatch();
sw.Start();
NumType sum = 0;
foreach (var i in GeneratePrimesArray(limit)) {
sum += i;
}
sw.Stop();
Console.Out.WriteLine("Array: ");
Console.Out.WriteLine("\tSum: " + sum);
Console.Out.WriteLine("\tTime elapsed: " + sw.Elapsed.ToString());
}
public static IEnumerable<NumType> GeneratePrimesArray(NumType limit) {
if (limit <= 2) yield break;
NumType length = (limit - 1) / 2;
BitArray s = new BitArray((int)length, true);
int root = (int)Math.Sqrt(limit);
int i = 0;
int p = 3;
yield return 2;
while (p <= root) {
if (s.Get(i)) {
yield return p;
int j = (p * p - 3) / 2;
while (j < length) {
s.Set(j, false);
j += p;
}
}
i += 1;
p = 2 * i + 3;
}
while (i < length) {
if (s.Get(i)) {
yield return p;
}
i += 1;
p = 2 * i + 3;
}
}
}
}
Raw
segmented_primes.go
// Go implementation of a segmented Sieve of Eratosthenes. Much faster than the lazy
// evaluation variants above. Sieves first 10^12 primes in 1 hour (4 cores); 10^9
// candidates in the region of 10^18 takes about 2 minutes.
package main
import "flag"
import "fmt"
import "log"
import "math"
import "math/big"
import "time"
import "runtime"
var min_prime = flag.Int64("min", 0,
"Generate primes starting at this number.")
var max_prime = flag.Int64("max", 1000000,
"Generate primes up to and including this number.")
var processes = flag.Int("processes", 2, "Number of processes to use. "+
"Best depends on the processor itself and the other flags set.")
var sum_only = flag.Bool("sum_only", false, "Output only the sum to the main "+
"thread, not every individual prime.")
var vanilla = flag.Bool("vanilla", false,
"Skip the segmenting entirely and use one large array.")
// Fine tuning nobs; probably don't need to change these.
var segment_size = flag.Int64("segment_size", 0,
"Number of values in each segment to check.")
var segment_buffer = flag.Int("buffer", 0,
"Optional override for the size of channel buffers to use.")
var DEBUG = false
func main() {
flag.Parse()
if *max_prime < *min_prime {
log.Fatal("Max must be >= min.")
}
segment_min := int64(math.Sqrt(float64(*max_prime)))
if *segment_size == 0 {
*segment_size = segment_min * 2
} else if segment_min > *segment_size {
fmt.Printf("WARNING: Low segment size detected. Recommend setting "+
"--segment_size to at least %v (currently %v)\n",
segment_min, *segment_size)
}
if *sum_only {
size := big.NewInt(*segment_size)
max := big.NewInt(*max_prime)
size.Mul(size, max)
if size.Cmp(big.NewInt((1<<62)-1)) > 0 {
log.Fatalln("ERROR: Numbers too large for --sum_only to be safe.")
}
}
if *segment_buffer == 0 {
*segment_buffer = int(*segment_size / 20)
}
runtime.GOMAXPROCS(*processes)
t := time.Now()
Run(*min_prime, *max_prime)
fmt.Println("Took", time.Since(t))
}
func Run(min int64, max int64) {
// Schedule worker that fills a channel with segments
// Each segment does its work then fills a channel with the results
// Main thread then consumes segments; consumes primes from segment
primes := make(chan int64, *segment_buffer**processes)
go GenPrimes(min, max, primes)
sum := new(big.Int)
accumulator := int64(0)
for v := range primes {
if v > (1<<62)-1 {
AddInt(sum, v)
continue
}
if accumulator > (1<<62)-1 {
AddInt(sum, accumulator)
accumulator = 0
}
accumulator += v
}
AddInt(sum, accumulator)
fmt.Printf("Sum of primes from %v up to %v is %v\n", min, max, sum)
}
// Returns a slice of all primes in [3, max]. Uses a linear array, no segmenting.
func ArrayPrimes(max int64) []int64 {
if Even(max) {
max -= 1
}
primes := make([]int64, 0)
bits := NewBitArray(3, max)
for i := int64(3); i <= max; i += 2 {
if bits.Bit(i) == false {
// Prime!
primes = append(primes, i)
for next := i * i; next <= max; next += i * 2 {
bits.SetBit(next)
}
}
}
return primes
}
func GenPrimes(min int64, max int64, out chan<- int64) {
// Pre-calculate first sqrt(n) primes, for all the segments to reference.
var pre int64
if *vanilla {
pre = max
} else {
pre = int64(math.Sqrt(float64(max))) + 1
}
pre_primes := ArrayPrimes(pre)
if min <= 2 {
out <- 2
}
for _, v := range pre_primes {
if v >= min {
out <- v
}
}
segments := make(chan chan int64, *processes-1)
// Copy primes from the individual segments to 'out', closing it when every
// prime is in.
go func() {
start_time := time.Now()
last_time := start_time
num_segments := (max - min) / *segment_size
if num_segments == 0 {
num_segments = 1
}
done_segments := int64(0)
for s := range segments {
for v := range s {
out <- v
}
done_segments += 1
if time.Since(last_time) > time.Second*5 {
fmt.Printf("Elapsed %v Done segments %v/%v (%v%%)\n", time.Since(start_time),
done_segments, num_segments, 100.0*float64(done_segments)/float64(num_segments))
last_time = time.Now()
}
}
close(out)
}()
// Now setup the segments themselves for calculating the batches.
// Each outputs to its own channel, to keep ordering intact, although they may
// be processed in parallel.
next := pre + 1
if min > next {
next = min
}
for next <= max {
next_end := next + *segment_size
if max < next_end {
next_end = max
}
var c chan int64
if *sum_only {
c = make(chan int64, 1)
} else {
c = make(chan int64, *segment_buffer)
}
go GenSegment(next, next_end, pre_primes, c)
segments <- c
next = next_end + 1
}
close(segments)
}
// Generates primes in [min, max] and outputs them to out.
// 'primes' must contain all primes <= sqrt(max) (except 2).
func GenSegment(min int64, max int64, primes []int64, out chan<- int64) {
if Even(min) {
min += 1
}
if Even(max) {
max -= 1
}
bits := NewBitArray(min, max)
for _, p := range primes {
start := p * p
if start > max {
break
}
// Find first value that is a multiple of p, >= min, and odd.
if start < min {
start = min - (min % p)
if start < min {
start += p
}
if Even(start) {
start += p
}
}
for start <= max {
bits.SetBit(start)
start += 2 * p
}
}
// Output the primes remaining.
sum := int64(0)
sum_only := *sum_only
for i := min; i <= max; i += 2 {
if bits.Bit(i) == false {
sum += i
if !sum_only {
out <- i
}
}
}
if sum_only {
out <- sum
}
close(out)
}
func Even(num int64) bool {
return num&1 == 0
}
func AddInt(z *big.Int, x int64) {
x_big := big.NewInt(x)
z.Add(z, x_big)
}
// A structure for efficiently addressing odd bits. Similar performance to a
// bool array, except requires much less space.
type BitArray struct {
Min int64
Max int64
data []byte
}
func NewBitArray(min int64, max int64) *BitArray {
if Even(min) || Even(max) {
log.Fatalf("min and max must be odd; got (%v, %v)\n", min, max)
}
elements := ((max - min) / 2) + 1
data := make([]byte, (elements+7)>>3)
return &BitArray{min, max, data}
}
func (b *BitArray) Bit(index int64) bool {
if DEBUG && (index < b.Min || index > b.Max) {
log.Fatalf("Index out of bounds! %v outside [%v, %v]\n", index, b.Min, b.Max)
}
offset := (index - b.Min) / 2
offset_byte := offset >> 3
offset_bit := uint(offset & 0x7)
return b.data[offset_byte]&(1<<offset_bit) > 0
}
func (b *BitArray) SetBit(index int64) {
if DEBUG && (index < b.Min || index > b.Max) {
log.Fatalf("Index out of bounds! %v outside [%v, %v]\n", index, b.Min, b.Max)
}
offset := (index - b.Min) / 2
offset_byte := offset >> 3
offset_bit := uint(offset & 0x7)
b.data[offset_byte] |= 1 << offset_bit
}
Raw
v1.py
import heapq
import itertools
from time import clock
# Returns an iterator for generating all primes < stop
def primes(stop=None):
class primeIterObject:
def __init__(self, p, it):
self.car = p*p
self.cdr = map(lambda x:x*p, it)
def next(self):
self.car = next(self.cdr)
def __eq__(self, other):
return self.car == other.car
def __ne__(self, other):
return self.car != other.car
def __lt__(self, other):
return self.car < other.car
def __gt__(self, other):
return self.car > other.car
def __le__(self, other):
return self.car <= other.car
def __ge__(self, other):
return self.car >= other.car
# Find all primes by iterating though the candidates and checking
# if each is a known composite (multiple of some prime). If not, build
# a composite iterator based on it and add it to the heap of iterators
# to check against.
# Based on http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf
heap = [] # Set of composite iterators, one for each prime
def incHeap(heap):
front = heap[0]
front.next()
heapq.heapreplace(heap, front)
def betterIterator(n):
cur = n
cycle = itertools.cycle([2, 4])
if (n+1) % 3 != 0:
next(cycle)
while True:
yield cur
cur += next(cycle)
itmaker = betterIterator
it = itmaker(5)
cur = next(it)
yield 2
heap.append(primeIterObject(2, itmaker(5)))
yield 3
heap.append(primeIterObject(3, itmaker(5)))
while stop is None or cur < stop:
while heap[0].car < cur:
incHeap(heap)
if heap[0].car == cur:
while heap[0].car == cur:
incHeap(heap)
else:
yield cur
if stop is None or cur*cur <= stop:
it2 = itmaker(cur)
next(it2)
heapq.heappush(heap, primeIterObject(cur, it2))
cur = next(it)
# Runs the function described by the string s, then prints out how long
# it took to run.
def time(s):
t = clock()
eval(s)
print(s + " took " + str(clock() - t))
Raw
v2.py
import heapq
import itertools
from time import clock
# Returns an iterator for generating all primes < stop.
def primes(stop):
if stop is None:
print("This algorithm doesn't support unbounded ranges")
raise StopIteration
# Find all primes by iterating though the candidates and checking
# if each is a known composite (multiple of some prime). If not, build
# a composite iterator based on it and add it to the heap of iterators
# to check against.
# Based loosely on http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf.
md = {}
def addToMultiDict(d, key, value):
if key in d:
d[key].append(value)
else:
d[key] = [value]
# This iterator skips all multiples of 2 and 3
# Can apply an optional multiplier
def betterIterator(n, mult = 1):
cycle = itertools.cycle([2, 4])
if (n+1) % 3 != 0:
next(cycle)
while True:
yield n * mult
n += next(cycle)
itmaker = betterIterator
it = itmaker(5)
cur = next(it)
yield 2
yield 3
while cur < stop:
if cur in md:
allIters = md[cur]
del md[cur]
for i in allIters:
n = next(i)
addToMultiDict(md, n, i)
else:
yield cur
if stop is None or cur*cur < stop:
it2 = itmaker(cur, cur)
n = next(it2)
addToMultiDict(md, n, it2)
cur = next(it)
# Runs the function described by the string s, then prints out how long
# it took to run.
def time(s):
t = clock()
eval(s)
print(s + " took " + str(clock() - t))
Raw
v3.cs
using System;
using System.Collections.Generic;
using System.Diagnostics;
namespace PrimeGen
{
using NumType = Int64;
public class MultiDictionary<TKey, TVal> : Dictionary<TKey, List<TVal>> {
public void Add(TKey key, TVal val) {
List<TVal> container;
if (!this.TryGetValue(key, out container)) {
container = new List<TVal>();
base.Add(key, container);
}
container.Add(val);
}
}
class PrimeGen
{
static void Main(string[] args) {
NumType limit = 0;
if (args.Length < 1 || !NumType.TryParse(args[0], out limit)) {
Console.Error.WriteLine(
"Usage: " + System.AppDomain.CurrentDomain.FriendlyName +
" {limit}");
return;
}
Console.Out.WriteLine("Using limit of " + limit + ":");
Stopwatch sw = new Stopwatch();
sw.Start();
NumType sum = 0;
foreach (var i in GeneratePrimes(limit)) {
sum += i;
}
sw.Stop();
Console.Out.WriteLine("Sum: " + sum);
Console.Out.WriteLine("Time elapsed: " + sw.Elapsed.ToString());
}
// Returns candidate composites starting at first. If 'multiplier' is
// specified, all returned candidates are multiplied by 'multiplier'
// before being returned.
public static IEnumerator<NumType> Candidates(NumType first,
NumType multiplier = 1) {
NumType[] steps = {2, 4};
int nextStepIndex = 0;
if ((first + 1) % 3 != 0) {
nextStepIndex = 1;
}
while (true) {
yield return first * multiplier;
first += steps[nextStepIndex];
nextStepIndex = 1 - nextStepIndex;
}
}
// Returns a generator for prime numbers < limit
public static IEnumerable<NumType> GeneratePrimes(NumType limit) {
if (limit <= 2) yield break;
yield return 2;
if (limit == 3) yield break;
yield return 3;
var generators = new MultiDictionary<NumType, IEnumerator<NumType>>();
var it = Candidates(5);
it.MoveNext();
NumType cur = it.Current;
while (cur < limit) {
if (generators.ContainsKey(cur)) {
var all = generators[cur];
generators.Remove(cur);
foreach (var i in all) {
i.MoveNext();
var n = i.Current;
generators.Add(n, i);
}
} else {
yield return cur;
if (cur * cur < limit) {
var it2 = Candidates(cur, cur);
it2.MoveNext();
var n = it2.Current;
generators.Add(n, it2);
}
}
it.MoveNext();
cur = it.Current;
}
}
}
}
Raw
v4.cs
using System;
using System.Collections.Generic;
using System.Diagnostics;
namespace PrimeGen
{
using NumType = Int64;
public class MultiDictionary<TKey, TVal> : Dictionary<TKey, SimpleContainer<TVal>> {
public void Add(TKey key, TVal val) {
SimpleContainer<TVal> container;
if (!this.TryGetValue(key, out container)) {
container = new SimpleContainer<TVal>();
base.Add(key, container);
}
container.Add(val);
}
}
public class SimpleContainer<T> {
T[] array;
int count;
public void Clear() {
count = 0;
}
public void Add(T value) {
if (array == null) {
array = new T[6];
}
if (count >= array.Length) {
Array.Resize(ref array, array.Length * 2);
}
array[count] = value;
count++;
}
public int Count {
get {
return count;
}
}
public T this[int index] {
get {
if (index < count) {
return array[index];
} else {
throw new ArgumentOutOfRangeException("index");
}
}
}
}
class PrimeGen
{
static void Main(string[] args) {
NumType limit = 0;
if (args.Length < 1 || !NumType.TryParse(args[0], out limit)) {
Console.Error.WriteLine(
"Usage: " + System.AppDomain.CurrentDomain.FriendlyName +
" {limit}");
return;
}
Console.Out.WriteLine("Using limit of " + limit + ":");
Stopwatch sw = new Stopwatch();
sw.Start();
NumType sum = 0;
foreach (var i in GeneratePrimes(limit)) {
sum += i;
}
sw.Stop();
Console.Out.WriteLine("Sum: " + sum);
Console.Out.WriteLine("Time elapsed: " + sw.Elapsed.ToString());
}
// Returns candidate composites starting at first. If 'multiplier' is
// specified, all returned candidates are multiplied by 'multiplier'
// before being returned.
public static IEnumerator<NumType> Candidates(NumType first,
NumType multiplier = 1) {
NumType[] steps = {2, 4};
int nextStepIndex = 0;
if ((first + 1) % 3 != 0) {
nextStepIndex = 1;
}
while (true) {
yield return first * multiplier;
first += steps[nextStepIndex];
nextStepIndex = 1 - nextStepIndex;
}
}
// Returns a generator for prime numbers < limit
public static IEnumerable<NumType> GeneratePrimes(NumType limit) {
if (limit <= 2) yield break;
yield return 2;
if (limit == 3) yield break;
yield return 3;
var generators = new MultiDictionary<NumType, IEnumerator<NumType>>();
var it = Candidates(5);
it.MoveNext();
NumType cur = it.Current;
while (cur < limit) {
if (generators.ContainsKey(cur)) {
var all = generators[cur];
generators.Remove(cur);
for (int i = 0; i < all.Count; ++i) {
all[i].MoveNext();
var n = all[i].Current;
generators.Add(n, all[i]);
}
} else {
yield return cur;
if (cur * cur < limit) {
var it2 = Candidates(cur, cur);
it2.MoveNext();
var n = it2.Current;
generators.Add(n, it2);
}
}
it.MoveNext();
cur = it.Current;
}
}
}
}
Raw
v5.cs
using System;
using System.Collections.Generic;
using System.Diagnostics;
namespace PrimeGen
{
using NumType = Int64;
public struct SimpleContainer<T> {
T[] array;
int count;
public void Clear() {
count = 0;
}
public void Add(T value) {
if (array == null) {
array = new T[6];
}
if (count >= array.Length) {
Array.Resize(ref array, array.Length * 2);
}
array[count] = value;
count++;
}
public int Count {
get {
return count;
}
}
public T this[int index] {
get {
if (index < count) {
return array[index];
} else {
throw new ArgumentOutOfRangeException("index");
}
}
}
}
// This class represents a list of elements [firstIndex, firstIndex + size).
// When Increment is called, the available window shifts forward by one.
// Each element in the set can hold multiple values, and every element is
// initially guaranteed to be empty.
// It is tailored to only allow values not divisible by 2 or 3
public class MultiCyclicArray<T> {
private SimpleContainer<T>[] array;
private NumType firstIndex;
private NumType firstIndexCount;
private int offset;
public MultiCyclicArray(int arraySize, NumType firstIndex) {
// Approximate the number of elements in a window not
// divisible by 2 or 3.
arraySize = arraySize / 3 + 2;
this.array = new SimpleContainer<T>[arraySize];
for (int i = 0; i < arraySize; ++i) {
array[i] = new SimpleContainer<T>();
}
this.firstIndex = firstIndex;
firstIndexCount = firstIndex / 3 + 1;
this.offset = 0;
}
public SimpleContainer<T> this[NumType index] {
get {
return Get(index);
}
}
public SimpleContainer<T> Get(NumType index) {
return array[TranslateIndex(index)];
}
public void Add(NumType index, T value) {
array[TranslateIndex(index)].Add(value);
}
public void Clear(NumType index) {
array[TranslateIndex(index)].Clear();
}
// Takes an index from the user's perspective, and returns the raw array
// index it represents.
private int TranslateIndex(NumType emulatedIndex) {
Debug.Assert(emulatedIndex % 6 == 1 || emulatedIndex % 6 == 5);
NumType emulatedIndexCount = emulatedIndex / 3 + 1;
if (emulatedIndexCount < firstIndexCount ||
emulatedIndexCount >= firstIndexCount + array.Length) {
throw new ArgumentException("index outside of current window!");
}
return GetRawIndex((int) (emulatedIndexCount - firstIndexCount));
}
// Takes an index into the current window, and applies the offsetting
// to get an index into the base array.
private int GetRawIndex(int unoffsettedIndex) {
var index = unoffsettedIndex + offset;
if (index >= array.Length) {
index -= array.Length;
}
return index;
}
// Shifts the window right by num elements.
public void Increment(int num) {
if (num != 2 && num != 4) {
throw new ArgumentException("Only supports incrementing to next number not divisible by 2 or 3!",
"num");
}
array[offset].Clear();
firstIndex += num;
firstIndexCount += 1;
offset = offset + 1;
if (offset >= array.Length) {
offset -= array.Length;
}
}
}
class PrimeGen
{
static void Main(string[] args) {
NumType limit = 0;
if (args.Length < 1 || !NumType.TryParse(args[0], out limit)) {
Console.Error.WriteLine(
"Usage: " + System.AppDomain.CurrentDomain.FriendlyName +
" {limit}");
return;
}
Console.Out.WriteLine("Using limit of " + limit + ":");
Stopwatch sw = new Stopwatch();
sw.Start();
NumType sum = 0;
foreach (var i in GeneratePrimes(limit)) {
sum += i;
}
sw.Stop();
Console.Out.WriteLine("Sum: " + sum);
Console.Out.WriteLine("Time elapsed: " + sw.Elapsed.ToString());
}
// Returns candidate composites starting at first. If 'multiplier' is
// specified, all returned candidates are multiplied by 'multiplier'
// before being returned.
public static int MaxStep = 4;
public static IEnumerator<NumType> Candidates(NumType first,
NumType multiplier = 1) {
NumType[] steps = {2, 4};
int nextStepIndex = 0;
if ((first + 1) % 3 != 0) {
nextStepIndex = 1;
}
while (true) {
yield return first * multiplier;
first += steps[nextStepIndex];
nextStepIndex = 1 - nextStepIndex;
}
}
// Returns a generator for prime numbers < limit
public static IEnumerable<NumType> GeneratePrimes(NumType limit) {
if (limit <= 2) yield break;
yield return 2;
if (limit == 3) yield break;
yield return 3;
// Needs to support values up to the maximum prime (<= sqrt(limit)),
// which can be incremented MaxStep times to the next potential
// composite.
var generators = new MultiCyclicArray<IEnumerator<NumType>>(((int) Math.Sqrt(limit))*MaxStep + 1, 5);
var queuedGenerators = new Queue<IEnumerator<NumType>>();
var it = Candidates(5);
it.MoveNext();
NumType cur = it.Current;
while (cur < limit) {
var all = generators[cur];
var count = all.Count;
if (count > 0) {
for (int i = 0; i < count; ++i) {
var it2 = all[i];
it2.MoveNext();
var n = it2.Current;
generators.Add(n, it2);
}
} else {
yield return cur;
if (cur * cur < limit) {
var it2 = Candidates(cur, cur);
it2.MoveNext();
queuedGenerators.Enqueue(it2);
}
}
NumType old = cur;
it.MoveNext();
cur = it.Current;
generators.Increment((int) (cur - old));
if (queuedGenerators.Count > 0 &&
queuedGenerators.Peek().Current == cur) {
var it2 = queuedGenerators.Dequeue();
generators.Add(it2.Current, it2);
}
}
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment