查找数字的最大素因数的算法

Question

计算数字的最大素因数的最佳方法是什么？

我认为最有效的方法如下：

1. 找到能完全除尽的最低素数
2. 检查除法结果是否是素数
3. 如果不是，则找到下一个最低素数
4. 转到2。

我基于这个假设，因为计算小数更容易主要原因。 这大约是对的吗？ 我还应该研究哪些其他方法？

编辑：我现在意识到，如果有超过 2 个素数因子在起作用，我的方法是徒劳的，因为当结果是其他两个素数的乘积时，步骤 2 会失败，因此需要递归算法。

再次编辑：现在我意识到这仍然有效，因为最后找到的素数必须是最大的素数，因此对步骤 2 中的非素数结果的任何进一步测试都会导致更小的素数。

Answer 1

我知道这不是一个快速的解决方案。 希望发布为更容易理解的缓慢解决方案。

 public static long largestPrimeFactor(long n) {

        // largest composite factor must be smaller than sqrt
        long sqrt = (long)Math.ceil(Math.sqrt((double)n));

        long largest = -1;

        for(long i = 2; i <= sqrt; i++) {
            if(n % i == 0) {
                long test = largestPrimeFactor(n/i);
                if(test > largest) {
                    largest = test;
                }
            }
        }

        if(largest != -1) {
            return largest;
        }

        // number is prime
        return n;
    } 

Answer 2

下面的 C++ 算法不是最好的算法，但它适用于十亿以下的数字，而且速度相当快

#include <iostream>
using namespace std;

// ------ is_prime ------
// Determines if the integer accepted is prime or not
bool is_prime(int n){
    int i,count=0;
    if(n==1 || n==2)
      return true;
    if(n%2==0)
      return false;
    for(i=1;i<=n;i++){
    if(n%i==0)
        count++;
    }
    if(count==2)
      return true;
    else
      return false;
 }
 // ------ nextPrime -------
 // Finds and returns the next prime number
 int nextPrime(int prime){
     bool a = false;
     while (a == false){
         prime++;
         if (is_prime(prime))
            a = true;
     }
  return prime;
 }
 // ----- M A I N ------
 int main(){

      int value = 13195;
      int prime = 2;
      bool done = false;

      while (done == false){
          if (value%prime == 0){
             value = value/prime;
             if (is_prime(value)){
                 done = true;
             }
          } else {
             prime = nextPrime(prime);
          }
      }
        cout << "Largest prime factor: " << value << endl;
 }

Answer 3

这是我快速计算最大素因数的方法。
它基于修改后的x不包含非素数因子的事实。 为了实现这一点，我们一旦找到一个因子就除x。 那么剩下的就是返回最大的因数了。 它已经是最佳状态了。

代码（哈斯克尔）：

f max' x i | i > x = max'
           | x `rem` i == 0 = f i (x `div` i) i  -- Divide x by its factor
           | otherwise = f max' x (i + 1)  -- Check for the next possible factor

g x = f 2 x 2

Answer 4

在 C++ 中使用递归计算数字的最大素因数。 代码的工作原理解释如下：

int getLargestPrime(int number) {
    int factor = number; // assumes that the largest prime factor is the number itself
    for (int i = 2; (i*i) <= number; i++) { // iterates to the square root of the number till it finds the first(smallest) factor
        if (number % i == 0) { // checks if the current number(i) is a factor
            factor = max(i, number / i); // stores the larger number among the factors
            break; // breaks the loop on when a factor is found
        }
    }
    if (factor == number) // base case of recursion
        return number;
    return getLargestPrime(factor); // recursively calls itself
}

Answer 5

#python implementation
import math
n = 600851475143
i = 2
factors=set([])
while i<math.sqrt(n):
   while n%i==0:
       n=n/i
       factors.add(i)
   i+=1
factors.add(n)
largest=max(factors)
print factors
print largest

Answer 6

我正在使用继续将数字除以当前质因数的算法。

我在 python 3 中的解决方案：

def PrimeFactor(n):
    m = n
    while n%2==0:
        n = n//2
    if n == 1:         # check if only 2 is largest Prime Factor 
        return 2
    i = 3
    sqrt = int(m**(0.5))  # loop till square root of number
    last = 0              # to store last prime Factor i.e. Largest Prime Factor
    while i <= sqrt :
        while n%i == 0:
            n = n//i       # reduce the number by dividing it by it's Prime Factor
            last = i
        i+=2
    if n> last:            # the remaining number(n) is also Factor of number 
        return n
    else:
        return last
print(PrimeFactor(int(input()))) 

输入：10
输出：5

输入：600851475143
输出：6857

Answer 7

这是我在 c# 中的尝试。 最后打印出来的是该数字的最大素因数。 我检查了一下，它有效。

