Prime Factorization

Introduction

Given a positive integer \(N \in \mathbf{N}_{+}\), the task is to quickly find one of its non-trivial divisors.

Consider the naive algorithm: divisors come in pairs, and all divisors of \(N\) can be divided into two intervals: \([2, \sqrt N]\) and \([\sqrt N+1,N)\). By iterating through all numbers in \([2, \sqrt N]\) and using division, we can find at least two divisors. The time complexity of this method is \(O(\sqrt N)\).

When \(N\ge10^{18}\), the runtime of this algorithm becomes unacceptable, and we desire a more efficient algorithm. One idea is to use a randomized approach: guess whether a number is a divisor of \(N\). With good luck, we can solve the problem in \(O(1)\) time complexity. However, for \(N\ge10^{18}\), the probability of a successful guess is \(\frac{1}{10^{18}}\), and the expected number of guesses is \(10^{18}\). If we restrict guessing to \([2,\sqrt N]\), the success probability increases. We hope to optimize the guessing process.

Naive Algorithm

The simplest algorithm is to iterate through \([2, \sqrt N]\).

C++Python

vector<int> breakdown(int N) {
  vector<int> result;
  for (int i = 2; i * i <= N; i++) {
    if (N % i == 0) {  // If i divides N, then i is a prime factor of N.
      while (N % i == 0) N /= i;
      result.push_back(i);
    }
  }
  if (N != 1) {  // After the operations, N is left as a prime number.
    result.push_back(N);
  }
  return result;
}

def breakdown(N):
    result = []
    for i in range(2, int(sqrt(N)) + 1):
        if N % i == 0:  # If i divides N, then i is a prime factor of N.
            while N % i == 0:
                N //= i
            result.append(i)
    if N != 1:  # After the operations, N is left as a prime number.
        result.append(N)
    return result

We can prove that all elements in result are exactly the prime factors of N.

Proof that result contains all prime factors of \(N\)

First, consider the changes in N. At the end of the loop for i, after executing while(N % i == 0) N /= i, i no longer divides N. Moreover, each time a factor is removed, it is guaranteed that N still divides \(N\). These two points ensure that at the start of the loop for i, N is a factor of \(N\) and is not divisible by any integer less than i.

Next, prove that all elements in result are factors of \(N\). When the loop reaches i, the condition for storing i in result is N % i == 0, which means i divides N. Since N is a factor of \(N\), it follows that i is a factor of \(N\). At the end of the loop for i, if N is not 1, it is also stored in result. According to the previous discussion, it must also be a factor of \(N\).

Then, prove that all elements in result are prime. Suppose there exists a composite number \(K\) in result. Then there must exist an i not exceeding \(\sqrt K\) such that i is a factor of K. Such a \(K\) cannot be stored as some i during the loop because the first paragraph shows that when the loop reaches \(K\), N is not divisible by any i less than \(K\). Such a \(K\) also cannot be added after the loop ends because the loop termination condition is i * i > N, meaning all i not exceeding \(\sqrt K\) have been traversed. Moreover, as stated above, these i cannot divide the current N, i.e., \(K\).

Finally, prove that all prime factors of \(N\) must appear in result. Suppose \(p\) is a prime factor of \(N\) but does not appear in result. According to the above discussion, \(p\) cannot be any i that appeared during the loop. Let i be the last i before the loop exits. Then i is strictly less than \(p\), and after the loop, N is not divisible by any previous i. Hence, \(p\) divides N. Therefore, if the final N is greater than one, it must be prime, so N equals \(p\) and will be added to result at the end, contradicting the assumption.

It is worth noting that if a prime table is precomputed, the time complexity can be reduced from \(O(\sqrt N)\) to \(O(\frac {\sqrt{N}} {\ln N})\). Refer to Sieve Methods for more information on precomputation.

Example problem: CF 1445C

Pollard Rho Algorithm

Introduction

