CSharp - 算法 C#: Atkin的实现

  显示原文与译文双语对照的内容
97 3

我想知道这里有没有人能够实现Atkin的筛选。

我正在尝试实现它,但不能完全围绕着它。这是我目前所拥有的。

public class Atkin : IEnumerable<ulong>
{
 private readonly List<ulong> primes;
 private readonly ulong limit;
 public Atkin(ulong limit)
 {
 this.limit = limit;
 primes = new List<ulong>();
 }
 private void FindPrimes()
 {
 var isPrime = new bool[limit + 1];
 var sqrt = Math.Sqrt(limit);
 for (ulong x = 1; x <= sqrt; x++)
 for (ulong y = 1; y <= sqrt; y++)
 {
 var n = 4*x*x + y*y;
 if (n <= limit && (n % 12 == 1 || n % 12 == 5))
 isPrime[n] ^= true;
 n = 3*x*x + y*y;
 if (n <= limit && n % 12 == 7)
 isPrime[n] ^= true;
 n = 3*x*x - y*y;
 if (x> y && n <= limit && n % 12 == 11)
 isPrime[n] ^= true;
 }
 for (ulong n = 5; n <= sqrt; n++)
 if (isPrime[n])
 for (ulong k = n*n; k <= limit; k *= k)
 isPrime[k] = false;
 primes.Add(2);
 primes.Add(3);
 for (ulong n = 5; n <= limit; n++)
 if (isPrime[n])
 primes.Add(n);
 }
 public IEnumerator<ulong> GetEnumerator()
 {
 if (!primes.Any())
 FindPrimes();
 foreach (var p in primes)
 yield return p;
 }
 IEnumerator IEnumerable.GetEnumerator()
 {
 return GetEnumerator();
 }
}

我只是试图"翻译"列出在维基百科上列出的伪代码,但是它不能正常工作。所以要么我误解了一些东西要么做了错误的事情。或者很可能都是。

有一个 List的前 500个素,我用作测试,我的实现失败,在 40 ( 或者 41) 。

索引 [40 ] 处的值不同
预期:179
但是:175

你能找到我的错误,你是否有一个可以以共享的实现,或者者两者都是?

我正在使用的精确测试如下所示:

public abstract class AtkinTests
{
 [Test]
 public void GetEnumerator_FirstFiveHundredNumbers_AreCorrect()
 {
 var sequence = new Atkin(2000000);
 var actual = sequence.Take(500).ToArray();
 var expected = First500;
 CollectionAssert.AreEqual(expected, actual);
 }
 private static readonly ulong[] First500 = new ulong[]
 {
 2, 3, 5, 7, 11, 13, 17,.. .
 };
}
时间:原作者:0个回答

117 0

这个代码:

for (ulong k = n*n; k <= limit; k *= k)
 isPrime[k] = false;

似乎不是对这种伪代码的忠实翻译:

is_prime(k) ← false, k ∈ {n², 2n², 3n²,.. ., limit}

你的代码似乎会运行于n * n,n ^ 4,n ^ 8,等等 换句话说,每次不要添加n 平方。尝试这个:

ulong nSquared = n * n;
for (ulong k = nSquared; k <= limit; k += nSquared)
 isPrime[k] = false;
原作者:
80 4

这里的是另一个实现。它使用BitArray来节省内存。Parallel.For 需要. NET 框架 4.

static List<int> FindPrimesBySieveOfAtkins(int max)
{
//var isPrime = new BitArray((int)max+1, false); 
//Can't use BitArray because of threading issues.
 var isPrime = new bool[max + 1];
 var sqrt = (int)Math.Sqrt(max);
 Parallel.For(1, sqrt, x =>
 {
 var xx = x * x;
 for (int y = 1; y <= sqrt; y++)
 {
 var yy = y * y;
 var n = 4 * xx + yy;
 if (n <= max && (n % 12 == 1 || n % 12 == 5))
 isPrime[n] ^= true;
 n = 3 * xx + yy;
 if (n <= max && n % 12 == 7)
 isPrime[n] ^= true;
 n = 3 * xx - yy;
 if (x> y && n <= max && n % 12 == 11)
 isPrime[n] ^= true;
 }
 });
 var primes = new List<int>() { 2, 3 };
 for (int n = 5; n <= sqrt; n++)
 {
 if (isPrime[n])
 {
 primes.Add(n);
 int nn = n * n;
 for (int k = nn; k <= max; k += nn)
 isPrime[k] = false;
 }
 }
 for (int n = sqrt + 1; n <= max; n++)
 if (isPrime[n])
 primes.Add(n);
 return primes;
}
原作者:
122 1