namespace Problem_Prime
{
  class Program
  {
    static void Main(string[] args)
    {
      /*
       The prime factors of 13195 are 5, 7, 13 and 29.

      What is the largest prime factor of the number 600851475143 ?
       */
      long x = 600851475143;
      long y = 2;
      while (y < x)
      {
        if (x % y == 0)
        {
          // y is a factor of x, but is it prime
          if (IsPrime(y))
          {
            Console.WriteLine(y);
          }
          x /= y;
        }

        y++;

      }
      Console.WriteLine(y);
      Console.ReadLine();
    }
    static bool IsPrime(long number)
    {
      //check for evenness
      if (number % 2 == 0)
      {
        if (number == 2)
        {
          return true;
        }
        return false;
      }
      //don't need to check past the square root
      long max = (long)Math.Sqrt(number);
      for (int i = 3; i <= max; i += 2)
      {
        if ((number % i) == 0)
        {
          return false;
        }
      }
      return true;
    }

  }
}

Answer 8

这是作为生成器提供的相同函数@Triptych，它也被稍微简化了。

def primes(n):
    d = 2
    while (n > 1):
        while (n%d==0):
            yield d
            n /= d
        d += 1

然后可以使用以下方法找到最大素数：

n= 373764623
max(primes(n))

以及使用以下方法找到的因子列表：

list(primes(n))

Answer 9

我认为最好将所有可能的小于 n 的素数存储在某个地方，然后迭代它们以找到最大的除数。 您可以从 prime-numbers.org 获取素数。

当然我假设你的数字不会太大:)

Answer 10

#include<stdio.h>
#include<conio.h>
#include<math.h>
#include <time.h>