The naive algorithm for finding a non-trivial divisor has complexity \(O(p)=O(\sqrt N)\), where \(p\) is the smallest prime factor of \(N\). The Pollard-Rho algorithm introduced below is a randomized algorithm that can find a non-trivial divisor in expected complexity \(O(\sqrt p)=O(N^{1/4})\) (Note: A non-trivial divisor is not necessarily a prime factor).

The core idea is: for a random self-mapping \(f: \mathbb Z_p \rightarrow \mathbb Z_p\), starting from any point \(x_1\), iteratively compute \(x_n = f(x_{n-1})\). The sequence will enter a cycle in expected time \(O(\sqrt p)\). If we can find \(x_i \equiv x_j \pmod p\), then \(p\) divides \(\gcd(|x_i-x_j|, N)\), and this greatest common divisor is a non-trivial factor of \(N\).

To understand why the expected time to enter a cycle is \(O(\sqrt p)\), consider the birthday paradox.

Birthday Paradox

Ignoring birth years (assuming each year has 365 days), what is the minimum number of people in a room such that the probability of at least two people sharing a birthday is \(50\%\)?

Solution: Assume a year has \(n\) days, and there are \(k\) people in the room, labeled with integers \(1, 2,\dots, k\). Assume each person's birthday is uniformly distributed among the \(n\) days, and birthdays are independent.

Let event \(A\) be that all \(k\) people have distinct birthdays. The probability of \(A\) is

\[ P(A)=\prod_{i=0}^{k-1}\frac{n-i}{n} \]

The probability that at least two people share a birthday is \(P(\overline A)=1-P(A)\). According to the problem, \(P(\overline A)\ge\frac{1}{2}\), so

\[ P(A)=\prod_{i=0}^{k-1}\frac{n-i}{n} \le \frac{1}{2} \]

Using the inequality \(1+x\le \mathrm{e}^x\), we have

\[ P(A) \le \prod_{i=1}^{k-1}\exp\left({-\frac{i}{n}}\right)=\exp \left({-\frac{k(k-1)}{2n}}\right) \]

Thus,

\[ \exp\left({-\dfrac{k(k-1)}{2n}}\right) \le \frac{1}{2}\implies P(A) \le \frac{1}{2} \]

Substituting \(n=365\), we find \(k\geq 23\). Therefore, at least \(23\) people are needed in a room for the probability of at least two sharing a birthday to reach \(50\%\). This mathematical fact is counterintuitive, hence called a paradox.

When \(k>56\) and \(n=365\), the probability of at least two people sharing a birthday exceeds \(99\%\)¹. For a year with \(n\) days, when there are \(\frac{1}{2}(\sqrt{8n\ln 2+1}+1)\approx \sqrt{2n\ln 2}\) people in the room, the probability of at least two sharing a birthday is about \(50\%\).

Similarly, it can be computed that the expected number of people needed to first obtain a repeated birthday when uniformly randomly selecting a sequence of birthdays is also \(O(\sqrt n)\). Let this number be \(X\); then

\[ E(X) = \sum_{x=1}^{n+1}P(X\ge x+1) = \sum_{x=0}^n\frac{n!}{(n-x)!n^x} = \sqrt{\frac{\pi n}{2}}-\frac13+o(1). \]

This suggests that if we can randomly select a sequence of numbers, the expected sample size to observe a repetition is also \(O(\sqrt n)\).

Finding a Divisor Using Greatest Common Divisor

Actually constructing a random sequence modulo \(p\) is not feasible because \(p\) is what we need to find. Instead, we generate a pseudo-random sequence \(\{x_i\}\) using \(f(x)=(x^2+c)\bmod N\): randomly choose \(x_1\), then let \(x_2=f(x_1),\ x_3=f(x_2),\ \dots,\ x_i=f(x_{i-1})\), where \(c\in[1,N)\) is a randomly chosen constant.

This chosen function is easy to compute and often generates a fairly random sequence. However, it is not completely random. For example, let \(n=50,\ c=6,\ x_1=1\). The data generated by \(f(x)\) is

