【问题标题】:How can I test for primality?我如何测试素数?
【发布时间】:2010-10-12 06:21:09
【问题描述】:

我正在编写一个带有一些素数相关方法的小库。由于我已经完成了基础工作(又名工作方法),现在我正在寻找一些优化。 当然,互联网是一个很好的地方。然而,我偶然发现了一个舍入问题,我想知道如何解决这个问题。

在我用来测试一个数字的素数的循环中,搜索到 sqrt(n) 而不是 n/2 甚至 n - 1 会更有效。但是由于舍入问题,一些数字被跳过,因此一些素数是跳过!例如,第 10000 个素数应为:104729,但“优化”版本最终为:103811。

一些代码(我知道它可以进行更多优化,但我一次只能处理一件事):

/// <summary>
/// Method for testing the primality of a number e.g.: return IsPrime(29);
/// History:
/// 1. Initial version, most basic form of testing: m smaller then n -1
/// 2. Implemented m smaller then sqrt(n), optimization due to prime factoring
/// </summary>
/// <param name="test">Number to be tested on primality</param>
/// <returns>True if the number is prime, false otherwise</returns>
public static bool IsPrime(int test)
{
    // 0 and 1 are not prime numbers
    if (test == 0 || test == 1) return false;

    // 2 and 3 are prime numbers
    if (test == 2) return true;

    // all even numbers, save 2, are not prime
    if (test % 2 == 0) return false;

    double squared = Math.Sqrt(test);
    int flooredAndSquared = Convert.ToInt32(Math.Floor(squared));

    // start with 5, make increments of 2, even numbers do not need to be tested
    for (int idx = 3; idx < flooredAndSquared; idx++)
    {
        if (test % idx == 0)
        {
            return false;
        }
    }
    return true;
}

我知道平方部分失败了(或者我失败了),我也尝试了 Math.Ceiling,结果大致相同。

【问题讨论】:

  • 您的 for 循环似乎从 3 开始,并以 1 为增量;您的评论指出它从 5 开始并以 2 递增。
  • 这个问题似乎是题外话,因为它是关于数论的。试试 math.stackexchange.com。
  • squared 不是平方根结果的正确变量名。 “平方”是指提高到二次方;平方根提高到 1/2 次方。也许叫它sqrt_test 什么的。

标签: c# math primes


【解决方案1】:

重复模式操作会运行得很慢。使用eratosthenes网格按顺序获取素数列表。

/*
The Sieve Algorithm
http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
*/
numbers = new MyBitArray(limit, true);
for (long i = 2; i < limit; i++)
    if (numbers[i])
        for (long j = i * 2; j < limit; j += i)
            numbers[j] = false;
        }

public class MyBitArray: IDisposable
    {
        byte[] bytes;
        public MyBitArray(long limit, bool defaultValue = false)
        {
            long byteCount = (limit & 7) == 0 ? limit >> 3 : (limit >> 3) + 1;
            this.bytes = new byte[byteCount];
            for(long i = 0; i < byteCount; i++)
            {
                bytes[i] = (defaultValue == true ? (byte)0xFF : (byte)0x00);
            }
            this.limit = limit;
        }

        public MyBitArray(long limit, byte[] bytes)
        {
            this.limit = limit;
            this.bytes = bytes;
        }

        public bool this[long index]
        {
            get
            {
                return getValue(index);
            }
            set
            {
                setValue(index, value);
            }
        }
        
        bool getValue(long index)
        {
            if (index < 8)
            {
                return getBit(bytes[0], (byte)index);
            }

            long byteIndex = (index & 7) == 0 ? ((index >> 3) - 1) : index >> 3;
            byte bitIndex = (byte)(index & 7);
            return getBit(bytes[byteIndex], bitIndex);
        }
        void setValue(long index, bool value)
        {
            if (index < 8)
            {
                bytes[0] = setBit(bytes[0], (byte)index, value);
                return;
            }

            long byteIndex = (index & 7) == 0 ? (index >> 3) - 1 : index >> 3;
            byte bitIndex = (byte)(index & 7);
            
            bytes[byteIndex] = setBit(bytes[byteIndex], bitIndex, value);
        }

        bool getBit(byte byt, byte index)
        {
            return ((byt & (1 << index)) >> index) == 1;
        }

        byte setBit(byte byt, byte index, bool value)
        {
            return (byte)((byt & ~(1 << index)) + (value ? 1 << index : 0));
        }

        public void Dispose()
        {
            GC.Collect(2, GCCollectionMode.Optimized);
        }

        private long limit;
        public long Limit { get { return limit; } }
        public byte[] Bytes { get { return this.bytes; } } 
    }

