Program.cs

using System;
using System.Collections.Concurrent;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Runtime.CompilerServices;
using System.Threading.Tasks;
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs;
using BenchmarkDotNet.Running;

namespace bench
{
    public class Program
    {
        public static void Main()
        {
            BenchmarkRunner.Run<Primes>();
        }
        [SimpleJob(launchCount: 1, warmupCount: 2, targetCount: 5, runtimeMoniker: RuntimeMoniker.Net60)]
        public class Primes
        {
            private List<uint> _numbers;
            private static readonly ulong[] witnesses1 = new ulong[] { 9345883071009581737 };
            private const uint witnesses1Max = 341531;
            private static readonly ulong[] witnesses2 = new ulong[] { 15, 13393019396194701 };
            private const uint witnesses2Max = 716169301;
            private static readonly ulong[] witnesses3 = new ulong[] { 2, 7, 61 };

            [GlobalSetup]
            public void Setup()
            {
                const int limit = 50_000_000;
                _numbers = Enumerable.Range(1, limit).Select(i => (uint)i).ToList();
                Console.WriteLine(GetPrimeListWithParallel(_numbers).Count);
            }

            [Benchmark]
            public IList<uint> PrimeListClassic()
            {
                return GetPrimeListClassic(_numbers);
            }

            [Benchmark]
            public IList<uint> PrimeListWithParallelMillerRabin()
            {
                return GetPrimeListWithParallel(_numbers);
            }

            private static IList<uint> GetPrimeListClassic(IList<uint> numbers)
            {
                var primeList = new List<uint>();

                for (int i = 0; i < numbers.Count; i++)
                {
                    var currentItem = numbers[i];
                    if (IsPrime(currentItem)) primeList.Add(currentItem);
                }
                return primeList;
            }

            private static IList<uint> GetPrimeListWithParallel(IList<uint> numbers)
            {
                var primeNumbers = new ConcurrentBag<uint>();

                Parallel.ForEach(numbers, number =>
                {
                    if (IsPrime(number))
                    {
                        primeNumbers.Add(number);
                    }
                });

                return primeNumbers.ToList();
            }

            public static bool IsPrime(uint value)
            {
                if ((value & 1) == 0) return value == 2;
                if ((value % 3) == 0) return value == 3;
                if ((value % 5) == 0) return value == 5;
                if ((value % 7) == 0) return value == 7;
                if (value == 1) return false;
                if (value < 1000000)
                {
                    for (uint i = 11, skip = 2; (i <= ushort.MaxValue) && (i * i <= value); i += skip, skip = 6 - skip)
                    {
                        if (value % i == 0) return false;
                    }
                    return true;
                };
                if (((value % 11) == 0) || ((value % 13) == 0) || ((value % 17) == 0)
                    || ((value % 19) == 0) || ((value % 23) == 0) || ((value % 29) == 0)
                    || ((value % 31) == 0) || ((value % 37) == 0) || ((value % 41) == 0)
                    || ((value % 43) == 0) || ((value % 47) == 0) || ((value % 53) == 0)
                    || ((value % 59) == 0) || ((value % 61) == 0) || ((value % 67) == 0)
                    || ((value % 71) == 0) || ((value % 73) == 0) || ((value % 79) == 0)
                    || ((value % 83) == 0) || ((value % 89) == 0) || ((value % 97) == 0)
                    )
                {
                    return (value <= 97);
                }
                return IsPrimeMR(value);
            }

            public static bool IsPrimeMR(uint value)
            {
                if (value >= witnesses2Max) return InternalIsPrimeMR(value, witnesses3);
                else if (value >= witnesses1Max) return InternalIsPrimeMR(value, witnesses2);
                else return InternalIsPrimeMR(value, witnesses1);
            }

            private static bool InternalIsPrimeMR(uint value, ulong[] witnesses)
            {
                uint valLessOne = value - 1;
                uint d = valLessOne / 2;
                uint s = 1;
                while ((d & 1) == 0)
                {
                    d /= 2;
                    s++;
                }
                for (uint i = 0; i < witnesses.Length; i++)
                {
                    uint a;
                    ulong aL = witnesses[i];
                    if (aL >= value)
                    {
                        aL %= value;
                        if (aL < 2) continue;
                    }
                    a = (uint)aL;
                    if (a == valLessOne) continue;
                    uint x = InternalModPow(a, d, value);
                    if (x == 1) continue;
                    for (uint r = 1; (x != valLessOne) && (r < s); r++)
                    {
                        x = (uint)((x * (ulong)x) % value);
                        if (x == 1) return false;
                    }
                    if (x != valLessOne) return false;
                }
                return true;
            }
            public static uint ModPow(uint b, uint e, uint m)
            {
                if (b == 0 || m == 0) return 0;
                return InternalModPow(b, e, m);
            }
            private static readonly uint sqrtMax = (uint)Math.Sqrt(uint.MaxValue);
            private static uint InternalModPow(uint b, uint e, uint m)
            {
                switch (e)
                {
                    case 0: return 1;
                    case 1: return b % m;
                    case 2: return (b <= sqrtMax) ? (b * b) % m : (uint)((b * (ulong)b) % m);
                    default:
                        uint result = 1;
                        while (true)
                        {
                            if ((e & 1) != 0)
                            {
                                result = (uint)((result * (ulong)b) % m);
                                if (e == 1) return result;
                            }
                            e /= 2;
                            b = (uint)((b * (ulong)b) % m);
                        }
                }
            }
        }

    }
}