{{ message }}

Instantly share code, notes, and snippets.

# EricBurnett/array1.py

Created Mar 3, 2010
Evolution of an implementation of the Sieve of Eratosthenes
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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))
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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 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; } } } }
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 // 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< 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 }
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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))
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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))
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 using System; using System.Collections.Generic; using System.Diagnostics; namespace PrimeGen { using NumType = Int64; public class MultiDictionary : Dictionary> { public void Add(TKey key, TVal val) { List container; if (!this.TryGetValue(key, out container)) { container = new List(); 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 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 GeneratePrimes(NumType limit) { if (limit <= 2) yield break; yield return 2; if (limit == 3) yield break; yield return 3; var generators = new MultiDictionary>(); 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; } } } }
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 using System; using System.Collections.Generic; using System.Diagnostics; namespace PrimeGen { using NumType = Int64; public class MultiDictionary : Dictionary> { public void Add(TKey key, TVal val) { SimpleContainer container; if (!this.TryGetValue(key, out container)) { container = new SimpleContainer(); base.Add(key, container); } container.Add(val); } } public class SimpleContainer { 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 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 GeneratePrimes(NumType limit) { if (limit <= 2) yield break; yield return 2; if (limit == 3) yield break; yield return 3; var generators = new MultiDictionary>(); 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; } } } }
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 using System; using System.Collections.Generic; using System.Diagnostics; namespace PrimeGen { using NumType = Int64; public struct SimpleContainer { 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 { private SimpleContainer[] 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[arraySize]; for (int i = 0; i < arraySize; ++i) { array[i] = new SimpleContainer(); } this.firstIndex = firstIndex; firstIndexCount = firstIndex / 3 + 1; this.offset = 0; } public SimpleContainer this[NumType index] { get { return Get(index); } } public SimpleContainer 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 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 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>(((int) Math.Sqrt(limit))*MaxStep + 1, 5); var queuedGenerators = new Queue>(); 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); } } } } }