factor(long int n)
{
long int i,j;
while(n>=4)
 {
if(n%2==0) {  n=n/2;   i=2;   }

 else
 { i=3;
j=0;
  while(j==0)
  {
   if(n%i==0)
   {j=1;
   n=n/i;
   }
   i=i+2;
  }
 i-=2;
 }
 }
return i;
 }

 void main()
 { 
  clock_t start = clock();
  long int n,sp;
  clrscr();
  printf("enter value of n");
  scanf("%ld",&n);
  sp=factor(n);
  printf("largest prime factor is %ld",sp);

  printf("Time elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);
  getch();
 }

Answer 11

Python 迭代方法，从数字中删除所有质因数

def primef(n):
    if n <= 3:
        return n
    if n % 2 == 0:
        return primef(n/2)
    elif n % 3 ==0:
        return primef(n/3)
    else:
        for i in range(5, int((n)**0.5) + 1, 6):
            #print i
            if n % i == 0:
                return primef(n/i)
            if n % (i + 2) == 0:
                return primef(n/(i+2))
    return n

Answer 12

n = abs(number);
result = 1;
if (n mod 2 == 0) {
 result = 2;
 while (n mod 2 = 0) n /= 2;
}
for(i=3; i<sqrt(n); i+=2) {
 if (n mod i == 0) {
   result = i;
   while (n mod i = 0)  n /= i;
 }
}
return max(n,result)

有一些模测试是多余的，因为如果所有因子 2 和 3 都被删除，n 永远不能被 6 除。 您只能允许 i 为素数，这在其他几个答案中有所显示。

实际上，您可以在这里交织埃拉托色尼筛：

• 首先创建整数列表
  为sqrt(n)。
• 在 for 循环中标记所有倍数
  i 到新的 sqrt(n) 为止
  prime，并使用 while 循环代替。
• 将 i 设置为下一个质数
  列表。

另请参阅此问题。

Answer 13

在“James Wang”的网络上找到了这个解决方案

public static int getLargestPrime( int number) {

    if (number <= 1) return -1;

    for (int i = number - 1; i > 1; i--) {
        if (number % i == 0) {
            number = i;
        }
    }
    return number;
}

Answer 14

使用筛子的质因数：

#include <bits/stdc++.h>
using namespace std;
#define N 10001  
typedef long long ll;
bool visit[N];
vector<int> prime;

void sieve()
{
            memset( visit , 0 , sizeof(visit));
            for( int i=2;i<N;i++ )
            {
                if( visit[i] == 0)
                {
                    prime.push_back(i);
                    for( int j=i*2; j<N; j=j+i )
                    {
                        visit[j] = 1;
                    }
                }
            }   
}
void sol(long long n, vector<int>&prime)
{
            ll ans = n;
            for(int i=0; i<prime.size() || prime[i]>n; i++)
            {
                while(n%prime[i]==0)
                {
                    n=n/prime[i];
                    ans = prime[i];
                }
            }
            ans = max(ans, n);
            cout<<ans<<endl;
}
int main() 
{
           ll tc, n;
           sieve();

           cin>>n;
           sol(n, prime);

           return 0;
}

Answer 15

我猜，除了执行因式分解之外，没有直接的方法，就像上面的例子所做的那样，即

在迭代中，您确定数字N的“小”因子f，然后继续简化问题“找到具有候选因子 >=f 的 N':=N/f 的最大素因子”。

从一定大小的 f 开始，如果您对减少的 N' 进行素性测试，则预期搜索时间会更短，以防万一确认您的 N'< /em> 已经是初始 N 的最大质因数。

Answer 16

这是我在 Clojure 中的尝试。 只走素数？的赔率和素数因子的素数，即。 筛子。 使用惰性序列有助于在需要之前生成值。

(defn prime? 
  ([n]
    (let [oddNums (iterate #(+ % 2) 3)]
    (prime? n (cons 2 oddNums))))
  ([n [i & is]]
    (let [q (quot n i)
          r (mod n i)]
    (cond (< n 2)       false
          (zero? r)     false
          (> (* i i) n) true
          :else         (recur n is)))))

(def primes 
  (let [oddNums (iterate #(+ % 2) 3)]
  (lazy-seq (cons 2 (filter prime? oddNums)))))

;; Sieve of Eratosthenes
(defn sieve
  ([n] 
    (sieve primes n))
  ([[i & is :as ps] n]
    (let [q (quot n i)
          r (mod n i)]
    (cond (< n 2)       nil
          (zero? r)     (lazy-seq (cons i (sieve ps q)))
          (> (* i i) n) (when (> n 1) (lazy-seq [n]))
          :else         (recur is n)))))

(defn max-prime-factor [n]
  (last (sieve n)))

Answer 17

JavaScript 代码：

'option strict';

function largestPrimeFactor(val, divisor = 2) { 
    let square = (val) => Math.pow(val, 2);

    while ((val % divisor) != 0 && square(divisor) <= val) {
        divisor++;
    }

    return square(divisor) <= val
        ? largestPrimeFactor(val / divisor, divisor)
        : val;
}

使用示例：

let result = largestPrimeFactor(600851475143);

以下是代码示例：

Answer 18

类似于 Triptych 的答案但也不同。 在此示例中，未使用列表或字典。 代码是用 Ruby 编写的

更新

感谢 Cooper Wolfe 的建议，算法得到了增强和优化。

def largest_prime_factor(number)
  i = 2
  while number > 1
    if number % i == 0
      number /= i
    elsif i > Math.sqrt(number)
      i = number
    else
      i += 1
    end
  end
  return i
end

largest_prime_factor(600851475143)
# => 6857

Answer 19

实际上，有几种更有效的方法可以找到大数的因子（对于较小的数，试除法效果相当好）。

如果输入数字有两个非常接近其平方根的因子，则一种非常快的方法称为 Fermat因式分解。 它利用恒等式 N = (a + b)(a - b) = a^2 - b^2 并且易于理解和实现。 不幸的是，总体来说速度不是很快。

最著名的对 100 位以内的数字进行因式分解的方法是二次筛。 额外的好处是，该算法的一部分可以通过并行处理轻松完成。

我听说过的另一种算法是 Pollard 的 Rho 算法。 一般来说，它不如二次筛有效，但似乎更容易实现。

---

一旦您决定如何将数字拆分为两个因子，这是我能想到的查找数字最大质因子的最快算法：

创建一个最初存储数字本身的优先级队列。 每次迭代，您都会从队列中删除最大的数字，并尝试将其拆分为两个因素（当然，不允许 1 成为这些因素之一）。 如果此步骤失败，则该数字是素数，您就有了答案！ 否则，您将这两个因素添加到队列中并重复。

Answer 20

我的答案是基于三联画的，但改进了很多。 它基于以下事实：除了 2 和 3 之外，所有素数的形式都是 6n-1 或 6n+1。

var largestPrimeFactor;
if(n mod 2 == 0)
{
    largestPrimeFactor = 2;
    n = n / 2 while(n mod 2 == 0);
}
if(n mod 3 == 0)
{
    largestPrimeFactor = 3;
    n = n / 3 while(n mod 3 == 0);
}

multOfSix = 6;
while(multOfSix - 1 <= n)
{
    if(n mod (multOfSix - 1) == 0)
    {
        largestPrimeFactor = multOfSix - 1;
        n = n / largestPrimeFactor while(n mod largestPrimeFactor == 0);
    }

    if(n mod (multOfSix + 1) == 0)
    {
        largestPrimeFactor = multOfSix + 1;
        n = n / largestPrimeFactor while(n mod largestPrimeFactor == 0);
    }
    multOfSix += 6;
}

我最近写了一篇博客文章解释该算法的工作原理。

我敢说，一种不需要素性测试（并且不需要构造筛子）的方法会比使用这些方法的方法运行得更快。 如果是这样的话，这可能是这里最快的算法。

Answer 21

C 算法中的递归

可以是

1. 检查 n 是否是一个因子或 t
2. 检查 n 是否是素数。 如果是这样，记住n
3. 递增n
4. 重复，直到n> sqrt(t)

下面是 C 语言问题的（尾）递归解决方案的示例：

#include <stdio.h>
#include <stdbool.h>

bool is_factor(long int t, long int n){
    return ( t%n == 0);
}

bool is_prime(long int n0, long int n1, bool acc){
    if ( n1 * n1 > n0 || acc < 1 )
        return acc;
    else
        return is_prime(n0, n1+2, acc && (n0%n1 != 0));
}

int gpf(long int t, long int n, long int acc){
    if (n * n > t)
        return acc;
    if (is_factor(t, n)){
        if (is_prime(n, 3, true))
            return gpf(t, n+2, n);
        else
            return gpf(t, n+2, acc);
    }
    else
        return gpf(t, n+2, acc);
}

int main(int argc, char ** argv){
    printf("%d\n", gpf(600851475143, 3, 0));
    return 0;
}

该解决方案由三个函数组成。 一个用于测试候选因子是否是一个因子，另一个用于测试该因子是否是素数，最后一个用于将这两个因子组合在一起。

这里的一些关键思想是：

1- 在 sqrt(600851475143) 处停止递归

2- 仅测试奇数的因子性

3- 仅测试奇数的候选因子的素数

Answer 22

    //this method skips unnecessary trial divisions and makes 
    //trial division more feasible for finding large primes

    public static void main(String[] args) 
    {
        long n= 1000000000039L; //this is a large prime number 
        long i = 2L;
        int test = 0;

        while (n > 1)
        {
            while (n % i == 0)
            {
                n /= i;     
            }

            i++;

            if(i*i > n && n > 1) 
            {
                System.out.println(n); //prints n if it's prime
                test = 1;
                break;
            }
        }

        if (test == 0)  
            System.out.println(i-1); //prints n if it's the largest prime factor
    }

Answer 23

最简单的解决方案是一对相互递归函数。

第一个函数生成所有素数：

1. 从所有大于 1 的自然数的列表开始。
2. 删除所有非素数。 也就是说，没有质因数的数字（除了它们本身）。 见下文。

第二个函数按升序返回给定数字 n 的质因数。

1. 列出所有素数（见上文）。
2. 删除所有不是 n 因数的数字。

n 的最大质因数是第二个函数给出的最后一个数字。

该算法需要一个惰性列表或具有按需调用语义的语言（或数据结构）。

为了澄清起见，这里是 Haskell 中上述内容的一个（低效）实现：

import Control.Monad

-- All the primes
primes = 2 : filter (ap (<=) (head . primeFactors)) [3,5..]

-- Gives the prime factors of its argument
primeFactors = factor primes
  where factor [] n = []
        factor xs@(p:ps) n =
          if p*p > n then [n]
          else let (d,r) = divMod n p in
            if r == 0 then p : factor xs d
            else factor ps n

-- Gives the largest prime factor of its argument
largestFactor = last . primeFactors

使其更快只是更聪明地检测哪些数字是素数和/或 n 的因子，但是算法保持不变。

Answer 24

所有数字都可以表示为素数的乘积，例如：

102 = 2 x 3 x 17
712 = 2 x 2 x 2 x 89

您可以通过简单地从 2 开始并继续除以直到结果不是您的数字的倍数来找到这些数字：

712 / 2 = 356 .. 356 / 2 = 178 .. 178 / 2 = 89 .. 89 / 89 = 1

使用此方法您不必实际计算任何素数：它们都将是素数，基于这样的事实：您已经尽可能地将数字与所有前面的数字分解。

number = 712;
currNum = number;    // the value we'll actually be working with
for (currFactor in 2 .. number) {
    while (currNum % currFactor == 0) {
        // keep on dividing by this number until we can divide no more!
        currNum = currNum / currFactor     // reduce the currNum
    }
    if (currNum == 1) return currFactor;    // once it hits 1, we're done.
}

Answer 25

在我看来，给出的算法的第二步并不是那么有效的方法。 你没有合理的期望它是素数。

此外，之前提出埃拉托斯特尼筛法的答案是完全错误的。 我刚刚编写了两个程序来分解 123456789。一个是基于 Sieve，一个是基于以下内容：

1)  Test = 2 
2)  Current = Number to test 
3)  If Current Mod Test = 0 then  
3a)     Current = Current Div Test 
3b)     Largest = Test
3c)     Goto 3. 
4)  Inc(Test) 
5)  If Current < Test goto 4
6)  Return Largest

