`   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 = 2;`
` 119:          primes = 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 &lt;= 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 &lt;= 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`