这里,它使用一个名为CompartmentalisedParallel的类,它允许你执行并行操作,但是控制线程数。但是,每次更改或者者为每个线程创建一个单独的for时,都需要锁定。

using System;
using System.Collections;
using System.Collections.Generic;
using System.Threading.Tasks;
namespace PrimeGenerator
{
 public class Atkin : Primes
 {
 protected BitArray mbaPrimes;
 protected bool mbThreaded = true;
 public Atkin(int limit)
 : this(limit, true)
 {
 }
 public Atkin(int limit, bool threaded)
 : base(limit)
 {
 mbThreaded = threaded;
 if (mbaPrimes == null) FindPrimes();
 }
 public bool Threaded
 {
 get
 {
 return mbThreaded;
 }
 }
 public override IEnumerator<int> GetEnumerator()
 {
 yield return 2;
 yield return 3;
 for (int lsN = 5; lsN <= msLimit; lsN += 2)
 if (mbaPrimes[lsN]) yield return lsN;
 }
 private void FindPrimes()
 {
 mbaPrimes = new BitArray(msLimit + 1, false);
 int lsSQRT = (int)Math.Sqrt(msLimit);
 int[] lsSquares = new int[lsSQRT + 1];
 for (int lsN = 0; lsN <= lsSQRT; lsN++)
 lsSquares[lsN] = lsN * lsN;
 if (Threaded)
 {
 CompartmentalisedParallel.For<BitArray>(
 1, lsSQRT + 1, new ParallelOptions(),
 (start, finish) => { return new BitArray(msLimit + 1, false); },
 (lsX, lsState, lbaLocal) =>
 {
 int lsX2 = lsSquares[lsX];
 for (int lsY = 1; lsY <= lsSQRT; lsY++)
 {
 int lsY2 = lsSquares[lsY];
 int lsN = 4 * lsX2 + lsY2;
 if (lsN <= msLimit && (lsN % 12 == 1 || lsN % 12 == 5))
 lbaLocal[lsN] ^= true;
 lsN -= lsX2;
 if (lsN <= msLimit && lsN % 12 == 7)
 lbaLocal[lsN] ^= true;
 if (lsX> lsY)
 {
 lsN -= lsY2 * 2;
 if (lsN <= msLimit && lsN % 12 == 11)
 lbaLocal[lsN] ^= true;
 }
 }
 return lbaLocal;
 },
 (lbaResult, start, finish) =>
 {
 lock (mbaPrimes) 
 mbaPrimes.Xor(lbaResult);
 },
 -1
 );
 }
 else
 {
 for (int lsX = 1; lsX <= lsSQRT; lsX++)
 {
 int lsX2 = lsSquares[lsX];
 for (int lsY = 1; lsY <= lsSQRT; lsY++)
 {
 int lsY2 = lsSquares[lsY];
 int lsN = 4 * lsX2 + lsY2;
 if (lsN <= msLimit && (lsN % 12 == 1 || lsN % 12 == 5))
 mbaPrimes[lsN] ^= true;
 lsN -= lsX2;
 if (lsN <= msLimit && lsN % 12 == 7)
 mbaPrimes[lsN] ^= true;
 if (lsX> lsY)
 {
 lsN -= lsY2 * 2;
 if (lsN <= msLimit && lsN % 12 == 11)
 mbaPrimes[lsN] ^= true;
 }
 }
 }
 }
 for (int lsN = 5; lsN <lsSQRT; lsN += 2)
 if (mbaPrimes[lsN])
 {
 var lsS = lsSquares[lsN];
 for (int lsK = lsS; lsK <= msLimit; lsK += lsS)
 mbaPrimes[lsK] = false;
 }
 }
 }
}

和CompartmentalisedParallel类:

