Быстрый алгоритм определения НЕ простоты числа C#

Почему быстрый:

  • в 10 раз быстрее алгоритма Миллера-Рабина 
  • в 500 раз быстрее вероятностного теста
  • в 380 раз быстрее проверки на делимость перебором

MayBePrimeSpeed                :0,002691
IsMillerRabinPrimeSpeed     :0,027961
IsProbablePrimeSpeed         :1,393796
IsPrimeSpeed                         :1,031095
Positive Answer Count:
MayBePrimePos                   :8549/35007
IsMillerRabinPrimePos        :6763/35007
IsProbablePrimePos            :6762/35007

IsPrimePos                            :6763/35007


 Возьмем число 390, у нас есть 77 простых чисел, меньших 390. 
Четыре из них - 2,3,5,13 - не более, чем мусор, дающий лишние срабатывания алгоритма. В самом деле, если x%390==3, то x делится на 3, а значит, оно составное. 
Остальные 73 числа годятся.
Получается такая картина. Для простого x величина x%390 - некоторое число, взаимно простое с 390. Причём распределены эти остатки более-менее равномерно (на картинке они выглядят как красные столбики). Всего этих остатков 96.
Если x - простое число, то с вероятностью 73/96 остаток будет простым числом, и алгоритм честно скажет "скорее, простое". 
С вероятностью 23/96 -  25% - остаток окажется составным, и алгоритм ошибётся. Поэтому результат выполнения нужно перепроверять другими алгоритмами, которые будут медленнее данного.
Если x - составное, то вероятность того, что число объявят "скорее, простым" будет равна 77/390-1/ln(x) (очень мало, чем больше число тем меньше) (первое слагаемое - вероятность того, что остаток оказался простым, а второе - доля простых чисел в окрестности x).
Можно легко избавиться от ошибки первого вида: для этого надо в массив charr положить числа взаимно простые с 390. Тогда если алгоритм скажет "составное", то так и есть.
Можно вместо 390 взять N=30030=2*3*5*7*11*13. В таком случае отсеиваться, как составные, будет 4 числа из 5, а ложных отсеиваний простых чисел не будет вообще.

 Алгоритм вероятностного определения составности числа:
class MayBePrime
    {
        static bool[] charr = null;
        static int prime = 30030;//30030 для быстрого старта, ниже точность. Произведение первых простых чисел
        static public bool isMayBePrime(System.Numerics.BigInteger x)
        {
            if (charr == null)
            {
                charr = new bool[prime];
                for (int i = 0; i < prime; i++)
                {
                    if (NOD(prime, i) == 1)
                    {
                        charr[i] = true;
                    }
                    else
                    {
                        charr[i] = false;
                    }
                }
            }
            BigInteger n = x % prime;
            int n2 = 0;
            n2 = (int)n;
            if (charr[n2])
            {
                return true;
            }
            else
            {
                return false;
            }
        }
  static public int NOD(int a, int b)
        {
            if (a == 0) return b;
            if (b == 0) return a;
            if (a == b) return a;
            if (a == 1 || b == 1) return 1;
            if ((a % 2 == 0) && (b % 2 == 0)) return 2 * NOD(a / 2, b / 2);
            if ((a % 2 == 0) && (b % 2 != 0)) return NOD(a / 2, b);
            if ((a % 2 != 0) && (b % 2 == 0)) return NOD(a, b / 2);
            return NOD(b, Math.Abs(a - b));
        }
    }


За помощь в доработке и анализа алгоритма благодарность Mrrl


Алгоритмы с которыми производилось сравнение:


static public bool isPrime(int x)
        {
            //System.Numerics.BigInteger sqrt = System.Numerics.BigInteger.Pow(x, 1 / 2);
            for (int i = 2; i < x; i++)
            {
                if ((x % i) == 0) return false;
            }
            return true;
        }
