Evolution of an implementation of the Sieve of Eratosthenes
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)) |
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; | |
} | |
} | |
} | |
} |
// 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 | |
} |
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)) |
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)) |
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; | |
} | |
} | |
} | |
} |
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; | |
} | |
} | |
} | |
} |
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