using System;
using System.Threading.Tasks;
namespace PrimeGenerator
{
 public static class CompartmentalisedParallel
 {
 #region Int
 private static int[] CalculateCompartments(int startInclusive, int endExclusive, ref int threads)
 {
 if (threads == 0) threads = 1;
 if (threads == -1) threads = Environment.ProcessorCount;
 if (threads> endExclusive - startInclusive) threads = endExclusive - startInclusive;
 int[] liThreadIndexes = new int[threads + 1];
 liThreadIndexes[threads] = endExclusive - 1;
 int liIndexesPerThread = (endExclusive - startInclusive)/threads;
 for (int liCount = 0; liCount <threads; liCount++)
 liThreadIndexes[liCount] = liCount * liIndexesPerThread;
 return liThreadIndexes;
 }
 public static void For<TLocal>(
 int startInclusive, int endExclusive,
 ParallelOptions parallelOptions,
 Func<int, int, TLocal> localInit,
 Func<int, ParallelLoopState, TLocal, TLocal> body,
 Action<TLocal, int, int> localFinally,
 int threads)
 {
 int[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
 if (threads> 1)
 Parallel.For(
 0, threads, parallelOptions,
 (liThread, lsState) =>
 {
 TLocal llLocal = localInit(liThreadIndexes[liThread], liThreadIndexes[liThread + 1]);
 for (int liCounter = liThreadIndexes[liThread]; liCounter <liThreadIndexes[liThread + 1]; liCounter++)
 body(liCounter, lsState, llLocal);
 localFinally(llLocal, liThreadIndexes[liThread], liThreadIndexes[liThread + 1]);
 }
 );
 else
 {
 TLocal llLocal = localInit(startInclusive, endExclusive);
 for (int liCounter = startInclusive; liCounter <endExclusive; liCounter++)
 body(liCounter, null, llLocal);
 localFinally(llLocal, startInclusive, endExclusive);
 }
 }
 public static void For(
 int startInclusive, int endExclusive,
 ParallelOptions parallelOptions,
 Action<int, ParallelLoopState> body,
 int threads)
 {
 int[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
 if (threads> 1)
 Parallel.For(
 0, threads, parallelOptions,
 (liThread, lsState) =>
 {
 for (int liCounter = liThreadIndexes[liThread]; liCounter <liThreadIndexes[liThread + 1]; liCounter++)
 body(liCounter, lsState);
 }
 );
 else
 for (int liCounter = startInclusive; liCounter <endExclusive; liCounter++)
 body(liCounter, null);
 }
 public static void For(
 int startInclusive, int endExclusive,
 ParallelOptions parallelOptions,
 Action<int> body,
 int threads)
 {
 int[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
 if (threads> 1)
 Parallel.For(
 0, threads, parallelOptions,
 (liThread) =>
 {
 for (int liCounter = liThreadIndexes[liThread]; liCounter <liThreadIndexes[liThread + 1]; liCounter++)
 body(liCounter);
 }
 );
 else
 for (int liCounter = startInclusive; liCounter <endExclusive; liCounter++)
 body(liCounter);
 }
 public static void For(
 int startInclusive, int endExclusive,
 Action<int, ParallelLoopState> body,
 int threads)
 {
 For(startInclusive, endExclusive, new ParallelOptions(), body, threads);
 }
 public static void For(
 int startInclusive, int endExclusive,
 Action<int> body,
 int threads)
 {
 For(startInclusive, endExclusive, new ParallelOptions(), body, threads);
 }
 public static void For<TLocal>(
 int startInclusive, int endExclusive,
 Func<int, int, TLocal> localInit,
 Func<int, ParallelLoopState, TLocal, TLocal> body,
 Action<TLocal, int, int> localFinally,
 int threads)
 {
 For<TLocal>(startInclusive, endExclusive, new ParallelOptions(), localInit, body, localFinally, threads);
 }
 #endregion
 #region Long
 private static long[] CalculateCompartments(long startInclusive, long endExclusive, ref long threads)
 {
 if (threads == 0) threads = 1;
 if (threads == -1) threads = Environment.ProcessorCount;
 if (threads> endExclusive - startInclusive) threads = endExclusive - startInclusive;
 long[] liThreadIndexes = new long[threads + 1];
 liThreadIndexes[threads] = endExclusive - 1;
 long liIndexesPerThread = (endExclusive - startInclusive)/threads;
 for (long liCount = 0; liCount <threads; liCount++)
 liThreadIndexes[liCount] = liCount * liIndexesPerThread;
 return liThreadIndexes;
 }
 public static void For<TLocal>(
 long startInclusive, long endExclusive,
 ParallelOptions parallelOptions,
 Func<long, long, TLocal> localInit,
 Func<long, ParallelLoopState, TLocal, TLocal> body,
 Action<TLocal, long, long> localFinally,
 long threads)
 {
 long[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
 if (threads> 1)
 Parallel.For(
 0, threads, parallelOptions,
 (liThread, lsState) =>
 {
 TLocal llLocal = localInit(liThreadIndexes[liThread], liThreadIndexes[liThread + 1]);
 for (long liCounter = liThreadIndexes[liThread]; liCounter <liThreadIndexes[liThread + 1]; liCounter++)
 body(liCounter, lsState, llLocal);
 localFinally(llLocal, liThreadIndexes[liThread], liThreadIndexes[liThread + 1]);
 }
 );
 else
 {
 TLocal llLocal = localInit(startInclusive, endExclusive);
 for (long liCounter = startInclusive; liCounter <endExclusive; liCounter++)
 body(liCounter, null, llLocal);
 localFinally(llLocal, startInclusive, endExclusive);
 }
 }
 public static void For(
 long startInclusive, long endExclusive,
 ParallelOptions parallelOptions,
 Action<long, ParallelLoopState> body,
 long threads)
 {
 long[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
 if (threads> 1)
 Parallel.For(
 0, threads, parallelOptions,
 (liThread, lsState) =>
 {
 for (long liCounter = liThreadIndexes[liThread]; liCounter <liThreadIndexes[liThread + 1]; liCounter++)
 body(liCounter, lsState);
 }
 );
 else
 for (long liCounter = startInclusive; liCounter <endExclusive; liCounter++)
 body(liCounter, null);
 }
 public static void For(
 long startInclusive, long endExclusive,
 ParallelOptions parallelOptions,
 Action<long> body,
 long threads)
 {
 long[] liThreadIndexes = CalculateCompartments(startInclusive, endExclusive, ref threads);
 if (threads> 1)
 Parallel.For(
 0, threads, parallelOptions,
 (liThread) =>
 {
 for (long liCounter = liThreadIndexes[liThread]; liCounter <liThreadIndexes[liThread + 1]; liCounter++)
 body(liCounter);
 }
 );
 else
 for (long liCounter = startInclusive; liCounter <endExclusive; liCounter++)
 body(liCounter);
 }
 public static void For(
 long startInclusive, long endExclusive,
 Action<long, ParallelLoopState> body,
 long threads)
 {
 For(startInclusive, endExclusive, new ParallelOptions(), body, threads);
 }
 public static void For(
 long startInclusive, long endExclusive,
 Action<long> body,
 long threads)
 {
 For(startInclusive, endExclusive, new ParallelOptions(), body, threads);
 }
 public static void For<TLocal>(
 long startInclusive, long endExclusive,
 Func<long, long, TLocal> localInit,
 Func<long, ParallelLoopState, TLocal, TLocal> body,
 Action<TLocal, long, long> localFinally,
 long threads)
 {
 For<TLocal>(startInclusive, endExclusive, new ParallelOptions(), localInit, body, localFinally, threads);
 }
 #endregion
 }
}

素数基类:

using System.Collections;
using System.Collections.Generic;
namespace PrimeGenerator
{
 public abstract class Primes : IEnumerable<int>
 {
 protected readonly int msLimit;
 public Primes(int limit)
 {
 msLimit = limit;
 }
 public int Limit
 {
 get
 {
 return msLimit;
 }
 }
 public int Count
 {
 get
 {
 int liCount = 0;
 foreach (int liPrime in this)
 liCount++;
 return liCount;
 }
 }
 public abstract IEnumerator<int> GetEnumerator();
 IEnumerator IEnumerable.GetEnumerator()
 {
 return GetEnumerator();
 }
 }
}

通过执行以下操作来使用它:

 var lpPrimes = new Atkin(count, true);
 Console.WriteLine(lpPrimes.Count);
 Console.WriteLine(s.ElapsedMilliseconds);

但是,我发现Eratosthenes在所有情况下都更快,即使是在多线程模式下运行的四个核心 CPU,Atkin也会如此:

using System;
using System.Collections;
using System.Collections.Generic;
namespace PrimeGenerator
{
 public class Eratosthenes : Primes
 {
 protected BitArray mbaOddEliminated;
 public Eratosthenes(int limit)
 : base(limit)
 {
 if (mbaOddEliminated == null) FindPrimes();
 }
 public override IEnumerator<int> GetEnumerator()
 {
 yield return 2;
 for (int lsN = 3; lsN <= msLimit; lsN+=2)
 if (!mbaOddEliminated[lsN>>1]) yield return lsN;
 }
 private void FindPrimes()
 {
 mbaOddEliminated = new BitArray((msLimit>>1) + 1);
 int lsSQRT = (int)Math.Sqrt(msLimit);
 for (int lsN = 3; lsN <lsSQRT + 1; lsN += 2)
 if (!mbaOddEliminated[lsN>>1])
 for (int lsM = lsN*lsN; lsM <= msLimit; lsM += lsN<<1)
 mbaOddEliminated[lsM>>1] = true;
 }
 }
}

如果你让Atkin运行得更快,请让我知道 !

原作者:
...