Prime factorization

Introduction

Given a positive integer $N \in \mathbf{N}_{+}$, try to quickly find one of its non-trivial factors.

Consider the naive algorithm. Factors occur in pairs, so all factors of $N$ can be divided into $[2, \sqrt N]$ and $[\sqrt N+1,N)$. It is enough to enumerate the numbers in $[2, \sqrt N]$ and use division to find at least two factors. This method has time complexity $O(\sqrt N)$.

When $N\ge10^{18}$, the running time of this algorithm is unacceptable, so we seek a better algorithm. One idea is to guess randomly whether a number is a factor of $N$. A lucky guess can find the answer in $O(1)$ time, but for $N\ge10^{18}$ the probability of a successful guess is $\frac{1}{10^{18}}$, and the expected number of guesses is $10^{18}$. Guessing within $[2,\sqrt N]$ improves the success probability somewhat. We want to optimize this guessing process.

Naive algorithm

The simplest algorithm enumerates $[2, \sqrt N]$.

C++

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) {  // The operations have left a prime factor in N
    result.push_back(N);
  }
  return result;
}

Python

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:  # The operations have left a prime factor in N
        result.append(N)
    return result

We can prove that every element in result is a prime factor of N.

It is worth noting that if a prime table has already been prepared, the time complexity decreases from $O(\sqrt N)$ to $O(\frac {\sqrt{N}} {\ln N})$. See sieves for more information about precomputing such tables.

Example: CF 1445C

Pollard-Rho algorithm

Introduction

The brute-force complexity of finding a non-trivial factor is $O(p)=O(\sqrt N)$, where $p$ is the smallest prime factor of $N$. The Pollard-Rho algorithm introduced below is randomized and finds a non-trivial factor in expected time $O(\sqrt p)=O(N^{1/4})$ (Note: a non-trivial factor need not be prime).

The key idea is to take a random self-map $f: \mathbb Z_p \rightarrow \mathbb Z_p$, start from any point $x_1$, and iterate $x_n = f(x_{n-1})$. The sequence enters a cycle in expected time $O(\sqrt p)$. If we find $x_i \equiv x_j \pmod p$, then $p$ divides $\gcd(|x_i-x_j|, N)$, and this gcd is a non-trivial factor of $N$.

The birthday paradox provides intuition for the expected time $O(\sqrt p)$ to enter a cycle.

Birthday paradox

Ignoring birth years (and assuming 365 days in every year), how many people must a room contain for the probability that two of them share a birthday to reach $50\%$?

Solution: suppose a year has $n$ days and the room has $k$ people, numbered $1, 2,\dots, k$. Assume each birthday is uniformly distributed among the $n$ days and that birthdays are mutually independent.

Let event $k$ $A$ $A$ be that the people have distinct birthdays. Its probability is

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

Using the inequality $1+x\le \mathrm{e}^x$, we obtain

Therefore

Substituting $n=365$ gives $k\geq 23$. Thus a room needs at least $23$ people for the probability of a shared birthday to reach $50\%$. This mathematical fact is highly counterintuitive, so it is called a paradox.

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

Similarly, if birthdays are sampled independently and uniformly, the expected number of samples needed to obtain the first repeated birthday is also $O(\sqrt n)$. Let this number be $X$. Then

This suggests that if we can sample a sequence of numbers at random, the expected sample size before a repeated number appears is also $O(\sqrt n)$.

Finding a factor using the greatest common divisor

Constructing a random sequence modulo $p$ is impractical because $p$ is precisely what we seek. Instead, generate a pseudorandom sequence $f(x)=(x^2+c)\bmod N$ using $\{x_i\}$: choose $x_1$ randomly and set $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.

The chosen function is easy to compute and often generates a sufficiently random-looking sequence, but it is not truly random. For example, let $n=50,\ c=6,\ x_1=1$ $f(x)$, , and ; the values generated by are

After $x_4$, the values cycle among $31,17,45$.Arranging these values as in the figure produces a shape resembling $\rho$, which gives the algorithm its name.

More importantly, this function really is a self-map on $\mathbb Z_p$: if $x\equiv y\pmod p$, then $f(x)\equiv f(y)\pmod p$.

As a sequence obtained by repeatedly iterating a pseudorandom self-map on $\mathbb Z_p$, $\{x_n\bmod p\}$ repeats within expected time $O(\sqrt p)$. Once we observe $x_i\equiv x_j\pmod p$, we can obtain a non-trivial factor of $\gcd(|x_i-x_j|,N)$ from $N$. Since $p$ is unknown, we cannot detect the repetition directly; a simple test is that $\gcd(|x_i-x_j|,N)$ is strictly greater than one.