\[ 1, 7, 5, 31, 17, 45, 31, 17, 45, 31,\dots \]

It can be observed that after \(x_4\), the data cycles among \(31, 17, 45\). If these numbers are arranged as shown in the figure below, the image resembles the Greek letter \(\rho\), hence the algorithm's name.

pollard-rho

More importantly, such a function indeed provides a self-mapping on \(\mathbb Z_p\). That is, it satisfies: if \(x\equiv y\pmod p\), then \(f(x)\equiv f(y)\pmod p\).

Proof

If \(x\equiv y\pmod p\), then \(x^2+c\equiv y^2+c\pmod p\). Note that \(f(x)=x^2+c-k_xN\), where \(k_x\) is an integer depending on \(x\), and \(p|N\). Thus, \(f(x)=x^2+c\pmod p\), so \(f(x)=f(y)\pmod p\).

As a pseudo-random self-mapping on \(\mathbb Z_p\) iteratively generates the sequence, \(\{x_n\bmod p\}\) will repeat in expected time \(O(\sqrt p)\). Once we observe such a repetition \(x_i\equiv x_j\pmod p\), we can obtain a non-trivial factor of \(N\) via \(\gcd(|x_i-x_j|,N)\). Note that since \(p\) is unknown, we cannot directly detect repetition. A simple detection method is to check whether \(\gcd(|x_i-x_j|,N)\) is strictly greater than one.

This algorithm does not always succeed because \(\gcd(|x_i-x_j|,N)\) might equal \(N\). That is, \(x_i\equiv x_j\pmod N\). In this case, when \(\{x_n\bmod p\}\) first repeats, \(\{x_n\}\) also repeats simultaneously. We do not obtain a non-trivial factor. Moreover, once \(\{x_n\}\) enters a cycle, further iterations are meaningless as they will only repeat the cycle. The algorithm should output failure and restart with a different \(c\) in \(f(x)\).

Based on the above analysis, theoretically, any function \(f(x)\) satisfying \(\forall x \equiv y \pmod p, f(x) \equiv f(y) \pmod p\) and ensuring sufficient pseudo-randomness (e.g., certain polynomial functions) can be used here. In practice, \(f(x)=x^2+c\ (c\neq 0,-2)\) is primarily used.²

Implementation

We need to implement an algorithm that can quickly detect whether \(\{x_n\bmod p\}\) has repeated during iteration. Viewing \(f\) as edges on a directed graph with vertices \(\mathbb Z_p\), we essentially need a cycle detection algorithm, but with equality checking replaced by checking whether \(\gcd(|x_i-x_j|,N)\) is greater than one.

Floyd Cycle Detection

Assume two people are racing: A is fast, B is slow. After some time, A will definitely catch up with B, and the difference in distances traveled by A and B when they meet is a multiple of the cycle length.

Let \(a=f(0), b=f(f(0))\). Each iteration updates \(a=f(a), b=f(f(b))\). Check whether \(a\) and \(b\) become equal during updates; if so, a cycle is detected.

Each time, let \(d=\gcd(|x_i-x_j|,N)\). If \(1< d< N\), return \(d\) directly. If \(d=N\), it means \(\{x_i\}\) has formed a cycle. When a cycle forms, we cannot proceed further; return \(N\) itself, and in subsequent attempts, adjust the random constant \(c\) and restart the factorization.

Pollard-Rho Algorithm Based on Floyd Cycle Detection

C++Python

ll Pollard_Rho(ll N) {
  if (N == 4) return 2;  // Special case for 4 because we skip two steps initially.
  ll c = rand() % (N - 1) + 1;
  ll t = f(0, c, N);
  ll r = f(f(0, c, N), c, N);
  while (t != r) {
    ll d = gcd(abs(t - r), N);
    if (d > 1) return d;
    t = f(t, c, N);
    r = f(f(r, c, N), c, N);
  }
  return N;
}

import random