但是,我建议您使用更好的素数测试方法。对于 64 位数字,无论数字多大,它都会以毫秒为单位给出准确的结果。

public static bool IsPrime(ulong number)
{
    return number == 2 
        ? true 
        : (BigInterger.ModPow(2, number, number) == 2 
            ? (number & 1 != 0 && BinarySearchInA001567(number) == false) 
            : false)
}

public static bool BinarySearchInA001567(ulong number)
{
    // Is number in list?
    // todo: Binary Search in A001567 (https://oeis.org/A001567) below 2 ^ 64
    // Only 2.35 Gigabytes as a text file http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html
}

【讨论】:

    【解决方案2】:

    首先,素数从 2 开始。2 和 3 是素数。素数不能被 2 或 3 整除。其余素数的形式为 6k-1 和 6k+1。请注意,您应该检查直到 SQRT(input) 的数字。这种方法非常有效。希望对你有帮助。

    public class Prime {
    
        public static void main(String[] args) {
            System.out.format("%d is prime: %s.\n", 199, isPrime(199)); // Prime
            System.out.format("%d is prime: %s.\n", 198, isPrime(198)); // Not prime
            System.out.format("%d is prime: %s.\n", 104729, isPrime(104729)); // Prime
            System.out.format("%d is prime: %s.\n", 104727, isPrime(982443529)); // Prime
        }
    
        /**
         * Tells if a number is prime or not.
         *
         * @param input the input
         * @return If the input is prime or not
         */
        private boolean isPrime(long input) {
        if (input <= 1) return false; // Primes start from 2
        if (input <= 3) return true; // 2 and 3 are primes
        if (input % 2 == 0 || input % 3 == 0) return false; // Not prime if dividable by 2 or 3
        // The rest of the primes are in the shape of 6k-1 and 6k+1
        for (long i = 5; i <= Math.sqrt(input); i += 6) if (input % i == 0 || input % (i + 2) == 0) return false;
        return true;
        }
    
    }
    

    【讨论】:

    • FWIW,我已将其转换为 Excel 的 VBA。如果有人感兴趣,我已经在下面发布了。
    【解决方案3】:

    如果其他人感兴趣,这是我将 Mohammad 的上述程序转换为 VBA 的过程。我添加了一个检查以排除 1、0 和负数,因为它们都被定义为非质数。

    我只在 Excel VBA 中测试过:

    Function IsPrime(input_num As Long) As Boolean
        Dim i As Long
        If input_num < 2 Then '1, 0, and negative numbers are all defined as not prime.
            IsPrime = False: Exit Function
        ElseIf input_num = 2 Then
            IsPrime = True: Exit Function '2 is a prime
        ElseIf input_num = 3 Then
            IsPrime = True: Exit Function '3 is a prime.
        ElseIf input_num Mod 2 = 0 Then
            IsPrime = False: Exit Function 'divisible by 2, so not a prime.
        ElseIf input_num Mod 3 = 0 Then
            IsPrime = False: Exit Function 'divisible by 3, so not a prime.
        Else
            'from here on, we only need to check for factors where
            '6k ± 1 = square root of input_num:
            i = 5
            Do While i * i <= input_num
                If input_num Mod i = 0 Then
                    IsPrime = False: Exit Function
                ElseIf input_num Mod (i + 2) = 0 Then
                    IsPrime = False: Exit Function
                End If
                i = i + 6
            Loop
            IsPrime = True
        End If
    End Function
    

    【讨论】:

      【解决方案4】:

      您可能想查看Fermat's little theorem

      这是Algorithms by S. Dasgupta, C.H. Papadimitriou, and U.V. Vazirani 一书中的伪代码,其中 n 是您要测试素数的数字。

      Pick a positive integer a < n at random
      if a^n-1 is equivalent to 1 (mod n)
         return yes
      else
         return no
      

      实施费马定理应该比筛解法更快。然而,有一些卡迈克尔数通过了费马的检验并且不是素数。有解决方法。我建议咨询Section 1.3 in the fore mentioned book。这完全是关于素数测试,可能对您有所帮助。

      【讨论】:

      • 你想这样做几次才能有真正的信心。
      • 是的,但它足够快,你可以做到这一点。我编辑了答案以提及 Carmichael 的数字。
      • Soloway-Strassen 和 Miller-Rabin 素性检验几乎在各个方面都优于费马小定理;两者都可以简单地扩展到确定性(不仅仅是概率)测试,尽管运行时不是最佳的。不要为 FLT 烦恼。
      • @kquinn,您在技术上是正确的,但我提出 FLT 是因为它是 Miller-Rabin 的基础。我链接的书继续解释 FLT 的弱点,然后将其扩展到米勒拉宾。在我发布这个答案之前,我也没有看到关于 Miller-Rabin 的标记。
      • 什么是n(小写)和什么是N(大写)?
      【解决方案5】:

      这对于测试素数非常有效(vb.net)

      Dim rnd As New Random()
      Const one = 1UL
      
          Function IsPrime(ByVal n As ULong) As Boolean
              If n Mod 3 = 0 OrElse n Mod 5 = 0 OrElse n Mod 7 = 0 OrElse n Mod 11 = 0 OrElse n Mod 13 = 0 OrElse n Mod 17 = 0 OrElse n Mod 19 = 0 OrElse n Mod 23 = 0 Then
                 return false
              End If
      
              Dim s = n - one
      
              While s Mod 2 = 0
                  s >>= one
              End While
      
              For i = 0 To 10 - 1
                  Dim a = CULng(rnd.NextDouble * n + 1)
                  Dim temp = s
                  Dim m = Numerics.BigInteger.ModPow(a, s, n)
      
                  While temp <> n - one AndAlso m <> one AndAlso m <> n - one
                      m = (m * m) Mod n
                      temp = temp * 2UL
                  End While
      
                  If m <> n - one AndAlso temp Mod 2 = 0 Then
                      Return False
                  End If
              Next i
      
              Return True
          End Function
      

      【讨论】:

        【解决方案6】:

        我认为Prime numbers and primality testing 很有用,而且 AKS 算法听起来很有趣,即使与基于概率的测试相比它并不是特别实用。

        【讨论】:

          【解决方案7】:

          遗憾的是,我之前没有尝试过算法方法。但是如果你想有效地实现你的方法,我建议做一些缓存。创建一个数组来存储小于定义阈值的所有素数,填充此数组,并在其中搜索/使用它。

          在下面的例子中,找出一个数是否为素数在最佳情况下是 O(1)(即当该数小于或等于maxPrime,对于 64K 缓冲区为 821,461 时),并且是针对其他情况进行了一些优化(通过检查前 820,000 个中仅 64K 个数字的 mod——大约 8%)。

          (注意:不要将此答案视为“最佳”方法。更多的是关于如何优化实施的示例。)

          public static class PrimeChecker
          {
              private const int BufferSize = 64 * 1024; // 64K * sizeof(int) == 256 KB
          
              private static int[] primes;
              public static int MaxPrime { get; private set; }
          
              public static bool IsPrime(int value)
              {
                  if (value <= MaxPrime)
                  {
                      return Array.BinarySearch(primes, value) >= 0;
                  }
                  else
                  {
                      return IsPrime(value, primes.Length) && IsLargerPrime(value);
                  }
              }
          
              static PrimeChecker()
              {
                  primes = new int[BufferSize];
                  primes[0] = 2;
                  for (int i = 1, x = 3; i < primes.Length; x += 2)
                  {
                      if (IsPrime(x, i))
                          primes[i++] = x;
                  }
                  MaxPrime = primes[primes.Length - 1];
              }
          
              private static bool IsPrime(int value, int primesLength)
              {
                  for (int i = 0; i < primesLength; ++i)
                  {
                      if (value % primes[i] == 0)
                          return false;
                  }
                  return true;
              }
          
              private static bool IsLargerPrime(int value)
              {
                  int max = (int)Math.Sqrt(value);
                  for (int i = MaxPrime + 2; i <= max; i += 2)
                  {
                      if (value % i == 0)
                          return false;
                  }
                  return true;
              }
          }
          

          【讨论】:

          • 这种技术被称为memoization,以防万一有人想搜索它。
          【解决方案8】:

          正如 Mark 所说,Miller-Rabin 测试实际上是一个非常好的方法。另一个参考(带有伪代码)是关于它的Wikipedia article

          应该注意,虽然它是概率性的,但通过仅测试极少数情况,您可以确定一个数字对于 int(和近乎 long)范围内的数字是否是素数。有关详细信息,请参阅 this part of that Wikipedia articlethe Primality Proving reference

          我还建议阅读有关模幂运算的 this article,否则在尝试进行 Miller-Rabin 测试时您将处理非常非常大的数字...

          【讨论】:

            【解决方案9】:

            这是我为解决欧拉问题而编写的一个不错的函数:

            private static long IsPrime(long input)
            {
                if ((input % 2) == 0)
                {
                    return 2;
                }
                else if ((input == 1))
                {
                    return 1;
                }
                else
                {
                    long threshold = (Convert.ToInt64(Math.Sqrt(input)));
                    long tryDivide = 3;
                    while (tryDivide < threshold)
                    {
                        if ((input % tryDivide) == 0)
                        {
                            Console.WriteLine("Found a factor: " + tryDivide);
                            return tryDivide;
                        }
                        tryDivide += 2;
                    }
                    Console.WriteLine("Found a factor: " + input);
                    return -1;
                }
            }
            

            【讨论】:

            • 与 OP 相同的错误 - 这应该是“tryDivide
            • 点了,我已经调整回原来的答案了。
            • 抱歉,如果 sqrt 碰巧返回比需要的“少”一点(例如,n = 9、sqrt(n) == 2.99999999、floor -> 2、算法认为它是主要的......我有点搞混了。抱歉
            【解决方案10】:
            private static bool IsPrime(int number) {
                if (number <= 3)
                    return true;
                if ((number & 1) == 0)
                    return false;
                int x = (int)Math.Sqrt(number) + 1;
                for (int i = 3; i < x; i += 2) {
                    if ((number % i) == 0)
                        return false;
                }
                return true;
            }
            

            我不能更快地得到它......

            【讨论】:

            • 不错的尝试。查看我的答案以获取有关如何使其更快的示例。
            【解决方案11】:

            您的 for 循环应如下所示:

            for (int idx = 3; idx * idx <= test; idx++) { ... }
            

            这样,您就可以避免浮点计算。应该跑得更快,它会更准确。这就是为什么 for 条件语句只是一个布尔表达式:它使这样的事情成为可能。

            【讨论】:

            • 将其标记下来,因为它既不快也不准确。一个整数的平方根的底正是他想要的。需要 IEEE 浮点算法才能在整数上产生正确的结果,即 sqrt(25) 不能是 4.999999。而且它更慢,因为您在循环中引入了以前没有的乘法。最后,您还引入了一个错误,因为idx * idx 可能会溢出并产生负值,从而导致无限循环。考虑 test = 2147483647 和 idx = 46340 和 46341。
            【解决方案12】:

            我不知道这是否正是您想要的,但如果您真的关心速度,那么您应该研究测试素数的概率方法,而不是使用筛子。 Rabin-Miller 是 Mathematica 使用的概率素性测试。

            【讨论】:

              【解决方案13】:

              我在这里发布了一个使用筛子或埃拉托色尼计算素数的课程:

              Is the size of an array constrained by the upper limit of int (2147483647)?

              【讨论】:

              • Eratosthenes 的筛子非常快,但只有知道要测试的初选的上限在哪里,才能使用它。
              • 完全没有,看代码。它逐步扩展范围,因此仅受用于存储素数的数据类型的容量限制。在这种情况下很长,所以限制是 9223372036854775807。
              • 一点也不:您可以创建一个筛子近似值列表,并“尽可能多地”获取。你当然需要惰性求值,就像函数式语言一样。我已经在 C# 中使用 yield 语句编写了一个实现,据我所知它运行良好。不要随身携带笔记本电脑,所以我必须稍后再回来发布实际答案。 (如果你们都想要的话。)
              【解决方案14】:

              试试这个...

              if (testVal == 2) return true;
              if (testVal % 2 == 0) return false;
              
              for (int i = 3; i <= Math.Ceiling(Math.Sqrt(testVal)); i += 2)
              {
                 if (testVal % i == 0)
                     return false;
              }
              
              return true;
              

              我已经用过很多次了.. 不如筛子快.. 但它有效。

              【讨论】:

              • 我相信i &lt; Math.Ceiling(或i &lt;= Math.Floor)就足够了。你不需要i &lt;= Math.Ceiling
              【解决方案15】:

              我猜这是你的问题:

              for (int idx = 3; idx < flooredAndSquared; idx++)
              

              这应该是

              for (int idx = 3; idx <= flooredAndSquared; idx++)
              

              所以你不会得到平方数作为素数。此外,您可以使用“idx += 2”而不是“idx++”,因为您只需要测试奇数(正如您在上面的评论中所写...)。

              【讨论】:

                【解决方案16】:

                试试sieve of eratosthenes -- 这应该可以解决获取根和浮点问题。

                至于floor,您可以通过ceiling获得更好的服务。

                【讨论】:

                • 天花板会在我认为的素数平方上给你一个误报:)
                • 如果你爬到天花板上就不行。
                猜你喜欢
                • 1970-01-01
                • 2015-07-22
                • 2014-08-19
                • 1970-01-01
                • 1970-01-01
                • 1970-01-01
                • 2010-11-19
                • 1970-01-01
                相关资源
                最近更新 更多