The algorithm does not always succeed, because $\gcd(|x_i-x_j|,N)$ may equal $N$, meaning $x_i\equiv x_j\pmod N$. In that case, when $\{x_n\bmod p\}$ first repeats, $\{x_n\}$ also repeats. No non-trivial factor is obtained. Once $\{x_n\}$ enters a cycle, further iteration is pointless because the same cycle will repeat. The algorithm should report factorization failure and retry with a different choice of $f(x)$ in $c$.

According to the analysis above, any function $\forall x \equiv y \pmod p, f(x) \equiv f(y) \pmod p$ satisfying $f(x)$ and providing sufficient pseudorandomness (for example, certain polynomial functions) can be used here. In practice, $f(x)=x^2+c\ (c\neq 0,-2)$ is used most often.²

Implementation

The implementation must quickly determine during iteration whether $\{x_n\bmod p\}$ has repeated. Viewing $f$ as edges of a directed graph whose vertices are $\mathbb Z_p$, we need a cycle-detection algorithm, replacing equality testing with the test whether $\gcd(|x_i-x_j|,N)$ is greater than one.

Floyd cycle detection

Suppose two people race around a track, with A faster than B. After enough time A must meet B, and at the meeting point the distance run by A minus the distance run by B is a multiple of the track length.

Let $a=f(0),b=f(f(0))$ $a=f(a),b=f(f(b))$ $a$ $b$ and . At each step update and . If and become equal during the updates, a cycle has been found.

At each step set $d=\gcd(|x_i-x_j|,N)$ and test whether $1< d< N$; if so, return $d$. If $d=N$, then $\{x_i\}$ has formed a cycle, so stop and return $N$ itself, adjust the random constant $c$, and try factorization again.

ll Pollard_Rho(ll N) {
  if (N == 4) return 2;  // We initially advance by two steps, so 4 needs a special case
  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  # We initially advance by two steps, so 4 needs a special case
    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

The Floyd cycle-detection algorithm admits a constant-factor improvement. Brent's algorithm increases $k=1$ from $k$. In round $k$, A stays still while B moves forward $2^k$ steps. If B meets A, the cycle has been found; otherwise A jumps to B's position and the next round begins.

It can be proved³ that the number of calls to $f$ before finding a cycle is never greater than for Floyd detection. Tests in the original paper show that Brent detection reduces the average time by $24\%$ compared with Floyd detection.

Doubling optimization

Both Floyd and Brent cycle detection use $O(\sqrt p)$ $\gcd$ $\gcd$ iterations. However, computing a at every iteration slows the algorithm. Multiplicative accumulation can reduce the number of gcd computations.

In short, if $\gcd(a,N)>1$, then $\gcd(ab\bmod N,N)=\gcd(ab,N)>1$ for every $b\in\mathbb N_+$. Thus if $\gcd(\prod |x_i-x_j| \bmod N,N)>1$, some pair $(x_i,x_j)$ necessarily satisfies $\gcd(|x_i-x_j|,N)>1$. If the product becomes zero, factorization has failed; stop and return $N$ itself.

If a gcd is computed once every $k$ $\gcd$ $O(\sqrt p+k^{-1}\sqrt p\log N)$ $\log N$ $\gcd$ $k$ $\log N$ $O(\sqrt p)$ $k=128$ pairs, the complexity falls to , where is the cost of one gcd computation. When and have roughly the same order, the expected complexity is . Implementations usually choose .

Here is an implementation of Pollard-Rho using Brent cycle detection and the doubling optimization.

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 is zero, stop and factor again
      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 is zero, stop and factor again
            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$. Whether Floyd or Brent detection is used, the expected complexity without doubling is $O(\sqrt p\log N)$; with doubling, it is approximately $O(\sqrt p)$.

It is worth noting that the preceding analysis assumes a completely random self-map, whereas Pollard-Rho uses a pseudorandom function. Therefore it has no strict complexity analysis, although it is usually fast in practice.

Example: finding the largest prime factor of a number

Example: P4718 (https://www.luogu.com.cn/problem/P4718) Pollard-Rho algorithm

For a number $n$, use the Miller-Rabin algorithm to test primality and return immediately if it is prime. Otherwise, use Pollard-Rho to find a factor $p$, divide $n$ by $p$, and recursively factor $n$ and $p$. Use Miller-Rabin to identify prime factors and update max_factor to obtain the largest one. Because the data in this problem is very large, Floyd detection is insufficient, so the doubling optimization is used.

#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) {  // Test primality
  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 the 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 the number itself
      cout << "Prime\n";
    else
      cout << max_factor << '\n';
  }
  return 0;
}

References and links