这个版本比 Sieve 快 90 倍。

问题是，在现代处理器上，操作类型远不如操作数量重要，更不用说上面的算法可以在缓存中运行，而 Sieve 则不能。 筛子使用大量运算来剔除所有合数。

另请注意，我在确定因素时对其进行划分减少了必须测试的空间。

Answer 26

第 1 部分。 Pollard-Rho + 审判部门。

受到你的问题的启发，我决定在 Python 中实现我自己的因式分解（并找到最大素因数）版本。

据我所知，可能最简单但非常高效的因式分解算法是 Pollard 的 Rho算法。 它的运行时间最多为O(N^(1/4))，比O(N^(1/2))的时间快得多用于试除算法。 两种算法仅在复合（非质数）数的情况下才具有这些运行时间，这就是为什么应使用素性测试来过滤掉质数（不可因式分解）数。

我在代码中使用了以下算法：费马素性测试 ...，Pollard 的 Rho 算法 ..., 试除算法。 在运行 Pollard's Rho 之前使用 Fermat 素性测试来过滤掉素数。 Trial Division 被用作后备方案，因为 Pollard 的 Rho 在极少数情况下可能无法找到因子，特别是对于一些较小的数字。

显然，将一个数字完全分解为素因数排序列表后，最大的素因数将是该列表中的最后一个元素。 在一般情况下（对于任何随机数），除了完全因式分解数字之外，我不知道有任何其他方法可以找出最大素因数。

