PrimeSieve
1:
2: namespace Luschny.Math.Primes
3: {
4: using System;
5: using System.Collections;
6: using System.Collections.Generic;
7: using System.IO;
8: using System.Threading;
9: using Luschny.Math.MathUtils;
10:
11: using Xint = System.Numerics.BigInteger;
12:
13: using Xmath = Luschny.Math.MathUtils.Xmath;
14:
15: /// <summary>
16: /// This class implements a prime number sieve using
17: /// an algorithm given by Eratosthenes (276-194 B.C.)
18: /// Author: Peter Luschny
19: /// Version: 2008年09月12日
20: /// </summary>
21: public class PrimeSieve : IPrimeSieve
22: {
23: readonly int[] primes;
24: PositiveRange sieveRange;
25: PositiveRange primeRange;
26: int numberOfPrimes;
27:
28: /// <summary>
29: /// Initializes a new instance of the PrimeSieve.
30: /// </summary>
31: /// <param name="upTo">
32: /// The upper bound of the sieveRange to be sieved, this means,
33: /// the sieveRange [1,n] is searched for prime numbers.
34: /// </param>
35: public PrimeSieve(int upTo)
36: {
37: primes = new int[GetPiHighBound(upTo)];
38: sieveRange = new PositiveRange(1, upTo);
39:
40: numberOfPrimes = MakePrimeList(upTo);
41: primeRange = new PositiveRange(1, numberOfPrimes);
42: }
43:
44: /// <summary>
45: /// Prime number sieve, Eratosthenes (276-194 b. t.)
46: /// This implementation considers only multiples of primes
47: /// greater than 3, so the smallest value has to be mapped to 5.
48: /// Note: There is no multiplication operation in this function.
49: /// </summary>
50: /// <param name="composite">After execution of the function this
51: /// boolean array includes all composite numbers in [5,n]
52: /// disregarding multiples of 2 and 3.
53: /// </param>
54: private static void SieveOfEratosthenes(bool[] composite)
55: {
56: int d1 = 8;
57: int d2 = 8;
58: int p1 = 3;
59: int p2 = 7;
60: int s = 7;
61: int s2 = 3;
62: int n = 0;
63: int len = composite.Length;
64: bool toggle = false;
65:
66: while (s < len) // -- scan the sieve
67: {
68: if (!composite[n++]) // -- if a prime is found
69: { // -- cancel its multiples
70: int inc = p1 + p2;
71:
72: for (int k = s; k < len; k += inc)
73: {
74: composite[k] = true;
75: }
76:
77: for (int k = s + s2; k < len; k += inc)
78: {
79: composite[k] = true;
80: }
81: }
82:
83: if (toggle = !toggle)
84: {
85: s += d2;
86: d1 += 16;
87: p1 += 2;
88: p2 += 2;
89: s2 = p2;
90: }
91: else
92: {
93: s += d1;
94: d2 += 8;
95: p1 += 2;
96: p2 += 6;
97: s2 = p1;
98: }
99: }
100: }
101:
102: /// <summary>
103: /// Transforms the sieve of Eratosthenes into the
104: /// sequence of prime numbers.
105: /// </summary>
106: /// <param name="n">Upper bound of the sieve.</param>
107: /// <returns>Number of primes found.</returns>
108: private int MakePrimeList(int n)
109: {
110: bool[] composite = new bool[n / 3];
111:
112: SieveOfEratosthenes(composite);
113:
114: int[] primes = this.primes; // -- on stack for eff.
115: bool toggle = false;
116: int p = 5, i = 0, j = 2;
117:
118: primes[0] = 2;
119: primes[1] = 3;
120:
121: while (p <= n)
122: {
123: if (!composite[i++])
124: {
125: primes[j++] = p;
126: }
127: // -- never mind, it's ok.
128: p += (toggle = !toggle) ? 2 : 4;
129: }
130:
131: return j; // number of primes
132: }
133:
134: /// <summary>
135: /// GetPiHighBound estimates the number of primes <= n.
136: /// </summary>
137: /// <param name="n">The upper bound of the intervall under
138: /// consideration.</param>
139: /// <returns>a simple estimate of the number of primes <= n.
140: /// </returns>
141: private static int GetPiHighBound(long n)
142: {
143: if (n < 17) return 6;
144: return (int)System.Math.Floor(((double)n) / (System.Math.Log(n) - 1.5));
145: }
146:
147: /// <summary>
148: /// Get the n-th prime number.
149: /// </summary>
150: /// <param name="n">The index of the prime number.</param>
151: /// <returns>The n-th prime number.</returns>
152: public int GetNthPrime(int n)
153: {
154: // Handle potential under- or overflow
155: primeRange.ContainsOrFail(n);
156: return primes[n - 1];
157: }
158:
159: /// <summary>
160: /// Checks if a given number is prime.
161: /// </summary>
162: /// <param name="cand">The number to be checked.</param>
163: /// <returns>True iff the given number is prime.</returns>
164: public bool IsPrime(int cand)
165: {
166: sieveRange.ContainsOrFail(cand);
167: // The candidate is interpreted as an one point interval!
168: return (new PrimeCollection(this, cand)).IsPrime();
169: }
170:
171: /// <summary>
172: /// The default collection from the full sieve.
173: /// </summary>
174: /// <returns>The prime number collection.</returns>
175: public IPrimeCollection GetPrimeCollection()
176: {
177: return new PrimeCollection(this);
178: }
179:
180: /// <summary>
181: /// Gives the collection of the prime numbers in the given intervall.
182: /// </summary>
183: /// <param name="low">The lower bound of the collection interval.</param>
184: /// <param name="high">The higher bound of the collection interval.</param>
185: /// <returns>The collection of the prime numbers between low and high.</returns>
186: public IPrimeCollection GetPrimeCollection(int low, int high)
187: {
188: return new PrimeCollection(this, new PositiveRange(low, high));
189: }
190:
191: /// <summary>
192: /// Gives the collection of the prime numbers in the given range.
193: /// </summary>
194: /// <param name="range">The range of the collection.</param>
195: /// <returns>The prime number collection.</returns>
196: public IPrimeCollection GetPrimeCollection(PositiveRange range)
197: {
198: return new PrimeCollection(this, range);
199: }
200:
201: /// <summary>
202: /// Gives the Product of the prime numbers in the given sieveRange.
203: /// </summary>
204: /// <param name="low">The lower bound of the collection.</param>
205: /// <param name="high">The upper bound of the collection.</param>
206: /// <returns>The Product of the prime numbers between low and high.
207: /// </returns>
208: public Xint GetPrimorial(int low, int high)
209: {
210: return GetPrimorial(new PositiveRange(low, high));
211: }
212:
213: /// <summary>
214: /// Computes the Product of the prime numbers in the given sieveRange.
215: /// </summary>
216: /// <param name="range">The sieveRange of the enumeration.</param>
217: /// <returns>The Product of the prime numbers in the enumeration.
218: /// </returns>
219: public Xint GetPrimorial(PositiveRange range)
220: {
221: int start, size;
222: var pc = new PrimeCollection(this, range);
223: if (pc.GetSliceParameters(out start, out size))
224: {
225: return Xint.One;
226: }
227:
228: return Xmath.Product(primes, start, size);
229: }
230:
231: ////////////////////// Private Inner Class ///////////////////////
232:
233: /// <summary>
234: /// PrimeCollection
235: /// @author Peter Luschny
236: /// @version 2008年09月12日
237: /// </summary>
238: internal class PrimeCollection : IPrimeCollection
239: {
240: readonly PrimeSieve sieve;
241: readonly PositiveRange enumRange;
242: readonly PositiveRange primeRange;
243: readonly int start, end;
244: readonly bool isPrime;
245: int state, next, current;
246:
247: /// <summary>
248: /// Initializes the prime collection for the given sieve.
249: /// </summary>
250: /// <param name="sieve">The sieve, the prime numbers of which
251: /// are to be represented by a collection.</param>
252: public PrimeCollection(PrimeSieve sieve)
253: {
254: this.sieve = sieve;
255: end = sieve.numberOfPrimes - 1;
256:
257: enumRange = sieve.sieveRange;
258: primeRange = sieve.primeRange;
259: }
260:
261: /// <summary>
262: /// Initializes a collection over a subrange of the range
263: /// of the sieve. An exception is thrown, if the range given
264: /// is not contained in the range of the sieve.
265: /// </summary>
266: /// <param name="sieve">Prime number sieve to be used.</param>
267: /// <param name="enumRange">Range of collection.</param>
268: public PrimeCollection(PrimeSieve sieve, PositiveRange enumRange)
269: {
270: sieve.sieveRange.ContainsOrFail(enumRange);
271:
272: this.sieve = sieve;
273: this.enumRange = enumRange;
274:
275: int nops = sieve.numberOfPrimes;
276: start = IndexOf(enumRange.Min - 1, 0, nops - 1);
277: end = IndexOf(enumRange.Max, start, nops) - 1;
278:
279: if (end < start) //-- PrimeRange is empty.
280: {
281: end = -1;
282: primeRange = new PositiveRange(0, 0);
283: }
284: else
285: {
286: primeRange = new PositiveRange(start + 1, end + 1);
287: }
288: }
289:
290: /// <summary>
291: /// Initializes a collection consisting out of a single value.
292: /// An integer cand is prime iff there is a prime number
293: /// in the range (cand-1, cand].
294: /// </summary>
295: /// <param name="sieve">The sieve to be used.</param>
296: /// <param name="cand">The prime number candidate.</param>
297: public PrimeCollection(PrimeSieve sieve, int cand)
298: {
299: // Note, that this PrimeCollection is not made public
300: // via PrimeSieve. It is used only to test primality.
301: // It is a private constructor only for PrimeSieve.
302:
303: this.sieve = sieve;
304: primeRange = enumRange = null;
305:
306: start = IndexOf(cand - 1, 0, sieve.numberOfPrimes);
307: end = sieve.primes[start] == cand ? start + 1 : start;
308:
309: isPrime = start < end;
310: }
311:
312: /// <summary>
313: /// Affirms the primality of a collection of type (cand-1, cand].
314: /// </summary>
315: /// <returns>True if cand is prime.</returns>
316: public bool IsPrime()
317: {
318: return isPrime;
319: }
320:
321: /// <summary>
322: /// Provides an enumerator of the current prime number collection.
323: /// This enumerator is thread save.
324: /// </summary>
325: /// <returns>An enumerator of the current prime number collection.
326: /// </returns>
327: IEnumerator<int> IEnumerable<int>.GetEnumerator()
328: {
329: PrimeCollection result = this;
330:
331: // Make the Enumerator threadsave!
332: if (0 != Interlocked.CompareExchange(ref state, 1, 0))
333: {
334: result = new PrimeCollection(sieve, enumRange);
335: result.state = 1;
336: }
337: result.next = result.start;
338: return result;
339: }
340:
341: // Implementing the IEnumerator<int> interface requires implementing
342: // the nongeneric IEnumerator interface also. The Current property
343: // appears on both interfaces, and has different return types.
344:
345: /// <summary>
346: /// Provides an enumerator of the current prime number collection.
347: /// This enumerator is thread save.
348: /// </summary>
349: /// <returns>An enumerator of the current prime number
350: /// collection.</returns>
351: IEnumerator IEnumerable.GetEnumerator()
352: {
353: IEnumerable<int> enumerable = this;
354: return enumerable.GetEnumerator();
355: }
356:
357: /// <summary>
358: /// The next prime number in the collection.
359: /// </summary>
360: /// <returns> The current prime number in the collection.</returns>
361: object IEnumerator.Current
362: {
363: get { return sieve.primes[current]; }
364: }
365:
366: /// <summary>
367: /// The next prime number in the collection.
368: /// </summary>
369: /// <returns> The current prime number in the collection.</returns>
370: int IEnumerator<int>.Current
371: {
372: get { return sieve.primes[current]; }
373: }
374:
375: /// <summary>
376: /// Checks the current status of the enumerator.
377: /// </summary>
378: /// <returns>True iff there are more prime numbers to
379: /// be enumerated.</returns>
380: bool IEnumerator.MoveNext()
381: {
382: switch (state)
383: {
384: case 1:
385: if (next > end) goto case 2;
386: current = next++; return true;
387: case 2:
388: state = 2; return false;
389: default:
390: throw new InvalidOperationException();
391: }
392: }
393:
394: /// <summary>
395: /// Stop the enumerator before releasing resources.
396: /// </summary>
397: public void Dispose()
398: {
399: state = 2;
400: // This object will be cleaned up by the Dispose method.
401: // Therefore, a call to GC.SupressFinalize takes this object off
402: // the finalization queue and prevents finalization code for
403: // this object from executing a second time.
404: GC.SuppressFinalize(this);
405: }
406:
407: /// <summary>
408: /// Throws an NotImplementedException.
409: /// </summary>
410: void IEnumerator.Reset()
411: {
412: throw new NotImplementedException();
413: }
414:
415: /// <summary>
416: /// Identifies the index of a prime number.
417: /// Uses a (modified!) binary search.
418: /// </summary>
419: /// <param name="value">The prime number given.</param>
420: /// <param name="low">Lower bound for the index.</param>
421: /// <param name="high">Upper bound for the index.</param>
422: /// <returns>The index of the prime number.</returns>
423: private int IndexOf(int value, int low, int high)
424: {
425: int[] data = sieve.primes;
426:
427: while (low < high)
428: {
429: int mid = low + (high - low) / 2;
430:
431: if (data[mid] < value)
432: {
433: low = mid + 1;
434: }
435: else
436: {
437: high = mid;
438: }
439: }
440:
441: if (low >= data.Length)
442: {
443: return low;
444: }
445:
446: if (data[low] == value)
447: {
448: low++;
449: }
450:
451: return low;
452: }
453:
454: /// <summary>
455: /// Gives the prime numbers in the collection as an array.
456: /// </summary>
457: /// <returns>An array of prime numbers representing the collection.
458: /// </returns>
459: public int[] ToArray()
460: {
461: int primeCard = primeRange.Size();
462: int[] primeList = new int[primeCard];
463:
464: System.Array.Copy(sieve.primes, start, primeList, 0, primeCard);
465:
466: return primeList;
467: }
468:
469: /// <summary>
470: /// Gives the start and size of the prime range.
471: /// </summary>
472: /// <param name="begin">start of the prime range.</param>
473: /// <param name="size">Size of the prime range.</param>
474: /// <returns>Prime range is empty.</returns>
475: public bool GetSliceParameters(out int begin, out int size)
476: {
477: bool empty = 0 == primeRange.Max; // If the primeRange is empty...
478: begin = empty ? 0 : start;
479: size = empty ? 0 : primeRange.Size();
480: return empty;
481: }
482:
483: /// <summary>
484: /// Gets the number of primes in the current collection.
485: /// </summary>
486: /// <value>The number of primes in the current collection.
487: /// </value>
488: public int NumberOfPrimes
489: {
490: get
491: {
492: if (end == -1) return 0; // If primeRange is empty...
493: return primeRange.Size();
494: }
495: }
496:
497: /// <summary>
498: /// Gives the sieve range of the current collection.
499: /// </summary>
500: /// <value>Interval that was sieved.</value>
501: public PositiveRange SieveRange
502: {
503: get { return (PositiveRange)enumRange.Clone(); }
504: }
505:
506: /// <summary>
507: /// Gets the range of the indices of the prime numbers
508: /// in the collection.
509: /// </summary>
510: /// <value>Range of indices.</value>
511: public PositiveRange PrimeRange
512: {
513: get { return (PositiveRange) primeRange.Clone(); }
514: }
515:
516: /// <summary>
517: /// Writes the primes collection to a file.
518: /// </summary>
519: /// <param name="fileName">Name of the file to write to.
520: /// </param>
521: /// <returns>The full name of the file.</returns>
522: public string ToFile(string fileName)
523: {
524: FileInfo file = new FileInfo(@fileName);
525: StreamWriter primeReport = file.AppendText();
526: WritePrimes(primeReport);
527: primeReport.Close();
528: return file.FullName;
529: }
530:
531: /// <summary>
532: /// Writes the primes collection to a stream.
533: /// </summary>
534: /// <param name="file">The name of the file where the collection
535: /// is to be stored.</param>
536: private void WritePrimes(StreamWriter file)
537: {
538: file.WriteLine("SieveRange {0} : {1}",
539: enumRange.ToString(), enumRange.Size());
540: file.WriteLine("PrimeRange {0} : {1}",
541: primeRange.ToString(), primeRange.Size());
542: file.WriteLine("PrimeDensity {0:F3}",
543: (double)primeRange.Size() / (double)enumRange.Size());
544:
545: int primeOrdinal = start;
546: int primeLow = sieve.primes[start];
547: int lim = (primeLow / 100) * 100;
548:
549: foreach (int prime in this)
550: {
551: primeOrdinal++;
552: if (prime >= lim)
553: {
554: lim += 100;
555: file.Write(file.NewLine);
556: file.Write("<");
557: file.Write(primeOrdinal);
558: file.Write(".> ");
559: }
560: file.Write(prime);
561: file.Write(" ");
562: }
563: file.Write(file.NewLine);
564: file.Write(file.NewLine);
565: }
566:
567: } // endOfPrimeCollection
568: } } // endOfPrimeSieve
~bye