Divisors, GCD and Primes
The test everyone writes
Is 2,147,483,647 prime? It's the largest 32-bit integer, and it is. Here's the test most people write: try every divisor up to the square root, stopping once i × i passes n:
#include <bits/stdc++.h>
using namespace std;
bool isPrime(int n) { // the test everyone writes first
if (n < 2) return false;
for (int i = 2; i * i <= n; i++)
if (n % i == 0) return false;
return true;
}
int main() {
printf("%s %s\n", isPrime(97) ? "true" : "false", isPrime(2147483647) ? "true" : "false");
}
// true falsedef is_prime(n): # the test everyone writes first
if n < 2:
return False
i = 2
while i * i <= n:
if n % i == 0:
return False
i += 1
return True
print(is_prime(97), is_prime(2147483647))
# True Truepublic class Main {
static boolean isPrime(int n) { // the test everyone writes first
if (n < 2) return false;
for (int i = 2; i * i <= n; i++)
if (n % i == 0) return false;
return true;
}
public static void main(String[] args) {
System.out.println(isPrime(97) + " " + isPrime(2147483647));
}
}
// true falsefunction isPrime(n) { // the test everyone writes first
if (n < 2) return false;
for (let i = 2; i * i <= n; i++)
if (n % i === 0) return false;
return true;
}
console.log(isPrime(97), isPrime(2147483647));
// true truePython and JavaScript said yes. Java said no, after 10 seconds, and C++ said no at once. In 32 bits, i × i passes the largest int at i = 46,341 and wraps around to a negative number, so i * i <= n never turns false: no int is bigger than 2,147,483,647. Java ran on until i reached n itself and found that n divides n. In C++ the overflow is undefined behaviour, and the compiled program printed false without the wait. Python's ints don't overflow, and JavaScript's numbers are exact far past this size.
The idea under the test is right, and it's the idea of this lesson: divisors come in pairs, and one of each pair is at most √n. That bound finds every divisor of a number and tests it for primality; Euclid's algorithm finds the greatest common divisor of two numbers in a few dozen steps; and the sieve finds every prime below n at once.
Every divisor in √n steps
If d divides n, so does n / d, and the two sit on either side of √n: d × (n / d) = n, so both can't be bigger than √n. Checking d from 1 up to √n and keeping both halves of each pair finds every divisor. The test d <= n / d says the same as d * d <= n, and can't overflow:
#include <bits/stdc++.h>
using namespace std;
vector<long long> divisors(long long n) { // every divisor of n >= 1, in order
vector<long long> small, large;
for (long long d = 1; d <= n / d; d++) { // d <= n / d is d * d <= n, without the overflow
if (n % d != 0) continue;
small.push_back(d);
if (d != n / d) large.push_back(n / d); // a square's root pairs with itself: once
}
small.insert(small.end(), large.rbegin(), large.rend());
return small;
}
int main() {
for (long long n : {36LL, 97LL}) {
printf("%lld:", n);
for (long long d : divisors(n)) printf(" %lld", d);
printf("\n");
}
printf("720720: %zu divisors, d up to %lld\n", divisors(720720).size(), (long long) sqrt(720720.0));
}
// 36: 1 2 3 4 6 9 12 18 36
// 97: 1 97
// 720720: 240 divisors, d up to 848import math
def divisors(n): # every divisor of n >= 1, in order
small, large = [], []
d = 1
while d <= n // d: # d <= n / d is d * d <= n, without the overflow
if n % d == 0:
small.append(d)
if d != n // d: # a square's root pairs with itself: once
large.append(n // d)
d += 1
return small + large[::-1]
for n in (36, 97):
print(f"{n}:", " ".join(map(str, divisors(n))))
print(f"720720: {len(divisors(720720))} divisors, d up to {math.isqrt(720720)}")
# 36: 1 2 3 4 6 9 12 18 36
# 97: 1 97
# 720720: 240 divisors, d up to 848import java.util.*;
public class Main {
static List<Long> divisors(long n) { // every divisor of n >= 1, in order
List<Long> small = new ArrayList<>(), large = new ArrayList<>();
for (long d = 1; d <= n / d; d++) { // d <= n / d is d * d <= n, without the overflow
if (n % d != 0) continue;
small.add(d);
if (d != n / d) large.add(n / d); // a square's root pairs with itself: once
}
Collections.reverse(large);
small.addAll(large);
return small;
}
public static void main(String[] args) {
for (long n : new long[]{36, 97}) {
StringBuilder line = new StringBuilder(n + ":");
for (long d : divisors(n)) line.append(" ").append(d);
System.out.println(line);
}
System.out.println("720720: " + divisors(720720).size() + " divisors, d up to " + (long) Math.sqrt(720720));
}
}
// 36: 1 2 3 4 6 9 12 18 36
// 97: 1 97
// 720720: 240 divisors, d up to 848function divisors(n) { // every divisor of n >= 1, in order
const small = [], large = [];
for (let d = 1; d <= Math.floor(n / d); d++) { // d <= n / d is d * d <= n, without the overflow
if (n % d !== 0) continue;
small.push(d);
if (d !== n / d) large.push(n / d); // a square's root pairs with itself: once
}
return small.concat(large.reverse());
}
for (const n of [36, 97]) console.log(`${n}: ${divisors(n).join(" ")}`);
console.log(`720720: ${divisors(720720).length} divisors, d up to ${Math.floor(Math.sqrt(720720))}`);
// 36: 1 2 3 4 6 9 12 18 36
// 97: 1 97
// 720720: 240 divisors, d up to 84836 has nine divisors: four pairs and 6, the square root, which pairs with itself and is kept once. 97 has only 1 and itself, so it's prime. 720,720 has 240 divisors, found by trying d up to 848, where trying every number up to n would take 720,720 tries. O(√n) time.
The small halves come out in increasing order and the large halves in decreasing order, so the large ones are put after the small ones reversed, and the list is sorted with no sort.
The greatest common divisor: Euclid
The greatest common divisor of a and b, gcd(a, b), is the largest number that divides both: gcd(12, 18) = 6. Listing the divisors of both would work; Euclid's algorithm is far faster. It replaces the pair (a, b) by (b, a mod b) until the second number is 0, and the first is the answer.
Why it's right. Any number that divides a and b also divides a − q × b, which is a mod b. Any number that divides b and a mod b also divides q × b + a mod b, which is a. So (a, b) and (b, a mod b) have exactly the same common divisors, and the same greatest one. And gcd(a, 0) = a, because every number divides 0.
#include <bits/stdc++.h>
using namespace std;
long long gcd(long long a, long long b) { // Euclid: (a, b) becomes (b, a mod b) until b is 0
while (b != 0) {
long long r = a % b;
a = b;
b = r;
}
return a;
}
int steps(long long a, long long b) { // how many times the loop above runs
int k = 0;
for (; b != 0; k++) {
long long r = a % b;
a = b;
b = r;
}
return k;
}
int main() {
printf("gcd(1071, 462) = %lld in %d steps\n", gcd(1071, 462), steps(1071, 462));
printf("gcd(0, 5) = %lld\n", gcd(0, 5));
printf("gcd(1836311903, 1134903170) = %lld in %d steps\n", gcd(1836311903, 1134903170),
steps(1836311903, 1134903170));
}
// gcd(1071, 462) = 21 in 3 steps
// gcd(0, 5) = 5
// gcd(1836311903, 1134903170) = 1 in 44 stepsdef gcd(a, b): # Euclid: (a, b) becomes (b, a mod b) until b is 0
while b != 0:
a, b = b, a % b
return a
def steps(a, b): # how many times the loop above runs
k = 0
while b != 0:
a, b = b, a % b
k += 1
return k
print(f"gcd(1071, 462) = {gcd(1071, 462)} in {steps(1071, 462)} steps")
print(f"gcd(0, 5) = {gcd(0, 5)}")
print(f"gcd(1836311903, 1134903170) = {gcd(1836311903, 1134903170)} in {steps(1836311903, 1134903170)} steps")
# gcd(1071, 462) = 21 in 3 steps
# gcd(0, 5) = 5
# gcd(1836311903, 1134903170) = 1 in 44 stepspublic class Main {
static long gcd(long a, long b) { // Euclid: (a, b) becomes (b, a mod b) until b is 0
while (b != 0) {
long r = a % b;
a = b;
b = r;
}
return a;
}
static int steps(long a, long b) { // how many times the loop above runs
int k = 0;
for (; b != 0; k++) {
long r = a % b;
a = b;
b = r;
}
return k;
}
public static void main(String[] args) {
System.out.println("gcd(1071, 462) = " + gcd(1071, 462) + " in " + steps(1071, 462) + " steps");
System.out.println("gcd(0, 5) = " + gcd(0, 5));
System.out.println("gcd(1836311903, 1134903170) = " + gcd(1836311903L, 1134903170L) + " in "
+ steps(1836311903L, 1134903170L) + " steps");
}
}
// gcd(1071, 462) = 21 in 3 steps
// gcd(0, 5) = 5
// gcd(1836311903, 1134903170) = 1 in 44 stepsfunction gcd(a, b) { // Euclid: (a, b) becomes (b, a mod b) until b is 0
while (b !== 0) [a, b] = [b, a % b];
return a;
}
function steps(a, b) { // how many times the loop above runs
let k = 0;
for (; b !== 0; k++) [a, b] = [b, a % b];
return k;
}
console.log(`gcd(1071, 462) = ${gcd(1071, 462)} in ${steps(1071, 462)} steps`);
console.log(`gcd(0, 5) = ${gcd(0, 5)}`);
console.log(`gcd(1836311903, 1134903170) = ${gcd(1836311903, 1134903170)} in ${steps(1836311903, 1134903170)} steps`);
// gcd(1071, 462) = 21 in 3 steps
// gcd(0, 5) = 5
// gcd(1836311903, 1134903170) = 1 in 44 steps1,071 and 462 took 3 steps: 1,071 mod 462 = 147, 462 mod 147 = 21, 147 mod 21 = 0, so the answer is 21. Every two steps at least halve the smaller number, so the steps are O(log min(a, b)). The worst case is two neighbouring Fibonacci numbers, whose remainders shrink as slowly as possible: the largest pair below 2³¹, 1,836,311,903 and 1,134,903,170, took 44 steps. For more than two numbers, fold: gcd(a, b, c) = gcd(gcd(a, b), c), and a running gcd that reaches 1 can stop there.
The least common multiple is the smallest number both divide: lcm(4, 6) = 12. It's a × b / gcd(a, b), but written that way a × b can overflow before the division brings it back down. gcd(a, b) divides a exactly, so divide first:
#include <bits/stdc++.h>
using namespace std;
long long gcd(long long a, long long b) { return b == 0 ? a : gcd(b, a % b); }
long long lcm(long long a, long long b) { // divide first: gcd divides a exactly
return a / gcd(a, b) * b;
}
int main() {
printf("%lld %lld %lld\n", lcm(4, 6), lcm(21, 6), lcm(60000, 90000));
}
// 12 42 180000def gcd(a, b):
return a if b == 0 else gcd(b, a % b)
def lcm(a, b): # divide first: gcd divides a exactly
return a // gcd(a, b) * b
print(lcm(4, 6), lcm(21, 6), lcm(60000, 90000))
# 12 42 180000public class Main {
static long gcd(long a, long b) { return b == 0 ? a : gcd(b, a % b); }
static long lcm(long a, long b) { // divide first: gcd divides a exactly
return a / gcd(a, b) * b;
}
public static void main(String[] args) {
System.out.println(lcm(4, 6) + " " + lcm(21, 6) + " " + lcm(60000, 90000));
}
}
// 12 42 180000const gcd = (a, b) => (b === 0 ? a : gcd(b, a % b));
function lcm(a, b) { // divide first: gcd divides a exactly
return a / gcd(a, b) * b;
}
console.log(lcm(4, 6), lcm(21, 6), lcm(60000, 90000));
// 12 42 180000lcm(60,000, 90,000) = 180,000 fits easily in 32 bits, but 60,000 × 90,000 = 5,400,000,000 doesn't (Mistakes that cost an interview shows what multiplying first prints).
The GCD of two strings. A string x divides s if s is x repeated, like "AB" divides "ABABAB". Greatest Common Divisor of Strings asks for the longest x that divides both s and t. If some x builds both, then s + t and t + s are the same string, x repeated; and if s + t equals t + s, both are built from the prefix whose length is gcd(len(s), len(t)). So one comparison decides whether there is an answer, and Euclid gives its length:
#include <bits/stdc++.h>
using namespace std;
string gcdOfStrings(const string& s, const string& t) {
if (s + t != t + s) return ""; // no common building block
return s.substr(0, __gcd(s.size(), t.size())); // the block is gcd(lengths) long
}
int main() {
printf("[%s] [%s] [%s]\n", gcdOfStrings("ABCABC", "ABC").c_str(), gcdOfStrings("ABABAB", "ABAB").c_str(),
gcdOfStrings("LEET", "CODE").c_str());
}
// [ABC] [AB] []import math
def gcd_of_strings(s, t):
if s + t != t + s:
return "" # no common building block
return s[:math.gcd(len(s), len(t))] # the block is gcd(lengths) long
print(f"[{gcd_of_strings('ABCABC', 'ABC')}] [{gcd_of_strings('ABABAB', 'ABAB')}] [{gcd_of_strings('LEET', 'CODE')}]")
# [ABC] [AB] []public class Main {
static int gcd(int a, int b) { return b == 0 ? a : gcd(b, a % b); }
static String gcdOfStrings(String s, String t) {
if (!(s + t).equals(t + s)) return ""; // no common building block
return s.substring(0, gcd(s.length(), t.length())); // the block is gcd(lengths) long
}
public static void main(String[] args) {
System.out.println("[" + gcdOfStrings("ABCABC", "ABC") + "] [" + gcdOfStrings("ABABAB", "ABAB") + "] ["
+ gcdOfStrings("LEET", "CODE") + "]");
}
}
// [ABC] [AB] []const gcd = (a, b) => (b === 0 ? a : gcd(b, a % b));
function gcdOfStrings(s, t) {
if (s + t !== t + s) return ""; // no common building block
return s.slice(0, gcd(s.length, t.length)); // the block is gcd(lengths) long
}
console.log(`[${gcdOfStrings("ABCABC", "ABC")}] [${gcdOfStrings("ABABAB", "ABAB")}] [${gcdOfStrings("LEET", "CODE")}]`);
// [ABC] [AB] []"ABCABC" and "ABC" give "ABC"; "ABABAB" and "ABAB" give "AB", the gcd of 6 and 4 being 2; "LEET" and "CODE" give the empty string, because "LEETCODE" isn't "CODELEET".
Is n prime?
n is prime if it's at least 2 and nothing from 2 to √n divides it: a divisor above √n would pair with one below. The same loop as the first one, with the bound that can't overflow:
#include <bits/stdc++.h>
using namespace std;
bool isPrime(int n) {
if (n < 2) return false; // 0 and 1 aren't prime
for (int i = 2; i <= n / i; i++) // i <= n / i: the same as i * i <= n, and it can't overflow
if (n % i == 0) return false;
return true;
}
int main() {
for (int n : {1, 2, 97, 2147483647, 2147483646})
printf("%d %s\n", n, isPrime(n) ? "prime" : "not prime");
}
// 1 not prime
// 2 prime
// 97 prime
// 2147483647 prime
// 2147483646 not primedef is_prime(n):
if n < 2:
return False # 0 and 1 aren't prime
i = 2
while i <= n // i: # i <= n / i: the same as i * i <= n, and it can't overflow
if n % i == 0:
return False
i += 1
return True
for n in (1, 2, 97, 2147483647, 2147483646):
print(n, "prime" if is_prime(n) else "not prime")
# 1 not prime
# 2 prime
# 97 prime
# 2147483647 prime
# 2147483646 not primepublic class Main {
static boolean isPrime(int n) {
if (n < 2) return false; // 0 and 1 aren't prime
for (int i = 2; i <= n / i; i++) // i <= n / i: the same as i * i <= n, and it can't overflow
if (n % i == 0) return false;
return true;
}
public static void main(String[] args) {
for (int n : new int[]{1, 2, 97, 2147483647, 2147483646})
System.out.println(n + " " + (isPrime(n) ? "prime" : "not prime"));
}
}
// 1 not prime
// 2 prime
// 97 prime
// 2147483647 prime
// 2147483646 not primefunction isPrime(n) {
if (n < 2) return false; // 0 and 1 aren't prime
for (let i = 2; i <= Math.floor(n / i); i++) // i <= n / i: the same as i * i <= n
if (n % i === 0) return false;
return true;
}
for (const n of [1, 2, 97, 2147483647, 2147483646])
console.log(n, isPrime(n) ? "prime" : "not prime");
// 1 not prime
// 2 prime
// 97 prime
// 2147483647 prime
// 2147483646 not prime2,147,483,647 is prime in all four now, after 46,339 tries. 0 and 1 aren't prime, and 2 is the only even one.
Every prime above 3 is one more or one less than a multiple of 6: 6k, 6k + 2 and 6k + 4 are even, and 6k + 3 is a multiple of 3. So after 2 and 3, only 6k − 1 and 6k + 1 need trying:
#include <bits/stdc++.h>
using namespace std;
bool isPrime(int n) { // 2 and 3, then only 6k - 1 and 6k + 1
if (n < 4) return n >= 2; // 2 and 3 are prime; 0 and 1 aren't
if (n % 2 == 0 || n % 3 == 0) return false;
for (int i = 5; i <= n / i; i += 6) // i is 6k - 1, and i + 2 is 6k + 1
if (n % i == 0 || n % (i + 2) == 0) return false;
return true;
}
int main() {
for (int n : {1, 2, 97, 2147483647, 2147483646})
printf("%d %s\n", n, isPrime(n) ? "prime" : "not prime");
}
// 1 not prime
// 2 prime
// 97 prime
// 2147483647 prime
// 2147483646 not primedef is_prime(n): # 2 and 3, then only 6k - 1 and 6k + 1
if n < 4:
return n >= 2 # 2 and 3 are prime; 0 and 1 aren't
if n % 2 == 0 or n % 3 == 0:
return False
i = 5
while i <= n // i: # i is 6k - 1, and i + 2 is 6k + 1
if n % i == 0 or n % (i + 2) == 0:
return False
i += 6
return True
for n in (1, 2, 97, 2147483647, 2147483646):
print(n, "prime" if is_prime(n) else "not prime")
# 1 not prime
# 2 prime
# 97 prime
# 2147483647 prime
# 2147483646 not primepublic class Main {
static boolean isPrime(int n) { // 2 and 3, then only 6k - 1 and 6k + 1
if (n < 4) return n >= 2; // 2 and 3 are prime; 0 and 1 aren't
if (n % 2 == 0 || n % 3 == 0) return false;
for (int i = 5; i <= n / i; i += 6) // i is 6k - 1, and i + 2 is 6k + 1
if (n % i == 0 || n % (i + 2) == 0) return false;
return true;
}
public static void main(String[] args) {
for (int n : new int[]{1, 2, 97, 2147483647, 2147483646})
System.out.println(n + " " + (isPrime(n) ? "prime" : "not prime"));
}
}
// 1 not prime
// 2 prime
// 97 prime
// 2147483647 prime
// 2147483646 not primefunction isPrime(n) { // 2 and 3, then only 6k - 1 and 6k + 1
if (n < 4) return n >= 2; // 2 and 3 are prime; 0 and 1 aren't
if (n % 2 === 0 || n % 3 === 0) return false;
for (let i = 5; i <= Math.floor(n / i); i += 6) // i is 6k - 1, and i + 2 is 6k + 1
if (n % i === 0 || n % (i + 2) === 0) return false;
return true;
}
for (const n of [1, 2, 97, 2147483647, 2147483646])
console.log(n, isPrime(n) ? "prime" : "not prime");
// 1 not prime
// 2 prime
// 97 prime
// 2147483647 prime
// 2147483646 not primeThe same answers, with 15,448 divisions for 2,147,483,647 instead of 46,339, about a third. Still O(√n) per number: fine for one number, too slow for every number below a million, which is the sieve's job.
Every prime below n: the sieve
To find every prime below n, the Sieve of Eratosthenes starts with every number from 2 marked as a candidate. The smallest unmarked number is prime; cross out its multiples, then take the next unmarked number, and repeat.
Each prime p starts crossing at p × p, not at 2p. A smaller multiple, k × p with k below p, has a prime factor smaller than p (some prime factor of k), so it was crossed out already. And once p × p reaches n, p crosses out nothing, so the loop can stop: every composite below n has a prime factor no bigger than its square root, and has been crossed out by then.
#include <bits/stdc++.h>
using namespace std;
int countPrimes(int n) { // how many primes are below n
vector<bool> composite(max(n, 2), false);
for (long long p = 2; p * p < n; p++) { // a prime with p * p >= n crosses nothing out
if (composite[p]) continue;
for (long long m = p * p; m < n; m += p) // smaller multiples have a smaller prime factor
composite[m] = true;
}
int count = 0;
for (int k = 2; k < n; k++) count += !composite[k];
return count;
}
int main() {
printf("%d %d %d %d\n", countPrimes(10), countPrimes(100), countPrimes(1000000), countPrimes(5000000));
}
// 4 25 78498 348513def count_primes(n): # how many primes are below n
composite = [False] * max(n, 2)
p = 2
while p * p < n: # a prime with p * p >= n crosses nothing out
if not composite[p]:
for m in range(p * p, n, p): # smaller multiples have a smaller prime factor
composite[m] = True
p += 1
return sum(1 for k in range(2, n) if not composite[k])
print(count_primes(10), count_primes(100), count_primes(1000000), count_primes(5000000))
# 4 25 78498 348513public class Main {
static int countPrimes(int n) { // how many primes are below n
boolean[] composite = new boolean[Math.max(n, 2)];
for (long p = 2; p * p < n; p++) { // a prime with p * p >= n crosses nothing out
if (composite[(int) p]) continue;
for (long m = p * p; m < n; m += p) // smaller multiples have a smaller prime factor
composite[(int) m] = true;
}
int count = 0;
for (int k = 2; k < n; k++) if (!composite[k]) count++;
return count;
}
public static void main(String[] args) {
System.out.println(countPrimes(10) + " " + countPrimes(100) + " " + countPrimes(1000000) + " "
+ countPrimes(5000000));
}
}
// 4 25 78498 348513function countPrimes(n) { // how many primes are below n
const composite = new Uint8Array(Math.max(n, 2));
for (let p = 2; p * p < n; p++) { // a prime with p * p >= n crosses nothing out
if (composite[p]) continue;
for (let m = p * p; m < n; m += p) composite[m] = 1; // smaller multiples have a smaller prime factor
}
let count = 0;
for (let k = 2; k < n; k++) if (!composite[k]) count++;
return count;
}
console.log(countPrimes(10), countPrimes(100), countPrimes(1000000), countPrimes(5000000));
// 4 25 78498 348513There are 4 primes below 10, 25 below 100, 78,498 below a million and 348,513 below five million, the largest n that Count Primes allows. Each prime p crosses out about n / p numbers, and the sum over the primes is O(n log log n), close to linear. The price is memory: one cell per number.
Against trial division of every number, timed with the two functions above:
#include <bits/stdc++.h>
using namespace std;
bool isPrime(int n) {
if (n < 2) return false; // 0 and 1 aren't prime
for (int i = 2; i <= n / i; i++) // i <= n / i: the same as i * i <= n, and it can't overflow
if (n % i == 0) return false;
return true;
}
int countPrimes(int n) { // how many primes are below n
vector<bool> composite(max(n, 2), false);
for (long long p = 2; p * p < n; p++) { // a prime with p * p >= n crosses nothing out
if (composite[p]) continue;
for (long long m = p * p; m < n; m += p) // smaller multiples have a smaller prime factor
composite[m] = true;
}
int count = 0;
for (int k = 2; k < n; k++) count += !composite[k];
return count;
}
int main() {
int n = 10000000;
auto t = chrono::steady_clock::now();
int trial = 0;
for (int k = 2; k < n; k++) trial += isPrime(k);
double a = chrono::duration<double>(chrono::steady_clock::now() - t).count();
t = chrono::steady_clock::now();
int sieve = countPrimes(n);
double b = chrono::duration<double>(chrono::steady_clock::now() - t).count();
printf("below %d: trial division %d primes in %.3fs, the sieve %d in %.3fs\n", n, trial, a, sieve, b);
}
// below 10000000: trial division 664579 primes in 5.175s, the sieve 664579 in 0.038sdef is_prime(n):
if n < 2:
return False # 0 and 1 aren't prime
i = 2
while i <= n // i: # i <= n / i: the same as i * i <= n, and it can't overflow
if n % i == 0:
return False
i += 1
return True
def count_primes(n): # how many primes are below n
composite = [False] * max(n, 2)
p = 2
while p * p < n: # a prime with p * p >= n crosses nothing out
if not composite[p]:
for m in range(p * p, n, p): # smaller multiples have a smaller prime factor
composite[m] = True
p += 1
return sum(1 for k in range(2, n) if not composite[k])
import time
n = 1000000
start = time.perf_counter()
trial = sum(1 for k in range(2, n) if is_prime(k))
a = time.perf_counter() - start
start = time.perf_counter()
sieve = count_primes(n)
b = time.perf_counter() - start
print(f"below {n}: trial division {trial} primes in {a:.3f}s, the sieve {sieve} in {b:.3f}s")
# below 1000000: trial division 78498 primes in 8.203s, the sieve 78498 in 0.182spublic class Main {
static boolean isPrime(int n) {
if (n < 2) return false; // 0 and 1 aren't prime
for (int i = 2; i <= n / i; i++) // i <= n / i: the same as i * i <= n, and it can't overflow
if (n % i == 0) return false;
return true;
}
static int countPrimes(int n) { // how many primes are below n
boolean[] composite = new boolean[Math.max(n, 2)];
for (long p = 2; p * p < n; p++) { // a prime with p * p >= n crosses nothing out
if (composite[(int) p]) continue;
for (long m = p * p; m < n; m += p) // smaller multiples have a smaller prime factor
composite[(int) m] = true;
}
int count = 0;
for (int k = 2; k < n; k++) if (!composite[k]) count++;
return count;
}
public static void main(String[] args) {
int n = 10000000;
long t = System.nanoTime();
int trial = 0;
for (int k = 2; k < n; k++) if (isPrime(k)) trial++;
double a = (System.nanoTime() - t) / 1e9;
t = System.nanoTime();
int sieve = countPrimes(n);
double b = (System.nanoTime() - t) / 1e9;
System.out.printf("below %d: trial division %d primes in %.3fs, the sieve %d in %.3fs%n", n, trial, a, sieve, b);
}
}
// below 10000000: trial division 664579 primes in 10.731s, the sieve 664579 in 0.076sfunction isPrime(n) {
if (n < 2) return false; // 0 and 1 aren't prime
for (let i = 2; i <= Math.floor(n / i); i++) // i <= n / i: the same as i * i <= n
if (n % i === 0) return false;
return true;
}
function countPrimes(n) { // how many primes are below n
const composite = new Uint8Array(Math.max(n, 2));
for (let p = 2; p * p < n; p++) { // a prime with p * p >= n crosses nothing out
if (composite[p]) continue;
for (let m = p * p; m < n; m += p) composite[m] = 1; // smaller multiples have a smaller prime factor
}
let count = 0;
for (let k = 2; k < n; k++) if (!composite[k]) count++;
return count;
}
const n = 10000000;
let t = Date.now();
let trial = 0;
for (let k = 2; k < n; k++) if (isPrime(k)) trial++;
const a = (Date.now() - t) / 1000;
t = Date.now();
const sieve = countPrimes(n);
const b = (Date.now() - t) / 1000;
console.log(`below ${n}: trial division ${trial} primes in ${a.toFixed(3)}s, the sieve ${sieve} in ${b.toFixed(3)}s`);
// below 10000000: trial division 664579 primes in 10.918s, the sieve 664579 in 0.087sBelow ten million, C++ took 5.17 s by trial division and 0.038 s with the sieve; Java 10.73 s and 0.076 s; JavaScript 10.92 s and 0.087 s. Python stopped at a million, where trial division took 8.20 s and the sieve 0.182 s. Each found the same primes.
Cross Out the Multiples. Run the sieve on 2 to 100 yourself: pick each next prime, decide where its crossing starts, and decide when to stop.
Factorizing
Every n ≥ 2 is a product of primes in exactly one way. To find them, try d from 2 up: divide out every copy of d before moving on, so a d that still divides n is always prime. Stop at √n of what's left; if anything bigger than 1 remains, it's one last prime, larger than the root:
#include <bits/stdc++.h>
using namespace std;
string factorize(long long n) { // n >= 2 as primes with their powers
string out;
for (long long d = 2; d <= n / d; d++) {
int power = 0;
while (n % d == 0) { // divide out every copy of d
n /= d;
power++;
}
if (power) out += (out.empty() ? "" : " * ") + to_string(d) + (power > 1 ? "^" + to_string(power) : "");
}
if (n > 1) out += (out.empty() ? "" : " * ") + to_string(n); // what's left is a prime above the root
return out;
}
int main() {
for (long long n : {360LL, 97LL, 2147483646LL})
printf("%lld = %s\n", n, factorize(n).c_str());
}
// 360 = 2^3 * 3^2 * 5
// 97 = 97
// 2147483646 = 2 * 3^2 * 7 * 11 * 31 * 151 * 331def factorize(n): # n >= 2 as primes with their powers
parts = []
d = 2
while d <= n // d:
power = 0
while n % d == 0: # divide out every copy of d
n //= d
power += 1
if power:
parts.append(f"{d}^{power}" if power > 1 else str(d))
d += 1
if n > 1:
parts.append(str(n)) # what's left is a prime above the root
return " * ".join(parts)
for n in (360, 97, 2147483646):
print(f"{n} = {factorize(n)}")
# 360 = 2^3 * 3^2 * 5
# 97 = 97
# 2147483646 = 2 * 3^2 * 7 * 11 * 31 * 151 * 331public class Main {
static String factorize(long n) { // n >= 2 as primes with their powers
StringBuilder out = new StringBuilder();
for (long d = 2; d <= n / d; d++) {
int power = 0;
while (n % d == 0) { // divide out every copy of d
n /= d;
power++;
}
if (power > 0) out.append(out.length() > 0 ? " * " : "").append(d).append(power > 1 ? "^" + power : "");
}
if (n > 1) out.append(out.length() > 0 ? " * " : "").append(n); // what's left is a prime above the root
return out.toString();
}
public static void main(String[] args) {
for (long n : new long[]{360, 97, 2147483646L})
System.out.println(n + " = " + factorize(n));
}
}
// 360 = 2^3 * 3^2 * 5
// 97 = 97
// 2147483646 = 2 * 3^2 * 7 * 11 * 31 * 151 * 331function factorize(n) { // n >= 2 as primes with their powers
const parts = [];
for (let d = 2; d <= Math.floor(n / d); d++) {
let power = 0;
while (n % d === 0) { // divide out every copy of d
n /= d;
power++;
}
if (power) parts.push(power > 1 ? `${d}^${power}` : `${d}`);
}
if (n > 1) parts.push(`${n}`); // what's left is a prime above the root
return parts.join(" * ");
}
for (const n of [360, 97, 2147483646]) console.log(`${n} = ${factorize(n)}`);
// 360 = 2^3 * 3^2 * 5
// 97 = 97
// 2147483646 = 2 * 3^2 * 7 * 11 * 31 * 151 * 331360 = 2³ × 3² × 5. 2,147,483,646, one less than the prime above, factors as 2 × 3² × 7 × 11 × 31 × 151 × 331. O(√n) per number.
For many numbers up to a limit, sieve once more, keeping each number's smallest prime factor. Then factorizing x is dividing by spf[x] until 1 is left, one division per prime factor, at most log₂ x of them:
#include <bits/stdc++.h>
using namespace std;
const int N = 1000000;
vector<int> spf(N + 1, 0); // spf[x]: the smallest prime that divides x
void buildSpf() {
for (long long p = 2; p <= N; p++) {
if (spf[p] != 0) continue; // already has a smaller prime factor
for (long long m = p; m <= N; m += p)
if (spf[m] == 0) spf[m] = p;
}
}
string factorize(int x) { // x <= N: one division per prime factor
string out;
while (x > 1) {
int p = spf[x];
out += (out.empty() ? "" : " * ") + to_string(p);
x /= p;
}
return out;
}
int main() {
buildSpf();
for (int x : {360, 997, 720720}) printf("%d = %s\n", x, factorize(x).c_str());
}
// 360 = 2 * 2 * 2 * 3 * 3 * 5
// 997 = 997
// 720720 = 2 * 2 * 2 * 2 * 3 * 3 * 5 * 7 * 11 * 13N = 1000000
spf = [0] * (N + 1) # spf[x]: the smallest prime that divides x
for p in range(2, N + 1):
if spf[p] == 0: # nothing smaller divides p: it's prime
for m in range(p, N + 1, p):
if spf[m] == 0:
spf[m] = p
def factorize(x): # x <= N: one division per prime factor
parts = []
while x > 1:
parts.append(spf[x])
x //= spf[x]
return " * ".join(map(str, parts))
for x in (360, 997, 720720):
print(f"{x} = {factorize(x)}")
# 360 = 2 * 2 * 2 * 3 * 3 * 5
# 997 = 997
# 720720 = 2 * 2 * 2 * 2 * 3 * 3 * 5 * 7 * 11 * 13public class Main {
static final int N = 1000000;
static int[] spf = new int[N + 1]; // spf[x]: the smallest prime that divides x
static String factorize(int x) { // x <= N: one division per prime factor
StringBuilder out = new StringBuilder();
while (x > 1) {
int p = spf[x];
out.append(out.length() > 0 ? " * " : "").append(p);
x /= p;
}
return out.toString();
}
public static void main(String[] args) {
for (int p = 2; p <= N; p++) {
if (spf[p] != 0) continue; // already has a smaller prime factor
for (int m = p; m <= N; m += p)
if (spf[m] == 0) spf[m] = p;
}
for (int x : new int[]{360, 997, 720720}) System.out.println(x + " = " + factorize(x));
}
}
// 360 = 2 * 2 * 2 * 3 * 3 * 5
// 997 = 997
// 720720 = 2 * 2 * 2 * 2 * 3 * 3 * 5 * 7 * 11 * 13const N = 1000000;
const spf = new Int32Array(N + 1); // spf[x]: the smallest prime that divides x
for (let p = 2; p <= N; p++) {
if (spf[p] !== 0) continue; // already has a smaller prime factor
for (let m = p; m <= N; m += p) if (spf[m] === 0) spf[m] = p;
}
function factorize(x) { // x <= N: one division per prime factor
const parts = [];
while (x > 1) {
parts.push(spf[x]);
x /= spf[x];
}
return parts.join(" * ");
}
for (const x of [360, 997, 720720]) console.log(`${x} = ${factorize(x)}`);
// 360 = 2 * 2 * 2 * 3 * 3 * 5
// 997 = 997
// 720720 = 2 * 2 * 2 * 2 * 3 * 3 * 5 * 7 * 11 * 13720,720 came apart in 10 divisions, one per prime factor, and no x takes more than log₂ x of them; trial division needs about √x tries when x is prime. The table takes O(n log log n) to build and a cell per number, so it pays off for many numbers up to one limit; Sieve-Based Prime Factorization is this.
Spotting it
The clues: "divisors" or "factors", "prime", "coprime", "greatest common divisor" or "least common multiple", a fraction to reduce, a pattern that divides two strings, or many numbers to test or factor up to one limit (the sieve). Near misses:
- Power of Two asks whether one prime, 2, is the only factor, and a single bit operation answers it (Bit Manipulation).
- Pow(x, n) and Count Good Numbers multiply, not divide: fast power (Modular Arithmetic, Fast Power and Counting).
Mistakes that cost an interview
i * i <= n in 32 bits. The first test said 2,147,483,647 isn't prime, in 10 seconds in Java and at once in C++. Write i <= n / i, or square in 64 bits.
Counting a square's root twice. 36's divisors listed with 6 twice makes ten, not nine. Keep n / d only when it's a different number.
Calling 1 prime. It has one divisor, not two. Start the test with n < 2.
Multiplying first in an LCM. In 32 bits:
#include <bits/stdc++.h>
using namespace std;
int gcd(int a, int b) { return b == 0 ? a : gcd(b, a % b); }
int main() {
int a = 60000, b = 90000;
// a * b wraps first (multiplied as unsigned: an int overflow would be undefined behaviour)
printf("a * b / gcd: %d\n", (int) ((unsigned) a * b) / gcd(a, b));
printf("a / gcd * b: %d\n", a / gcd(a, b) * b);
}
// a * b / gcd: 36834
// a / gcd * b: 180000public class Main {
static int gcd(int a, int b) { return b == 0 ? a : gcd(b, a % b); }
public static void main(String[] args) {
int a = 60000, b = 90000;
System.out.println("a * b / gcd: " + a * b / gcd(a, b)); // a * b wraps first
System.out.println("a / gcd * b: " + a / gcd(a, b) * b);
}
}
// a * b / gcd: 36834
// a / gcd * b: 1800005,400,000,000 wrapped around 2³² before the division, and 36,834 came out. Divide by the gcd first.
Crossing out from p * p for every prime, in an int. Run for every prime below n, the start wraps negative once p passes 46,340. Count Primes at five million:
public class Main {
static int countPrimes(int n) {
boolean[] composite = new boolean[Math.max(n, 2)];
int count = 0;
for (int p = 2; p < n; p++) {
if (composite[p]) continue;
count++;
for (int m = p * p; m < n; m += p) composite[m] = true; // p * p wraps once p passes 46,340
}
return count;
}
public static void main(String[] args) {
System.out.println(countPrimes(5000000));
}
}It threw ArrayIndexOutOfBoundsException at index −2,146,737,495, in an array of length 5,000,000. Stop the outer loop once p × p reaches n, which also skips work, or compute the start in 64 bits.
Check yourself
Quiz
1 / 6Multiple choice
int count = 0;
for (int d = 1; d <= n / d; d++)
count++;Multiple choice
Multiple choice
Multiple choice
Multiple choice
Multiple choice
Practice
In this section, easy to hard:
- Finding All Divisors and All divisors of a Number: pairs up to √n, the root once.
- Greatest Common Divisor of 2 numbers: Euclid.
- Greatest Common Divisor of Strings: s + t against t + s, then gcd of the lengths.
- Primality Check: trial division with a bound that can't overflow.
- Count Primes: the sieve, from p × p.
- Prime Factorization: divide out each d up to √n.
- Sieve-Based Prime Factorization: the smallest prime factor table.
Next: Modular Arithmetic, Fast Power and Counting, where numbers get too big to hold and the answer is wanted mod 10⁹ + 7.