Skip to content

Instantly share code, notes, and snippets.

@EricBurnett EricBurnett/array1.py
Created Mar 3, 2010

Embed
What would you like to do?
Learn more about clone URLs
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
You can’t perform that action at this time.