def Pollard_Rho(N):
    if N == 4:
        return 2  # Special case for 4 because we skip two steps initially.
    c = random.randint(1, N - 1)
    t = f(0, c, N)
    r = f(f(0, c, N), c, N)
    while t != r:
        d = gcd(abs(t - r), N)
        if d > 1:
            return d
        t = f(t, c, N)
        r = f(f(r, c, N), c, N)
    return N

Brent Cycle Detection

Actually, Floyd's cycle detection algorithm can be improved in terms of constant factors. Brent's cycle detection starts with \(k=1\) and increments \(k\) each round. In the \(k\)-th round, A stays in place while B moves \(2^k\) steps forward. If B encounters A during this process, a cycle is detected. Otherwise, A teleports to B's position, and the next round begins.

It can be proven³ that the number of calls to \(f\) before detecting a cycle with Brent's method never exceeds that of Floyd's method. Tests in the original paper show that Brent's method reduces average time by about \(24\%\) compared to Floyd's method.

Multiplication Accumulation Optimization

Both Floyd and Brent cycle detection require \(O(\sqrt p)\) iterations. However, checking for a cycle using \(\gcd\) in each iteration slows down the algorithm. We can reduce the number of \(\gcd\) computations by accumulating products.

Simply put, if \(\gcd(a,N)>1\), then \(\gcd(ab\bmod N,N)=\gcd(ab,N)>1\) for any \(b\in\mathbb N_+\). That is, if we compute \(\gcd(\prod |x_i-x_j| \bmod N,N)>1\), then there must exist a pair \((x_i,x_j)\) such that \(\gcd(|x_i-x_j|,N)>1\). If the product becomes zero at any point, the factorization fails; exit and return \(N\) itself.

If we compute \(\gcd\) every \(k\) pairs, the algorithm complexity reduces to \(O(\sqrt p+k^{-1}\sqrt p\log N)\), where \(\log N\) is the cost of a single \(\gcd\) computation. Note that when \(k\) and \(\log N\) are roughly of the same order, we obtain an expected complexity of \(O(\sqrt p)\). In practice, \(k=128\) is commonly chosen.

Here is an implementation of Pollard-Rho using Brent cycle detection with multiplication accumulation optimization.

Implementation

C++Python

ll Pollard_Rho(ll x) {
  ll t = 0;
  ll c = rand() % (x - 1) + 1;
  ll s = t;
  int step = 0, goal = 1;
  ll val = 1;
  for (goal = 1;; goal <<= 1, s = t, val = 1) {
    for (step = 1; step <= goal; ++step) {
      t = f(t, c, x);
      val = val * abs(t - s) % x;
      // If val becomes 0, exit and restart factorization.
      if (!val) return x;
      if (step % 127 == 0) {
        ll d = gcd(val, x);
        if (d > 1) return d;
      }
    }
    ll d = gcd(val, x);
    if (d > 1) return d;
  }
}

from random import randint
from math import gcd


def Pollard_Rho(x):
    c = randint(1, x - 1)
    s = t = f(0, c, x)
    goal = val = 1
    while True:
        for step in range(1, goal + 1):
            t = f(t, c, x)
            val = val * abs(t - s) % x
            if val == 0:
                return x  # If val becomes 0, exit and restart factorization.
            if step % 127 == 0:
                d = gcd(val, x)
                if d > 1:
                    return d
        d = gcd(val, x)
        if d > 1:
            return d
        s = t
        goal <<= 1
        val = 1

Complexity

The expected number of iterations in Pollard-Rho is \(O(\sqrt p)\), where \(p\) is the smallest prime factor of \(N\). In implementations without multiplication accumulation optimization, whether using Floyd or Brent cycle detection, the expected complexity is \(O(\sqrt p\log N)\). With multiplication accumulation optimization, the expected complexity approximates \(O(\sqrt p)\).

It is worth mentioning that the above analysis is based on a completely random self-mapping function, but Pollard-Rho actually uses a pseudo-random function. Therefore, there is no strict complexity analysis, though it typically runs fast in practice.