作为我的代码中的一个示例，我正在分解 Pi 的前 190 小数位，代码在 1 秒内分解该数字，并显示最大的质因数，即 165 位（545位）的大小！

在线试用！

def is_fermat_probable_prime(n, *, trials = 32):
    # https://en.wikipedia.org/wiki/Fermat_primality_test
    import random
    if n <= 16:
        return n in (2, 3, 5, 7, 11, 13)
    for i in range(trials):
        if pow(random.randint(2, n - 2), n - 1, n) != 1:
            return False
    return True

def pollard_rho_factor(N, *, trials = 16):
    # https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
    import random, math
    for j in range(trials):
        i, stage, y, x = 0, 2, 1, random.randint(1, N - 2)
        while True:
            r = math.gcd(N, x - y)
            if r != 1:
                break
            if i == stage:
                y = x
                stage <<= 1
            x = (x * x + 1) % N
            i += 1
        if r != N:
            return [r, N // r]
    return [N] # Pollard-Rho failed

def trial_division_factor(n, *, limit = None):
    # https://en.wikipedia.org/wiki/Trial_division
    fs = []
    while n & 1 == 0:
        fs.append(2)
        n >>= 1
    d = 3
    while d * d <= n and limit is None or d <= limit:
        q, r = divmod(n, d)
        if r == 0:
            fs.append(d)
            n = q
        else:
            d += 2
    if n > 1:
        fs.append(n)
    return fs

def factor(n):
    if n <= 1:
        return []
    if is_fermat_probable_prime(n):
        return [n]
    fs = trial_division_factor(n, limit = 1 << 12)
    if len(fs) >= 2:
        return sorted(fs[:-1] + factor(fs[-1]))
    fs = pollard_rho_factor(n)
    if len(fs) >= 2:
        return sorted([e1 for e0 in fs for e1 in factor(e0)])
    return trial_division_factor(n)

def demo():
    import time, math
    # http://www.math.com/tables/constants/pi.htm
    # pi = 3.
    #     1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679
    #     8214808651 3282306647 0938446095 5058223172 5359408128 4811174502 8410270193 8521105559 6446229489 5493038196
    #     4428810975 6659334461 2847564823 3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273
    #     7245870066 0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652 1384146951 9415116094
    #     3305727036 5759591953 0921861173 8193261179 3105118548 0744623799 6274956735 1885752724 8912279381 8301194912
    #     9833673362 4406566430 8602139494 6395224737 1907021798 6094370277 0539217176 2931767523 8467481846 7669405132
    #
    # n = first 190 fractional digits of Pi
    n =   1415926535_8979323846_2643383279_5028841971_6939937510_5820974944_5923078164_0628620899_8628034825_3421170679_8214808651_3282306647_0938446095_5058223172_5359408128_4811174502_8410270193_8521105559_6446229489
    print('Number:', n)
    tb = time.time()
    fs = factor(n)
    print('All Prime Factors:', fs)
    print('Largest Prime Factor:', f'({math.log2(fs[-1]):.02f} bits, {len(str(fs[-1]))} digits)', fs[-1])
    print('Time Elapsed:', round(time.time() - tb, 3), 'sec')

if __name__ == '__main__':
    demo()

输出：

Number: 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489
All Prime Factors: [3, 71, 1063541, 153422959, 332958319, 122356390229851897378935483485536580757336676443481705501726535578690975860555141829117483263572548187951860901335596150415443615382488933330968669408906073630300473]
Largest Prime Factor: (545.09 bits, 165 digits) 122356390229851897378935483485536580757336676443481705501726535578690975860555141829117483263572548187951860901335596150415443615382488933330968669408906073630300473
Time Elapsed: 0.593 sec

---

第 2 部分。 椭圆曲线法。

一段时间后，我决定通过从头开始实施更先进的椭圆曲线分解方法（称为 ECM）来改进我的帖子，请参阅维基百科文章 Lenstra 椭圆曲线分解。

此方法比我的回答帖子第 1 部分中描述的 Pollard Rho 快得多。 椭圆 ECM 方法的时间复杂度为 O(exp[(Sqrt(2) + o(1)) Sqrt(ln p ln ln p)])，其中 p 表示一个数的最小质因数。 而 Pollard Rho 的时间复杂度为 O(Sqrt(p))。 因此，对于足够大的最小 P 因子，ECM 要快得多。

ECM方法的步骤：

1. 检查数字是否小于2^16，然后通过试除法将其分解 方法。 返回结果。

2. 检查数字是否可能是素数，因为我使用 费马测试 与32 次试验。 要克服卡迈克尔数，您可以使用米勒拉宾测试。 如果数字是素数，则将其作为唯一因子返回。

3. 随机生成曲线参数A, X, Y并从曲线方程Y^2 = X^3 + AX + B (mod N)导出B 。 检查曲线是否正常，值 4 * A ** 3 - 27 * B ** 2 应为非零。

4. 通过埃拉托斯特尼筛法生成小素数，这些素数低于我们的界限。 每个素数都会提升到某个小幂，这个提升的素数将被称为 K。当它小于某个 Bound2 时，我会进行幂提升，在我的例子中，它是 Sqrt(Bound)。

5. 计算椭圆点乘法P = k * P，其中 K 取自上一步，P 为 (X, Y)。 根据Wiki进行计算。

6. 点乘法使用模逆，计算GCD(SomeValue, N) 根据Wiki。 如果这个 GCD 不是 1，那么它给出 N 的非 1 因子，因此在这种情况下，我通过异常并从 ECM 分解算法返回该因子。

7. 如果 Bound 之前的所有素数都相乘并且没有给出任何因子，则使用另一条随机曲线和更大的 Bound 再次重新运行 ECM 分解算法（上面的 1.-6.）。 在我的代码中，我通过将 256 添加到旧边界来获取新边界。

我的 ECM 代码中使用了以下子算法，所有提到的子算法都有相应维基百科文章的链接： Trial除法分解、费马概率检验、埃拉托斯特尼筛法（素数生成器），欧几里得算法（计算最大公约数，GCD），扩展欧几里得算法（带有 Bezu 系数的 GCD），模乘逆，从右到左二进制指数（用于椭圆点乘法），椭圆曲线算术（点加法和乘法），Lenstra 椭圆曲线分解。

与我的答案的第 1 部分不同，我对 Pi 数字 的数字进行因式分解，在第 2 部分我创建一个由 n = Prime(24) * Prime(35) * Prime(37) * ... 等组成的特殊数字，这意味着数字是 24 的随机素数的乘积位、35 位、37 位等...这个自定义数字在视觉上更令人印象深刻，可以显示算法功能。

正如我的答案的第 1 部分一样，我还使用多种方法来比较速度，并用更简单的方法分解出更小的因素。 所以在下面的代码中我使用 Trial Division + Pollard Rho + 椭圆曲线法。

在下面的代码之后，请查看控制台输出，以找出我的代码输出的好东西。

在线尝试！

def is_fermat_probable_prime(n, *, trials = 32):
    # https://en.wikipedia.org/wiki/Fermat_primality_test
    import random
    if n <= 16:
        return n in (2, 3, 5, 7, 11, 13)
    for i in range(trials):
        if pow(random.randint(2, n - 2), n - 1, n) != 1:
            return False
    return True

def pollard_rho_factor(N, *, trials = 8, limit = None):
    # https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
    import random, math
    if N <= 1:
        return []
    if is_fermat_probable_prime(N):
        return [N]
    for j in range(trials):
        i, stage, y, x = 0, 2, 1, random.randint(1, N - 2)
        while True:
            r = math.gcd(N, x - y)
            if r != 1:
                break
            if i == stage:
                y = x
                stage <<= 1
            x = (x * x + 1) % N
            i += 1
            if limit is not None and i >= limit:
                return [N] # Pollard-Rho failed
        if r != N:
            return sorted(pollard_rho_factor(r, trials = trials, limit = limit) +
                pollard_rho_factor(N // r, trials = trials, limit = limit))
    return [N] # Pollard-Rho failed

def trial_division_factor(n, *, limit = None):
    # https://en.wikipedia.org/wiki/Trial_division
    fs = []
    while n & 1 == 0:
        fs.append(2)
        n >>= 1
    d = 3
    while d * d <= n and limit is None or d <= limit:
        q, r = divmod(n, d)
        if r == 0:
            fs.append(d)
            n = q
        else:
            d += 2
    if n > 1:
        fs.append(n)
    return fs
    
def ecm_factor(N0, *, verbose = False):
    # https://en.wikipedia.org/wiki/Lenstra_elliptic-curve_factorization
    import math, random, time; gmpy2 = None
    #import gmpy2
    def GenPrimes_SieveOfEratosthenes(end):
        # https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
        composite = [False] * (end // 2) # Allocate for odd numbers only
        for p in range(3, int(end ** 0.5 + 3), 2):
            if composite[p >> 1]:
                continue
            for i in range(p * p, end, p * 2):
                composite[i >> 1] = True
        yield 2
        for p in range(3, end, 2):
            if not composite[p >> 1]:
                yield p
    def GCD(a, b):
        # https://en.wikipedia.org/wiki/Euclidean_algorithm
        # return math.gcd(a, b)
        while b != 0:
            a, b = b, a % b
        return a
    def EGCD(a, b):
        # https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
        if gmpy2 is None:
            ro, r, so, s = a, b, 1, 0
            while r != 0:
                ro, (q, r) = r, divmod(ro, r)
                so, s = s, so - q * s
            return ro, so, (ro - so * a) // b
        else:
            return tuple(map(int, gmpy2.gcdext(a, b)))
    def ModularInverse(a, n):
        # https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
        # return pow(a, -1, n)
        g, s, t = EGCD(a, n)
        if g != 1:
            raise ValueError(a)
        return s % n
    def EllipticCurveAdd(N, A, B, X0, Y0, X1, Y1):
        # https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication
        if X0 == X1 and Y0 == Y1:
            # Double
            l = ((3 * X0 ** 2 + A) * ModularInverse(2 * Y0, N)) % N
            x = (l ** 2 - 2 * X0) % N
            y = (l * (X0 - x) - Y0) % N
        else:
            # Add
            l = ((Y1 - Y0) * ModularInverse(X1 - X0, N)) % N
            x = (l ** 2 - X0 - X1) % N
            y = (l * (X0 - x) - Y0) % N
        return x, y
    def EllipticCurveMul(N, A, B, X, Y, k):
        # https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
        assert k >= 2, k
        k -= 1
        BX, BY = X, Y
        while k != 0:
            if k & 1:
                X, Y = EllipticCurveAdd(N, A, B, X, Y, BX, BY)
            BX, BY = EllipticCurveAdd(N, A, B, BX, BY, BX, BY)
            k >>= 1
        return X, Y
    
    bound_start = 1 << 9
    
    def Main(N, *, bound = bound_start, icurve = 0):
        def NextFactorECM(x):
            return Main(x, bound = bound + bound_start, icurve = icurve + 1)
        def PrimePow(p, *, bound2 = int(bound ** 0.5 + 1.01)):
            mp = p
            while True:
                mp *= p
                if mp >= bound2:
                    return mp // p
        
        if N <= 1:
            return []
            
        if N < (1 << 16):
            fs = trial_division_factor(N)
            if verbose and len(fs) >= 2:
                print('Factors from TrialDiv:', fs, flush = True)
            return fs
            
        if is_fermat_probable_prime(N):
            return [N]
        
        if verbose:
            print(f'Curve {icurve:>4},  bound 2^{math.log2(bound):>7.3f}', flush = True)
        
        while True:
            X, Y, A = [random.randrange(N) for i in range(3)]
            B = (Y ** 2 - X ** 3 - A * X) % N
            if 4 * A ** 3 - 27 * B ** 2 == 0:
                continue
            break
        
        for p in GenPrimes_SieveOfEratosthenes(bound):
            k = PrimePow(p)
            try:
                X, Y = EllipticCurveMul(N, A, B, X, Y, k)
            except ValueError as ex:
                g = GCD(ex.args[0], N)
                assert g > 1, g
                if g != N:
                    if verbose:
                        print(f'Factor from ECM: {g} ({math.log2(g):.1f} bits)', flush = True)
                    return sorted(NextFactorECM(g) + NextFactorECM(N // g))
                else:
                    return NextFactorECM(N)
        
        return NextFactorECM(N)
    
    return Main(N0)

def factor(n):
    if n <= 1:
        return []
    if is_fermat_probable_prime(n):
        return [n]
    fs1 = trial_division_factor(n, limit = 1 << 12)
    fs1, n2 = fs1[:-1], fs1[-1]
    print(len(fs1), 'factors from TrialDivision')
    fs2 = pollard_rho_factor(n2, limit = 1 << 17)
    fs2, n3 = fs2[:-1], fs2[-1]
    print(len(fs2), 'factors from PollardRho')
    fs3 = ecm_factor(n3, verbose = True)
    print(len(fs3), 'factors from ECM')
    return sorted(fs1 + fs2 + fs3)

def demo():
    import time, math, random
    def Prime(bits):
        while True:
            x = random.randrange(1 << (bits - 1), 1 << bits)
            if is_fermat_probable_prime(x):
                return x
    
    # http://www.math.com/tables/constants/pi.htm
    # pi = 3.
    #     1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679
    #     8214808651 3282306647 0938446095 5058223172 5359408128 4811174502 8410270193 8521105559 6446229489 5493038196
    #     4428810975 6659334461 2847564823 3786783165 2712019091 4564856692 3460348610 4543266482 1339360726 0249141273
    #     7245870066 0631558817 4881520920 9628292540 9171536436 7892590360 0113305305 4882046652 1384146951 9415116094
    #     3305727036 5759591953 0921861173 8193261179 3105118548 0744623799 6274956735 1885752724 8912279381 8301194912
    #     9833673362 4406566430 8602139494 6395224737 1907021798 6094370277 0539217176 2931767523 8467481846 7669405132
    #
    # n =   1415926535_8979323846_2643383279_5028841971_6939937510_5820974944_5923078164_0628620899_8628034825_3421170679_8214808651_3282306647_0938446095_5058223172_5359408128_4811174502_8410270193_8521105559_6446229489_5493038196_4428810975_6659334461_2847564823_3786783165_2712019091_4564856692_3460348610_4543266482_1339360726_0249141273_7245870066_0631558817_4881520920_9628292540_9171536436_7892590360_0113305305_4882046652_1384146951_9415116094_3305727036_5759591953_0921861173_8193261179_3105118548_0744623799_6274956735_1885752724_8912279381_8301194912_9833673362_4406566430_8602139494_6395224737_1907021798_6094370277_0539217176_2931767523_8467481846_7669405132
    #       141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132
    
    n = Prime(9) * Prime(10) * Prime(11) * Prime(20) * Prime(24) * Prime(35) * Prime(37) * Prime(38) * Prime(120)
    print('Number:', n)
    tb = time.time()
    fs = factor(n)
    print('All Prime Factors:', fs)
    print('Largest Prime Factor:', f'({math.log2(fs[-1]):.02f} bits, {len(str(fs[-1]))} digits)', fs[-1])
    if len(fs) >= 2:
        print('2nd-largest Prime Factor:', f'({math.log2(fs[-2]):.02f} bits, {len(str(fs[-2]))} digits)', fs[-2])
    print('Time Elapsed:', round(time.time() - tb, 3), 'sec')

if __name__ == '__main__':
    demo()

控制台输出：

Number: 3020823358956369790763854998578637168366763837218991014777892420353187988302225517459334041

3 factors from TrialDivision
2 factors from PollardRho

Curve    0,  bound 2^  9.000
Curve    1,  bound 2^ 10.000
Curve    2,  bound 2^ 10.585
Factor from ECM: 20028139561 (34.2 bits)
Curve    3,  bound 2^ 11.000
Curve    4,  bound 2^ 11.322
Curve    5,  bound 2^ 11.585
Curve    6,  bound 2^ 11.807
Curve    7,  bound 2^ 12.000
Curve    8,  bound 2^ 12.170
Factor from ECM: 96583780901 (36.5 bits)
Curve    9,  bound 2^ 12.322
Curve   10,  bound 2^ 12.459
Curve   11,  bound 2^ 12.585
Curve   12,  bound 2^ 12.700
Curve   13,  bound 2^ 12.807
Curve   14,  bound 2^ 12.907
Curve   15,  bound 2^ 13.000
Curve   16,  bound 2^ 13.087
Curve   17,  bound 2^ 13.170
Factor from ECM: 239171423261 (37.8 bits)
4 factors from ECM

All Prime Factors: [397, 1021, 1459, 754333, 16156687, 20028139561,
    96583780901, 239171423261, 905908369146483365552973334921981697]
Largest Prime Factor: (119.45 bits, 36 digits) 905908369146483365552973334921981697
2nd-largest Prime Factor: (37.80 bits, 12 digits) 239171423261

Time Elapsed: 17.156 sec

Answer 27

首先计算一个存储素数的列表，例如 2 3 5 7 11 13 ...

每次对一个数进行素因式分解时，使用 Triptych 的实现，但迭代这个素数列表而不是自然整数。

Answer 28

对于 Java：

对于 int 值：

public static int[] primeFactors(int value) {
    int[] a = new int[31];
    int i = 0, j;
    int num = value;
    while (num % 2 == 0) {
        a[i++] = 2;
        num /= 2;
    }
    j = 3;
    while (j <= Math.sqrt(num) + 1) {
        if (num % j == 0) {
            a[i++] = j;
            num /= j;
        } else {
            j += 2;
        }
    }
    if (num > 1) {
        a[i++] = num;
    }
    int[] b = Arrays.copyOf(a, i);
    return b;
}

对于 long 值：

static long[] getFactors(long value) {
    long[] a = new long[63];
    int i = 0;
    long num = value;
    while (num % 2 == 0) {
        a[i++] = 2;
        num /= 2;
    }
    long j = 3;
    while (j <= Math.sqrt(num) + 1) {
        if (num % j == 0) {
            a[i++] = j;
            num /= j;
        } else {
            j += 2;
        }
    }
    if (num > 1) {
        a[i++] = num;
    }
    long[] b = Arrays.copyOf(a, i);
    return b;
}

Answer 29

这可能并不总是更快，但更乐观的是您找到一个大素数：

1. N 是您的数字
2. 如果它是素数，则 return(N)
3. 计算素数直到 < code>Sqrt(N)
4. 按降序排列素数（最大的在前）
   • 如果N 能被素数整除，则Return(Prime)

编辑：在步骤 3 中，您可以使用埃拉托斯特尼筛法或阿特金斯筛法或任何您喜欢的方法，但筛子本身不会找到最大的素因数。 （这就是为什么我不会选择 SQLMenace 的帖子作为官方答案......）

Answer 30

这是我所知道的最好的算法（在 Python 中）

def prime_factors(n):
    """Returns all the prime factors of a positive integer"""
    factors = []
    d = 2
    while n > 1:
        while n % d == 0:
            factors.append(d)
            n /= d
        d = d + 1

    return factors


pfs = prime_factors(1000)
largest_prime_factor = max(pfs) # The largest element in the prime factor list

在最坏的情况下（当输入是素数时），上述方法的运行时间为 O(n) 。

编辑：
下面是 O(sqrt(n)) 版本，如评论中所建议。 这是代码。

def prime_factors(n):
    """Returns all the prime factors of a positive integer"""
    factors = []
    d = 2
    while n > 1:
        while n % d == 0:
            factors.append(d)
            n /= d
        d = d + 1
        if d*d > n:
            if n > 1: factors.append(n)
            break
    return factors


pfs = prime_factors(1000)
largest_prime_factor = max(pfs) # The largest element in the prime factor list