static class Amicable
    {
        public static BigInteger Pow2(BigInteger exp)
        {
            BigInteger res = 1;
 
            for (BigInteger i = 0; i < exp; i++)
                res *= 2;
 
            return res;
        }
 
        public static BigInteger ExpMod(BigInteger value, BigInteger exp, BigInteger mod)
        {
            BigInteger ullResult = 1;
            BigInteger ullValue = value;
 
            while (exp > 0)
            {
                if (exp % 2 != 0)
                { // odd
                    ullResult *= ullValue;
                    ullResult %= mod;
                }
 
                ullValue *= ullValue;
                ullValue %= mod;
                exp /= 2;
            }
 
            return (ullResult);
        }
 
        static public bool IsMillerRabinPrime(BigInteger prime)
        {
            // randomWitness = witness that the "prime" is NOT composite
            // 1 < randomWitness < prime-1
            long totalWitness = 15;
            BigInteger[] randomWitness = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 51001, 351011 };
            BigInteger primeMinusOne = prime - 1;
            BigInteger oddMultiplier;
            long bitLength;
            long i, j;
            BigInteger result;
            // find oddMultiplier, defined as "prime - 1 = (2^B) * M"
            // get bitLength (B) and find the oddMultiplier (M)
 
            // init value multiplier
            oddMultiplier = primeMinusOne;
 
            bitLength = 0; // reset
            while (oddMultiplier % 2 == 0)
            {
                oddMultiplier /= 2;
                bitLength++;
            }
            for (i = 0; i < totalWitness; i++)
            {
                if (randomWitness[i] == prime)
                    return true;
 
                // init value of result = (randomWitness ^ oddMultiplier) mod prime
                result = ExpMod(randomWitness[i], oddMultiplier, prime);
 
                // is y = 1 ?
                if (result == 1)
                    continue; // maybe prime
 
                // is y = prime-1 ?
                if (result == primeMinusOne)
                    continue; // maybe prime
 
                // loop to get AT LEAST one "result = primeMinusOne"
                for (j = 1; j <= bitLength; j++)
                {
                    // square of result
                    result = ExpMod(result, 2, prime);
 
                    // is result = primeMinusOne ?
                    if (result == primeMinusOne)
                        break; // maybe prime
                }
 
                if (result != primeMinusOne)
                    return false; // COMPOSITE
            }
 
            // it may be pseudoprime/prime
            return true;
        }
 
    }
    public static class BigIntegerExtensions
    {
        public static bool IsProbablePrime(this BigInteger source, int certainty)
        {
            if (source == 2 || source == 3)
                return true;
            if (source < 2 || source % 2 == 0)
                return false;
 
            BigInteger d = source - 1;
            int s = 0;
 
            while (d % 2 == 0)
            {
                d /= 2;
                s += 1;
            }
 
            // There is no built-in method for generating random BigInteger values.
            // Instead, random BigIntegers are constructed from randomly generated
            // byte arrays of the same length as the source.
            RandomNumberGenerator rng = RandomNumberGenerator.Create();
            byte[] bytes = new byte[source.ToByteArray().LongLength];
            BigInteger a;
 
            for (int i = 0; i < certainty; i++)
            {
                do
                {
                    // This may raise an exception in Mono 2.10.8 and earlier.
                    // http://bugzilla.xamarin.com/show_bug.cgi?id=2761
                    rng.GetBytes(bytes);
                    a = new BigInteger(bytes);
                }
                while (a < 2 || a >= source - 2);
 
                BigInteger x = BigInteger.ModPow(a, d, source);
                if (x == 1 || x == source - 1)
                    continue;
 
                for (int r = 1; r < s; r++)
                {
                    x = BigInteger.ModPow(x, 2, source);
                    if (x == 1)
                        return false;
                    if (x == source - 1)
                        break;
                }
 
                if (x != source - 1)
                    return false;
            }
 
            return true;
        }






Программа сравнения скорости:


using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Numerics;
using System.IO;
using System.Security.Cryptography;
 
namespace ConsoleApplication1
{
    class Program
    {
        static int kk = 0;
        static DateTime dt1;
        static DateTime dt2;
        static TimeSpan ts;
 
        static void StartTimer()
        {
            dt1 = DateTime.Now;
        }
        static void EndTimer()
        {
            dt2 = DateTime.Now;
            ts = dt2 - dt1;
        }
 
 
 
        static void Main(string[] args)
        {
            string[] primeliststr = File.ReadAllLines(@"C:\2 .txt");
            List<BigInteger> primelist = new List<BigInteger>();
            int kk = 0;
            foreach (string item in primeliststr)
            {
                if (item!="")
                {
                    primelist.Add(int.Parse(item));
                }
                kk++;
                if (kk>5000)
                {
                    break;
                }
            }
            int primescount = primelist.Count ;
 
            var MayBePrimeSpeed = new List<int>();
            var IsMillerRabinPrimeSpeed = new List<int>();
            var IsProbablePrimeSpeed = new List<int>();
            var IsPrimeSpeed = new List<int>();
            bool MayBePrimeResult;
            bool IsMillerRabinPrimeResult;
            bool IsProbablePrimeResult;
            bool IsPrimeResult;
 
            int MayBePrimePos = 0;
            int IsMillerRabinPrimePos=0;
            int IsProbablePrimePos = 0;
            int IsPrimePos = 0;
 
            int MayBePrimeNeg = 0;
            int IsMillerRabinPrimeNeg = 0;
            int IsProbablePrimeNeg = 0;
            int IsPrimeNeg = 0;
 
            var m = primelist.Count*2;
            for (int i = 0; i < m ; i++)
            {
                primelist.Add(13 * i);//добавляем заведомо составные числа
            }
 
            for (int i = 0; i < m; i++)
            {
                primelist.Add(13 ^ 10 * i);//добавляем заведомо составные большие числа  
            }
 
            Random rnd = new Random(DateTime.Now.Millisecond);
            for (int i = 0; i < m; i++)
            {
                primelist.Add(rnd.Next(1000) ^ 20 * i);//добавляем cлучайные большие числа 
            }
 
            for (int i = 10; i < primelist.Count ; i++)
   {
                StartTimer();
                MayBePrimeResult = MayBePrime.isMayBePrime(primelist[i]);
                EndTimer();
                MayBePrimeSpeed.Add(ts.Milliseconds);
                if (MayBePrimeResult)
                {
                    MayBePrimePos++;
                }
                else
                {
                    MayBePrimeNeg++;
                }
 
                StartTimer();
                IsMillerRabinPrimeResult = MayBePrime.IsMillerRabinPrime(primelist[i]);
                EndTimer();
                IsMillerRabinPrimeSpeed.Add(ts.Milliseconds);
                if (IsMillerRabinPrimeResult)
                {
                    IsMillerRabinPrimePos++;
                }
                else
                {
                    IsMillerRabinPrimeNeg++;
                }
 
                StartTimer();
                IsProbablePrimeResult = MayBePrime.IsProbablePrime(primelist[i], 10);
                EndTimer();
                IsProbablePrimeSpeed.Add(ts.Milliseconds);
                if (IsProbablePrimeResult)
                {
                    IsProbablePrimePos++;
                }
                else
                {
                    IsProbablePrimeNeg++;
                }
 
                StartTimer();
                IsPrimeResult = MayBePrime.isPrime(primelist[i]);
                EndTimer();
                IsPrimeSpeed.Add(ts.Milliseconds);
                if (IsPrimeResult)
                {
                    IsPrimePos++;
                }
                else
                {
                    IsPrimeNeg++;
                }
 
                if (i%150==0)
                {
                    Console.Clear();
                    Console.WriteLine("Current Position :{0}", primelist[i].ToString());
 
                    Console.WriteLine("MayBePrimeSpeed                :{0}", MayBePrimeSpeed.Average());
                    Console.WriteLine("IsMillerRabinPrimeSpeed        :{0}", IsMillerRabinPrimeSpeed.Average());
                    Console.WriteLine("IsProbablePrimeSpeed           :{0}", IsProbablePrimeSpeed.Average());
                    Console.WriteLine("IsPrimeSpeed                   :{0}", IsPrimeSpeed.Average());
 
                    Console.WriteLine("Positive Answer Count:");
 
                    Console.WriteLine("MayBePrimePos                :{0}/{1}", MayBePrimePos, primelist.Count);
                    Console.WriteLine("IsMillerRabinPrimePos        :{0}/{1}", IsMillerRabinPrimePos, primelist.Count);
                    Console.WriteLine("IsProbablePrimePos           :{0}/{1}", IsProbablePrimePos, primelist.Count);
                    Console.WriteLine("IsPrimePos                   :{0}/{1}", IsPrimePos, primelist.Count);
                }
   }
            Console.Clear();
            Console.WriteLine("MayBePrimeSpeed                :{0}", MayBePrimeSpeed.Average());
            Console.WriteLine("IsMillerRabinPrimeSpeed        :{0}", IsMillerRabinPrimeSpeed.Average());
            Console.WriteLine("IsProbablePrimeSpeed           :{0}", IsProbablePrimeSpeed.Average());
            Console.WriteLine("IsPrimeSpeed                   :{0}", IsPrimeSpeed.Average());
 
            Console.WriteLine("Positive Answer Count:");
 
            Console.WriteLine("MayBePrimePos                :{0}/{1}", MayBePrimePos, primelist.Count);
            Console.WriteLine("IsMillerRabinPrimePos        :{0}/{1}", IsMillerRabinPrimePos, primelist.Count);
            Console.WriteLine("IsProbablePrimePos           :{0}/{1}", IsProbablePrimePos, primelist.Count);
            Console.WriteLine("IsPrimePos                   :{0}/{1}", IsPrimePos, primelist.Count);
 
            Console.WriteLine("Primes count: {0}", IsPrimePos);
            Console.ReadLine();
        }
     }
}
 








































Комментарии








Отправить комментарий




















Популярные сообщения