Example: Finding the Largest Prime Factor of a Number

Example problem: P4718【Template】Pollard-Rho Algorithm

For a number \(n\), first use the Miller Rabin primality test to check if it is prime. If so, return directly. Otherwise, use Pollard-Rho to find a factor \(p\), divide \(n\) by \(p\), and recursively factorize \(n\) and \(p\). Use Miller Rabin to check for prime factors and update max_factor to find the largest prime factor. Due to the large data size of this problem, Floyd cycle detection is insufficient; the multiplication accumulation optimization method is used here.

Implementation

#include <algorithm>
#include <cstdlib>
#include <ctime>
#include <iostream>

using namespace std;
using ll = long long;
using ull = unsigned long long;

int t;
ll max_factor, n;

ll gcd(ll a, ll b) {
  if (b == 0) return a;
  return gcd(b, a % b);
}

ll bmul(ll a, ll b, ll m) {  // Fast multiplication
  ull c = (ull)a * (ull)b - (ull)((long double)a / m * b + 0.5L) * (ull)m;
  if (c < (ull)m) return c;
  return c + m;
}

ll qpow(ll x, ll p, ll mod) {  // Fast exponentiation
  ll ans = 1;
  while (p) {
    if (p & 1) ans = bmul(ans, x, mod);
    x = bmul(x, x, mod);
    p >>= 1;
  }
  return ans;
}

bool Miller_Rabin(ll p) {  // Primality test
  if (p < 2) return false;
  if (p == 2) return true;
  if (p == 3) return true;
  ll d = p - 1, r = 0;
  while (!(d & 1)) ++r, d >>= 1;  // Make d odd
  for (ll k = 0; k < 10; ++k) {
    ll a = rand() % (p - 2) + 2;
    ll x = qpow(a, d, p);
    if (x == 1 || x == p - 1) continue;
    for (int i = 0; i < r - 1; ++i) {
      x = bmul(x, x, p);
      if (x == p - 1) break;
    }
    if (x != p - 1) return false;
  }
  return true;
}

ll Pollard_Rho(ll x) {
  ll s = 0, t = 0;
  ll c = (ll)rand() % (x - 1) + 1;
  int step = 0, goal = 1;
  ll val = 1;
  for (goal = 1;; goal *= 2, s = t, val = 1) {  // Doubling optimization
    for (step = 1; step <= goal; ++step) {
      t = (bmul(t, t, x) + c) % x;
      val = bmul(val, abs(t - s), x);
      if ((step % 127) == 0) {
        ll d = gcd(val, x);
        if (d > 1) return d;
      }
    }
    ll d = gcd(val, x);
    if (d > 1) return d;
  }
}

void fac(ll x) {
  if (x <= max_factor || x < 2) return;
  if (Miller_Rabin(x)) {              // If x is prime
    max_factor = max(max_factor, x);  // Update answer
    return;
  }
  ll p = x;
  while (p >= x) p = Pollard_Rho(x);  // Use this algorithm
  while ((x % p) == 0) x /= p;
  fac(x), fac(p);  // Continue factoring x and p
}

int main() {
  cin >> t;
  while (t--) {
    srand((unsigned)time(NULL));
    max_factor = 0;
    cin >> n;
    fac(n);
    if (max_factor == n)  // The largest prime factor is itself
      cout << "Prime\n";
    else
      cout << max_factor << '\n';
  }
  return 0;
}

References and Links

https://en.wikipedia.org/wiki/Birthday_problem#Reverse_problem ↩
Menezes, Alfred J.; van Oorschot, Paul C.; Vanstone, Scott A. (2001). Handbook of Applied Cryptography. Section 3.11 and 3.12. ↩
Brent, R. P. (1980), An improved Monte Carlo factorization algorithm, BIT Numerical Mathematics, 20(2): 176–184, doi:10.1007/BF01933190 ↩

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用