17 — Bits and mathematics

A large share of algorithm problems turn out to be arithmetic in disguise. Testing whether a number is a power of two, finding the one element without a partner, enumerating every subset of a small set, computing $\binom{n}{r}$ modulo a prime, picking a uniform random line from a log file of unknown length, deciding whether two line segments cross: each has a short, exact answer that rests on a mathematical fact rather than on a data structure. This chapter collects those facts and the code that uses them.

The chapter has five parts. Bit manipulation comes first, starting with how integers are actually represented, because Python differs from C and Java in ways that break code ported naively. Then number theory (primes, GCD, modular arithmetic, exponentiation), combinatorics, randomized sampling, and finally the small amount of computational geometry that appears in interviews and in practice: orientation, segment intersection and convex hulls.

Every implementation below was checked against a brute-force or standard-library reference on thousands of random inputs. Timings were measured with Python 3.14 on one Mac with timeit, taking the best of several repeats. Exact numbers will differ elsewhere; the ratios are what matter.

How integers are represented

Two’s complement

Hardware stores a fixed-width signed integer in two’s complement. With $w$ bits, the pattern $b_{w-1} \dots b_1 b_0$ means

$-b_{w-1} 2^{w-1} + \sum_{i=0}^{w-2} b_i 2^i$

The top bit carries a negative weight. For 8 bits:

pattern     unsigned   signed (two's complement)
00000000        0          0
00000001        1          1
01111111      127        127      largest positive
10000000      128       -128      most negative
11111010      250         -6
11111111      255         -1

Three properties follow and are used constantly:

• Negation is “invert and add one”: $-x = \lnot x + 1$. Equivalently $\lnot x = -x - 1$.
• Addition is the same circuit for signed and unsigned values. The hardware adds bit patterns modulo $2^w$; only the interpretation differs.
• The range is asymmetric: $-2^{w-1}$ to $2^{w-1} - 1$. Negating the most negative value overflows back to itself in C or Java, which is the root of the INT_MIN / -1 edge case in “Divide Two Integers” below.

Python’s integers have no width

A Python int is arbitrary precision. It never overflows; 2**100 + 1 is exact. For bitwise operations, Python behaves as if every integer were written in two’s complement with an infinite number of sign bits extending to the left: a non-negative number has infinitely many leading 0s, a negative number infinitely many leading 1s. That model explains every surprise:

x = -6
print(bin(x), x & 0xFF, bin(x & 0xFF))
print(~5, -5 >> 1, -1 >> 100, 5 >> 1)
print((-6).bit_length(), (6).bit_length())
print(format(-6 & 0xFFFFFFFF, '032b'))

-0b110 250 0b11111010
-6 -3 -1 2
3 3
11111111111111111111111111111010

• bin(-6) prints a minus sign and the magnitude, not a bit pattern. To see the pattern a 32-bit machine would hold, mask it: x & 0xFFFFFFFF. Masking with a positive number cuts the infinite row of 1s down to $w$ bits and yields a non-negative result.
• ~5 is -6. NOT flips infinitely many leading 0s into 1s, so ~x == -x - 1 always.
• >> on a negative number is an arithmetic shift: it floors. -5 >> 1 is -3 (floor of $-2.5$), and shifting -1 right any distance stays -1. Java’s >>> (logical shift) does not exist because there is no fixed width to shift zeros into.
• << never overflows. 1 << 64 is a 65-bit number, not 0.
• bit_length() ignores the sign.

Going the other way, from a 32-bit unsigned pattern back to a signed value, needs an explicit sign fix:

def to_signed32(u):
    u &= 0xFFFFFFFF
    return u - (1 << 32) if u & 0x80000000 else u

print(to_signed32(0xFFFFFFFA), to_signed32(2**31), to_signed32(2**31 - 1))

-6 -2147483648 2147483647

This matters in any problem specified for 32-bit integers. The usual pattern is: mask every intermediate result with & 0xFFFFFFFF, then convert at the end. Without the mask, a loop that relies on overflow to terminate (for example, adding two integers using only XOR and AND for the carry) runs forever on negative inputs in Python, because the carry never falls off the top.

Division is the other trap. Python’s // floors and % takes the sign of the divisor; C, C++ and Java truncate toward zero:

print(-7 // 2, -7 % 2, int(-7 / 2))

-4 1 -3

int(a / b) truncates but goes through a float, which loses precision past $2^{53}$. For exact truncating division of large integers, compute on absolute values and reapply the sign.

Bitwise operators

print(12 & 10, 12 | 10, 12 ^ 10, ~12, 12 << 2, 12 >> 2)

8 14 6 -13 48 3

With $12 = 1100_2$ and $10 = 1010_2$: AND keeps the shared $1000_2$, OR gives $1110_2$, XOR gives $0110_2$.

The single-bit idioms, with i counting from 0 at the least significant bit:

Precedence differs between languages. In Python, arithmetic binds tighter than shifts, shifts tighter than &, then ^, then |, and all bitwise operators bind tighter than comparisons, so x & 1 == 0 means (x & 1) == 0. In C, C++ and Java, == binds tighter than &, so the same text means x & (1 == 0), which is always 0. Code translated from Python to C breaks silently on this. Both languages agree that 1 << n - 1 is 1 << (n - 1), which surprises people in the other direction. Parenthesizing bit expressions is cheap insurance.

XOR has the algebra that makes most tricks work. It is associative and commutative, $x \oplus x = 0$, and $x \oplus 0 = x$. So XOR of a multiset cancels every value that appears an even number of times, regardless of order.

The two key identities: x & (x - 1) and x & -x

Subtracting 1 from $x$ flips the lowest set bit to 0 and every 0 below it to 1:

x       = 1011000
x - 1   = 1010111
x & x-1 = 1010000     lowest 1 cleared

And because $-x = \lnot x + 1$, negation flips every bit above the lowest set bit and leaves that bit (and the zeros below it) unchanged:

x       = ...0001011000
-x      = ...1110101000
x & -x  =          1000     only the lowest 1 survives

print(bin(0b1011000 & -0b1011000), 88 & -88)

0b1000 8

x & -x is called the lowbit. It drives the Fenwick tree (binary indexed tree), where the index arithmetic i += i & -i and i -= i & -i walks the implicit tree; see chapter 11. It also gives the index of the lowest set bit as (x & -x).bit_length() - 1.

Power of two

A positive power of two has exactly one set bit, so clearing its lowest set bit leaves zero:

def is_power_of_two(n):
    return n > 0 and n & (n - 1) == 0

print([n for n in range(-4, 70) if is_power_of_two(n)])

[1, 2, 4, 8, 16, 32, 64]

The n > 0 guard is required: 0 & -1 is 0, so without it zero passes. Negative numbers fail correctly in Python, but in C, INT_MIN & (INT_MIN - 1) is undefined behavior, so the guard matters there too.

Related checks: a power of four is a power of two whose single bit sits at an even position, n & 0x55555555 != 0 for 32-bit inputs, or equivalently n % 3 == 1. The next power of two at or above n (for n >= 1) is 1 << (n - 1).bit_length().

Counting set bits

The number of 1 bits is called the population count or Hamming weight. Three ways to compute it in Python:

def popcount_kernighan(n):
    count = 0
    while n:
        n &= n - 1
        count += 1
    return count

print(popcount_kernighan(0b1011000), (0b1011000).bit_count())

3 3

Kernighan’s loop clears the lowest set bit each iteration, so it runs once per set bit rather than once per bit position: $O(\text{popcount})$ iterations instead of $O(w)$. It is commonly called Brian Kernighan’s algorithm because it appears in The C Programming Language; Peter Wegner published it in 1960. It only terminates for non-negative n in Python, since a negative number has infinitely many 1 bits. Mask first (n & 0xFFFFFFFF) if the problem means a 32-bit pattern.

int.bit_count() (Python 3.10 and later) counts the 1 bits of abs(n) in C, using the hardware popcount instruction where available. Measured over 10,000 random 64-bit integers:

The interpreted loop runs about 32 iterations per random 64-bit number and is over 100 times slower than bit_count. In Python, the loop is worth knowing for its idea (the x & (x - 1) step reappears in many problems) but not for production counting.

A classic use is Counting Bits: produce the popcount of every integer from 0 to $n$ in $O(n)$, using the recurrence bits[i] = bits[i & (i - 1)] + 1 (or bits[i >> 1] + (i & 1)). Each entry depends on a smaller one already computed. Hamming Distance between two integers is (a ^ b).bit_count(): XOR marks the positions where they differ.

XOR tricks

Single Number

Every element appears twice except one. XOR everything together; pairs cancel and the singleton remains. $O(n)$ time, $O(1)$ space, which beats the obvious hash-set approach on memory:

def single_number(nums):
    acc = 0
    for x in nums:
        acc ^= x
    return acc

The trick depends on “appears twice” exactly (or any even count). If every other element appears three times, XOR no longer cancels. The general fix counts, for each bit position, how many numbers have that bit set; the count modulo 3 gives the singleton’s bit. In Python that needs care with negatives: either sum each of 32 bit positions and convert the result with to_signed32, or use the two-variable state machine ones = (ones ^ x) & ~twos; twos = (twos ^ x) & ~ones, which works directly on Python’s infinite two’s complement.

Missing Number

The array holds $n$ distinct values from $0..n$ with one missing. XOR every index $0..n$ with every value; each present value cancels its own index, leaving the missing one:

def missing_number(nums):
    acc = len(nums)
    for i, x in enumerate(nums):
        acc ^= i ^ x
    return acc

acc starts at len(nums) because the indices only go up to $n-1$ and the value $n$ also needs a partner. The arithmetic alternative, $n(n+1)/2 - \sum$, is equally correct in Python; in fixed-width languages the sum can overflow while the XOR cannot.

Two non-repeating elements

Every element appears twice except two, $a$ and $b$. XOR of everything is $a \oplus b$, which is non-zero because $a \ne b$. Any set bit of $a \oplus b$ is a position where $a$ and $b$ differ. Take the lowest one with x & -x, split the array by that bit, and XOR each half: the pairs still cancel within their half, and $a$ and $b$ land in different halves.

def two_singles(nums):
    xor_all = 0
    for x in nums:
        xor_all ^= x
    low = xor_all & -xor_all          # a bit where the two differ
    a = 0
    for x in nums:
        if x & low:
            a ^= x
    return sorted((a, xor_all ^ a))

print(single_number([4, 1, 2, 1, 2]), missing_number([3, 0, 1]), two_singles([1, 2, 1, 3, 2, 5]))

4 2 [3, 5]

nums: 1 2 1 3 2 5          xor_all = 3 ^ 5 = 011 ^ 101 = 110,  low = 010
bit 1 set:    2 3 2        -> 3
bit 1 clear:  1 1 5        -> 5

This works unchanged for negative numbers in Python: x & -x of a negative xor_all still isolates its lowest set bit.

Swapping with XOR, and why not to

a ^= b      a = a0 ^ b0
b ^= a      b = b0 ^ a0 ^ b0 = a0
a ^= b      a = a0 ^ b0 ^ a0 = b0

It swaps without a temporary. It is a curiosity rather than a technique, for three reasons. First, it fails when both names refer to the same storage: swapping a[i] with a[j] when i == j zeroes the element.

a = [5, 9]; i = j = 0
a[i] ^= a[j]; a[j] ^= a[i]; a[i] ^= a[j]
print(a)

[0, 9]

Second, it is slower. Each step depends on the previous one, so a CPU cannot overlap them, whereas a register swap is often free after renaming; compilers recognize the temporary-variable swap and emit the best code. In Python, a, b = b, a measured 6 ns against 30 ns for the three XOR statements. Third, it only works for integers. Use tuple assignment.

Reversing bits

“Reverse Bits” asks for the 32-bit unsigned integer whose bits are the input’s in reverse order. The direct loop shifts bits out of n from the bottom and into result from the bottom, so the first bit out ends up on top after 32 steps:

def reverse_bits(n, width=32):
    result = 0
    for _ in range(width):
        result = (result << 1) | (n & 1)
        n >>= 1
    return result

print(reverse_bits(0b1011), format(reverse_bits(0b1011), '032b'))

3489660928 11010000000000000000000000000000

The width argument matters: the input’s leading zeros become the output’s trailing zeros, so the loop must run 32 times, not n.bit_length() times.

A branch-free alternative swaps adjacent bits, then adjacent pairs, then nibbles, bytes and half-words, five steps of masks and shifts:

def reverse_bits_swap(n):
    n = ((n >> 1) & 0x55555555) | ((n & 0x55555555) << 1)
    n = ((n >> 2) & 0x33333333) | ((n & 0x33333333) << 2)
    n = ((n >> 4) & 0x0F0F0F0F) | ((n & 0x0F0F0F0F) << 4)
    n = ((n >> 8) & 0x00FF00FF) | ((n & 0x00FF00FF) << 8)
    n = ((n >> 16) & 0x0000FFFF) | ((n & 0x0000FFFF) << 16)
    return n

print(reverse_bits(43261596), reverse_bits_swap(43261596))

964176192 964176192

0x55555555 is 0101...01, 0x33333333 is 0011...0011, and so on; each line exchanges two halves of every block of size $2^k$. Because Python’s left shift does not discard high bits, the masks must be applied before shifting left, as written, or the result grows past 32 bits. In Python the simplest correct version is int(format(n, '032b')[::-1], 2). If the function is called many times, a 256-entry table of reversed bytes reduces the work to four lookups.

Bitmasks as sets

An integer with $n$ bits is a subset of $\{0, 1, \dots, n-1\}$: bit $i$ set means element $i$ is in. Set operations become single machine instructions:

Enumerating all subsets

The integers $0$ to $2^n - 1$ are exactly the $2^n$ subsets:

def subsets(items):
    n = len(items)
    return [[items[i] for i in range(n) if mask >> i & 1] for mask in range(1 << n)]

print(subsets(["a", "b", "c"]))

[[], ['a'], ['b'], ['a', 'b'], ['c'], ['a', 'c'], ['b', 'c'], ['a', 'b', 'c']]

$O(n 2^n)$ time to materialize, which is unavoidable since the output has that size. Practical up to about $n = 20$ (a million subsets). This is the iterative alternative to the backtracking subset generator in chapter 15.

Enumerating the submasks of a mask

To visit every subset of a given mask m (not of the full set), use sub = (sub - 1) & m. Subtracting 1 clears the lowest set bit of sub and sets every bit below it; ANDing with m keeps only bits that belong to m. The result is the next smaller submask:

def submasks(mask):
    sub = mask
    while sub:
        yield sub
        sub = (sub - 1) & mask
    yield 0

print([format(s, '04b') for s in submasks(0b1011)])

['1011', '1010', '1001', '1000', '0011', '0010', '0001', '0000']

Iterating over all masks and all their submasks costs $3^n$ in total, not $4^n$: each (mask, submask) pair corresponds to choosing, for every element independently, one of three states (outside the mask, in the mask but not the submask, or in both). This bound is what makes “DP over subsets of subsets” feasible for $n \approx 15$.

DP over subsets

Bitmasks index DP tables whose state is “which items have been used”. The Held-Karp Traveling Salesman algorithm stores dp[mask][last], the cheapest path visiting exactly the cities in mask and ending at last, in $O(2^n n^2)$ time (chapter 13). Assignment problems (“give each of $n$ workers one of $n$ jobs”) use dp[mask] with the number of workers already assigned equal to mask.bit_count(). The DP chapter covers these in depth (chapter 14); the bit-level operations needed are the ones in the table above.

Bitmasks are limited by $2^n$ growth, not by integer width: Python integers can hold a 1,000-bit set fine, and &/| on them are fast, which makes Python ints a reasonable bitset for things like reachability rows. What cannot scale is enumerating all subsets.

Divide Two Integers

The problem: divide two 32-bit signed integers without *, / or %, truncating toward zero, and return $2^{31} - 1$ on overflow. Repeated subtraction is $O(\text{quotient})$, up to two billion steps. Shifts make it $O(\log)$: this is binary long division. For each bit position from high to low, check whether the divisor shifted to that position still fits into what remains of the dividend; if so, subtract it and set that quotient bit.

INT_MAX, INT_MIN = 2**31 - 1, -2**31

def divide(dividend, divisor):
    if dividend == INT_MIN and divisor == -1:
        return INT_MAX                   # the only overflowing case
    negative = (dividend < 0) != (divisor < 0)
    a, b = abs(dividend), abs(divisor)
    quotient = 0
    for shift in range(31, -1, -1):
        if (a >> shift) >= b:            # b << shift fits in a
            a -= b << shift
            quotient |= 1 << shift
    return -quotient if negative else quotient

print(divide(10, 3), divide(7, -3), divide(-2**31, -1), divide(-2**31, 1))

3 -2 2147483647 -2147483648

Testing (a >> shift) >= b instead of a >= (b << shift) avoids overflow in languages with fixed-width integers. Working on absolute values gives truncation toward zero, which is what the problem specifies and what Python’s // does not do. In C or Java, abs(INT_MIN) itself overflows, so implementations there either use 64-bit intermediates or work entirely in negative numbers. Tested against truncating division on 100,000 random pairs, including the boundary values.

Prime numbers

Sieve of Eratosthenes

To find every prime up to $n$, start with all numbers marked prime and, for each prime $p$ found in order, cross out its multiples. Two observations make it fast:

• Crossing out can start at $p^2$. Any smaller multiple $kp$ with $k < p$ has a prime factor smaller than $p$ and was already crossed out.
• The outer loop can stop at $\sqrt{n}$. A composite $\le n$ has a prime factor $\le \sqrt{n}$.

For $n = 30$:

start:  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
p=2:    cross 4 6 8 10 12 14 16 18 20 22 24 26 28 30
p=3:    cross 9 15 21 27                (6, 12, ... already gone)
p=5:    cross 25                         (start at 25 = 5*5)
p=7:    7*7 = 49 > 30, stop
primes: 2 3 5 7 11 13 17 19 23 29

import math

def sieve(n):
    """All primes <= n."""
    if n < 2:
        return []
    is_prime = bytearray([1]) * (n + 1)
    is_prime[0] = is_prime[1] = 0
    for p in range(2, math.isqrt(n) + 1):
        if is_prime[p]:
            is_prime[p*p::p] = bytes(len(range(p*p, n + 1, p)))
    return [i for i, flag in enumerate(is_prime) if flag]

print(sieve(50))

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

Complexity. Crossing out multiples of $p$ costs about $n/p$. Summed over primes $p \le \sqrt{n}$, that is $n \sum 1/p$, and the sum of reciprocals of primes up to $x$ grows as $\ln \ln x$ (Mertens), giving $O(n \log \log n)$ time. $\log \log n$ is under 3 for any $n$ that fits in memory, so the sieve is linear in practice. Space is $O(n)$: one byte per number in the bytearray, or one bit with a packed bitset. Storing only odd numbers halves both.

Implementation note. The slice assignment is_prime[p*p::p] = bytes(...) crosses out all multiples of p in one C-level operation. Compared with an inner Python for loop:

Both grow roughly linearly. A large share of the slice version’s remaining time is the final list comprehension that collects the primes; the crossing out itself is cheap.

Linear sieve. A variant (often attributed to Gries and Misra, 1978) crosses out each composite exactly once, by its smallest prime factor, in $O(n)$. It also yields the smallest prime factor of every number, which allows factoring any $m \le n$ in $O(\log m)$. In Python the standard sieve with slice assignment is faster in wall-clock time; the linear sieve is worth using when the smallest-prime-factor table is the goal. (A plain sieve can also record smallest prime factors, by writing $p$ into each unmarked multiple.)

Segmented sieve

For primes in a window $[L, R]$ with $R$ large (say $10^{12}$) but $R - L$ modest, a full sieve to $R$ is impossible. Every composite in the window has a prime factor $\le \sqrt{R}$, so: sieve the base primes up to $\sqrt{R}$, then for each base prime cross out its multiples inside the window only, starting from the first multiple $\ge \max(p^2, L)$.

def segmented_sieve(lo, hi):
    """All primes in [lo, hi]."""
    base = sieve(math.isqrt(hi))
    lo = max(lo, 2)
    if lo > hi:
        return []
    seg = bytearray([1]) * (hi - lo + 1)
    for p in base:
        start = max(p * p, (lo + p - 1) // p * p)
        seg[start - lo::p] = bytes(len(range(start, hi + 1, p)))
    return [lo + i for i, flag in enumerate(seg) if flag]

print(segmented_sieve(100, 160))

[101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157]

(lo + p - 1) // p * p is the smallest multiple of p at or above lo (ceiling division). The max(p * p, ...) term stops a base prime that lies inside the window from crossing out itself.

Memory is $O(\sqrt{R} + (R - L))$. Finding all primes in $[10^{12}, 10^{12} + 10^6]$ took 52.2 ms and returned 36,249 primes, close to the $10^6 / \ln 10^{12} \approx 36{,}191$ that the prime number theorem predicts. The same idea, processing $[0, n]$ in cache-sized blocks, is how fast sieves for very large $n$ are built: it keeps the working set in L1 or L2 cache.

Primality by trial division

For a single number, test divisors up to $\sqrt{n}$:

def is_prime_trial(n):
    if n < 2:
        return False
    if n % 2 == 0:
        return n == 2
    for d in range(3, math.isqrt(n) + 1, 2):
        if n % d == 0:
            return False
    return True

$O(\sqrt{n})$ divisions. Use math.isqrt, not int(n ** 0.5): the float square root is inexact for large n and can be off by one, which either skips the last divisor or does one extra. Trial division is fine up to about $10^{12}$ (500,000 odd divisors, 9.3 ms measured for the prime 999,999,999,989). For $10^{18}$ it took 9.64 s, and it grows by $\sqrt{10}$ for every extra digit.

Miller-Rabin

Fermat’s little theorem says that for prime $p$ and $a$ not divisible by $p$, $a^{p-1} \equiv 1 \pmod p$. Checking that for a few values of $a$ is the Fermat test, but some composites (Carmichael numbers, the smallest being $561 = 3 \cdot 11 \cdot 17$) satisfy it for every $a$ coprime to them. Miller-Rabin strengthens the test with a second fact: modulo a prime, the only square roots of 1 are $1$ and $-1$.

Write $n - 1 = d \cdot 2^s$ with $d$ odd. For a base $a$, compute $x = a^d \bmod n$, then square it $s - 1$ times. If $n$ is prime, the sequence $a^d, a^{2d}, \dots, a^{2^s d}$ ends at 1 (Fermat), and the first 1 must be preceded by $-1$ unless the sequence starts at 1. So $n$ is declared probably prime for base $a$ if $x = 1$ or some term equals $n - 1$. Otherwise $a$ is a witness that $n$ is composite, with certainty.

n = 221 = 13 * 17,  n - 1 = 220 = 55 * 2^2,  d = 55, s = 2

a = 174:  174^55 mod 221 = 47,  47^2 mod 221 = 220 = n-1   -> passes (a "strong liar")
a = 137:  137^55 mod 221 = 188, 188^2 mod 221 = 205        -> never hit 1 or n-1: composite

For a composite $n$, at most a quarter of bases are strong liars (Rabin, 1980), so $k$ random bases give error probability at most $4^{-k}$. For bounded $n$ there are fixed base sets that are proven to make the test deterministic. The first twelve primes, 2 through 37, are correct for all $n < 3.18 \times 10^{23}$, which covers every 64-bit integer:

MR_BASES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37)

def is_prime_mr(n):
    if n < 2:
        return False
    for p in MR_BASES:
        if n % p == 0:
            return n == p
    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2
        s += 1
    for a in MR_BASES:
        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(s - 1):
            x = x * x % n
            if x == n - 1:
                break
        else:
            return False             # a is a witness: definitely composite
    return True

The trial division by the base primes first handles small $n$ (where a base might equal or exceed $n$) and removes most composites cheaply. The built-in three-argument pow(a, d, n) is modular exponentiation in C.

The bound is exact, not conservative. The smallest composite that passes all twelve bases is

print(399165290221 * 798330580441)

318665857834031151167461

and the function above does return True for it. Beyond that bound, use random bases or a larger fixed set. A shorter seven-base set, {2, 325, 9375, 28178, 450775, 9780504, 1795265022}, found by Jim Sinclair, is also deterministic below $2^{64}$; with those bases, reduce a % n and skip a base that becomes 0.

Checked against the sieve for every $n \le 10^6$, against trial division for 3,000 random $n$ up to $10^{12}$, and on known strong pseudoprimes (2047, 3215031751, 3825123056546413051) and Carmichael numbers.

Miller-Rabin costs $O(k \log^3 n)$ with schoolbook multiplication ($k$ bases, $O(\log n)$ squarings, each $O(\log^2 n)$). It answers “is this prime”; it does not produce factors. For factoring 64-bit numbers, Pollard’s rho with Miller-Rabin as the primality check is the usual combination.

GCD, LCM and the Euclidean algorithm

The greatest common divisor rests on one identity: $\gcd(a, b) = \gcd(b, a \bmod b)$, because any common divisor of $a$ and $b$ also divides $a - qb$, and vice versa. Applying it until the second argument is 0 gives Euclid’s algorithm (Elements, Book VII, around 300 BC):

def gcd(a, b):
    while b:
        a, b = b, a % b
    return abs(a)

def lcm(a, b):
    return abs(a // gcd(a, b) * b) if a and b else 0

gcd(252, 105):  252 = 2*105 + 42
                105 = 2*42  + 21
                 42 = 2*21  + 0     -> 21

Complexity. $O(\log \min(a, b))$ iterations. Every two steps at least halve the larger argument. The worst case is consecutive Fibonacci numbers, where each quotient is 1 (Lamé, 1844), so the number of steps is at most five times the number of decimal digits of the smaller argument.

LCM. $\operatorname{lcm}(a, b) = |ab| / \gcd(a, b)$. Divide before multiplying (a // g * b) so that fixed-width languages do not overflow on the product.

In Python, use math.gcd and math.lcm, which accept any number of arguments (since 3.9) and run in C. The versions above match them on random inputs including negatives and zero.

Extended Euclid

The extended algorithm also finds integers $x, y$ with $ax + by = \gcd(a, b)$ (Bézout’s identity). It keeps, alongside each remainder, the coefficients that express that remainder as a combination of the original $a$ and $b$. Each step applies the same update $r_{new} = r_{old} - q r$ to the coefficients:

def extended_gcd(a, b):
    """Return (g, x, y) with a*x + b*y == g == gcd(a, b)."""
    old_r, r = a, b
    old_x, x = 1, 0
    old_y, y = 0, 1
    while r:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_x, x = x, old_x - q * x
        old_y, y = y, old_y - q * y
    return old_r, old_x, old_y

  q    r    x    y          invariant: r == 240*x + 46*y
     240    1    0
      46    0    1
  5   10    1   -5
  4    6   -4   21
  1    4    5  -26
  1    2   -9   47          240*(-9) + 46*47 = 2
  2    0   23 -120

Result: extended_gcd(240, 46) returns (2, -9, 47). The iterative form avoids Python’s recursion limit and is no harder to read once the invariant is stated.

Uses: modular inverses (next section), solving linear Diophantine equations $ax + by = c$ (solvable exactly when $\gcd(a, b)$ divides $c$), and the Chinese Remainder Theorem for combining congruences with coprime moduli.

Modular arithmetic

Problems whose answers are astronomically large usually ask for the result modulo a prime, almost always $10^9 + 7$ or $998244353$. Reduction commutes with addition, subtraction and multiplication:

$(a + b) \bmod m = ((a \bmod m) + (b \bmod m)) \bmod m$

$(a \cdot b) \bmod m = ((a \bmod m) \cdot (b \bmod m)) \bmod m$

so every intermediate value can be reduced and stays small. $10^9 + 7$ is chosen because it is prime (so every non-zero value has an inverse) and because the product of two reduced values, under $10^{18}$, fits in a signed 64-bit integer. $998244353 = 119 \cdot 2^{23} + 1$ is chosen because it supports the number-theoretic transform for fast polynomial multiplication.

In Python, integers never overflow, so reducing is about speed: arithmetic on numbers with millions of digits is slow. Reduce after every multiplication inside loops. Python’s % always returns a value in $[0, m)$ for positive $m$, even when the left side is negative; in C and Java, (a - b) % m can be negative and needs + m.

Division does not commute with reduction. $(a / b) \bmod m$ is not $(a \bmod m) / (b \bmod m)$. Instead, multiply by the modular inverse of $b$: the value $b^{-1}$ with $b \cdot b^{-1} \equiv 1 \pmod m$. It exists exactly when $\gcd(b, m) = 1$. Two ways to compute it:

• Extended Euclid gives $bx + my = 1$, so $x$ is the inverse. Works for any modulus coprime to $b$.
• Fermat’s little theorem gives $b^{m-2} \equiv b^{-1} \pmod m$ when $m$ is prime. One call to pow(b, m - 2, m).

Python 3.8 and later computes inverses directly with pow(b, -1, m), raising ValueError if none exists.

def mod_inverse(a, m):
    g, x, _ = extended_gcd(a % m, m)
    if g != 1:
        raise ValueError(f"{a} has no inverse mod {m}")
    return x % m

print(mod_inverse(3, 11), pow(3, -1, 11), (3 * 4) % 11)
try:
    mod_inverse(6, 9)
except ValueError as e:
    print(e)

4 4 1
6 has no inverse mod 9

The Fermat approach gives a wrong answer silently when the modulus is not prime, and when $b$ is a multiple of $m$ it returns 0. Extended Euclid reports the failure.

Fast exponentiation

Computing $a^e$ by repeated multiplication takes $e - 1$ steps, hopeless for $e = 10^{18}$. Binary exponentiation uses the binary digits of $e$: $a^{13} = a^{8} \cdot a^{4} \cdot a^{1}$ since $13 = 1101_2$. Square the base each step, and multiply it into the result when the current bit is set:

def power(base, exp, mod):
    result = 1
    base %= mod
    while exp:
        if exp & 1:
            result = result * base % mod
        base = base * base % mod
        exp >>= 1
    return result % mod

3^13 mod 1000
exp=13 (1101)  bit=1  base=  3  result=1     -> result = 3
exp= 6 (0110)  bit=0  base=  9  result=3
exp= 3 (0011)  bit=1  base= 81  result=3     -> result = 243
exp= 1 (0001)  bit=1  base=561  result=243   -> result = 323
answer: 323     (561 = 6561 mod 1000)

$O(\log e)$ multiplications: at most $\lfloor \log_2 e \rfloor$ squarings and as many multiplications. The final % mod handles mod == 1, where the answer must be 0 even for exp == 0. Python’s built-in pow(base, exp, mod) does the same thing in C and should be used in real code; the hand-written version matters because the same loop works for anything with an associative multiplication: matrices, polynomials, permutations, or function composition.

The idea also has a right-to-left counterpart (process the high bits first, squaring the result and multiplying by the base when a bit is set), which is what most libraries do.

Matrix exponentiation

A linear recurrence of order $k$ can be written as a $k \times k$ matrix acting on a vector of the last $k$ terms. For Fibonacci:

$\begin{pmatrix} F_{n+1} \\ F_n \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} F_n \\ F_{n-1} \end{pmatrix}, \qquad \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n = \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix}$

Raising the matrix to the $n$-th power with binary exponentiation gives $F_n$ in $O(k^3 \log n)$ arithmetic operations:

def mat_mult(A, B, mod):
    return [[sum(A[i][k] * B[k][j] for k in range(len(B))) % mod
             for j in range(len(B[0]))] for i in range(len(A))]

def mat_pow(M, exp, mod):
    n = len(M)
    result = [[int(i == j) for j in range(n)] for i in range(n)]   # identity
    while exp:
        if exp & 1:
            result = mat_mult(result, M, mod)
        M = mat_mult(M, M, mod)
        exp >>= 1
    return result

def fib_matrix(n, mod):
    if n == 0:
        return 0
    return mat_pow([[1, 1], [1, 0]], n, mod)[0][1]

def fib_iter(n, mod):
    a, b = 0, 1
    for _ in range(n):
        a, b = b, (a + b) % mod
    return a

print([fib_matrix(n, MOD) for n in range(12)])
print(fib_matrix(10**18, MOD))

[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89]
209783453

(MOD is $10^9 + 7$ throughout.) Measured against the iterative loop:

The iterative column grows tenfold per row; the matrix column grows by a few microseconds per row, one extra squaring per bit of $n$. With these list-based matrices the crossover is around $n = 1000$. Fast doubling uses the identities $F_{2k} = F_k (2F_{k+1} - F_k)$ and $F_{2k+1} = F_k^2 + F_{k+1}^2$, which are the matrix method with the redundant entries removed, and is about ten times faster again.

The technique generalizes: any recurrence $a_n = c_1 a_{n-1} + \dots + c_k a_{n-k}$ uses a companion matrix with the coefficients in the first row and an identity shifted below. Counting paths of length exactly $L$ in a graph is the $L$-th power of the adjacency matrix. For large $k$, the $k^3$ factor dominates and polynomial methods (Kitamasa) do better.

Without a modulus, $F_n$ has about $0.694n$ bits, so for large $n$ the cost is dominated by big-integer multiplication rather than the count of operations.

Combinatorics

Factorials and binomial coefficients

$\binom{n}{r} = \frac{n!}{r!(n-r)!}$ counts the $r$-element subsets of an $n$-element set. In Python, math.comb(n, r) computes it exactly, and math.factorial and math.perm cover the other cases. Exact values grow fast ($\binom{10^6}{5 \times 10^5}$ has about 300,000 digits), so contest and interview problems ask for $\binom{n}{r} \bmod p$.

The standard technique precomputes factorials and inverse factorials modulo a prime $p$ once, after which each query is three multiplications:

MOD = 10**9 + 7

class Binomial:
    def __init__(self, n_max, mod=MOD):
        self.mod = mod
        fact = [1] * (n_max + 1)
        for i in range(1, n_max + 1):
            fact[i] = fact[i - 1] * i % mod
        inv = [1] * (n_max + 1)
        inv[n_max] = pow(fact[n_max], mod - 2, mod)   # Fermat inverse
        for i in range(n_max, 0, -1):
            inv[i - 1] = inv[i] * i % mod             # 1/(i-1)! = i / i!
        self.fact, self.inv_fact = fact, inv

    def C(self, n, r):
        if r < 0 or r > n:
            return 0
        return self.fact[n] * self.inv_fact[r] % self.mod * self.inv_fact[n - r] % self.mod

B = Binomial(10**6)
print(B.C(10, 3), B.C(10**6, 500_000), math.comb(10**6, 500_000) % MOD)

120 996692777 996692777

Only one modular inverse is computed, for $n_{max}!$; the rest come from $\frac{1}{(i-1)!} = i \cdot \frac{1}{i!}$, walking down. Precomputation is $O(n_{max} + \log p)$ time and $O(n_{max})$ space; each query is $O(1)$.

Preconditions: $p$ must be prime and $n_{max} < p$, because otherwise $n!$ contains the factor $p$, is 0 modulo $p$, and has no inverse. For $n \ge p$ (small primes), use Lucas’s theorem: write $n$ and $r$ in base $p$, and $\binom{n}{r} \equiv \prod \binom{n_i}{r_i} \pmod p$. For a composite modulus, factor it and combine with the Chinese Remainder Theorem.

Pascal’s triangle

The recurrence $\binom{n}{r} = \binom{n-1}{r-1} + \binom{n-1}{r}$ (choose whether element $n$ is in the subset) builds every binomial coefficient by addition alone:

def pascal(rows):
    tri = [[1]]
    for _ in range(rows - 1):
        prev = tri[-1]
        tri.append([1] + [prev[i] + prev[i + 1] for i in range(len(prev) - 1)] + [1])
    return tri

for row in pascal(6):
    print(row)

[1]
[1, 1]
[1, 2, 1]
[1, 3, 3, 1]
[1, 4, 6, 4, 1]
[1, 5, 10, 10, 5, 1]

$O(n^2)$ time and space for $n$ rows, or $O(n)$ space for a single row by updating right to left. It needs no division, so it works modulo any number, prime or not. It is the right choice when $n$ is small (a few thousand) and the modulus is composite, or when many coefficients with small $n$ are needed. Row $n$ sums to $2^n$.

Permutations and combinations

“Selections with repetition” is stars and bars: placing $r$ identical items into $n$ bins is choosing positions for $n - 1$ dividers among $n + r - 1$ slots. Generating the actual permutations and combinations (not counting them) is covered in chapter 15; itertools.permutations and itertools.combinations do it lazily.

Catalan numbers

$C_n = \frac{1}{n+1}\binom{2n}{n} = \binom{2n}{n} - \binom{2n}{n+1}$

with recurrence $C_0 = 1$, $C_{n+1} = \sum_{i=0}^{n} C_i C_{n-i}$.

def catalan(n):
    return math.comb(2 * n, n) // (n + 1)

print([catalan(n) for n in range(11)])

[1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796]

The Catalan numbers count a remarkable number of structures, all equivalent under bijections:

• Balanced strings of $n$ pairs of parentheses (“Generate Parentheses” produces exactly $C_n$ strings).
• Structurally distinct binary trees with $n$ nodes (“Unique Binary Search Trees”).
• Full binary trees with $n + 1$ leaves, equivalently the ways to fully parenthesize a product of $n + 1$ factors.
• Triangulations of a convex polygon with $n + 2$ sides.
• Monotone lattice paths from $(0,0)$ to $(n,n)$ that never go above the diagonal.
• Permutations of $1..n$ that can be sorted by one stack.

The recurrence says it directly for trees: choose the root, then $i$ nodes go left and $n - 1 - i$ go right. Both the closed form and the recurrence were checked against brute-force counts of balanced strings and of binary trees. Growth is roughly $4^n / (n^{3/2}\sqrt{\pi})$, so $C_{36}$ already exceeds $2^{63}$. Modulo a prime, use the Binomial class: C(2n, n) * inverse(n + 1).

Fibonacci numbers

$F_0 = 0$, $F_1 = 1$, $F_n = F_{n-1} + F_{n-2}$. They count tilings of a $1 \times n$ strip with squares and dominoes ($F_{n+1}$ ways), which is the same as “Climbing Stairs” with steps of 1 or 2, and binary strings of length $n$ with no two adjacent 1s ($F_{n+2}$). Binet’s formula $F_n = (\varphi^n - \psi^n)/\sqrt{5}$ with $\varphi = (1+\sqrt5)/2$ is exact in real arithmetic but, computed in floats, gives wrong integers from about $n = 71$ on. For exact or modular values use the iterative loop for moderate $n$ and matrix exponentiation or fast doubling for huge $n$. The sequence modulo $m$ is periodic (the Pisano period), which sometimes helps for $F_n \bmod m$ with small $m$. The exponential naive recursion and its memoized fix are in the Big-O article.

Probability and sampling

Reservoir sampling

The problem: choose $k$ items uniformly at random from a stream whose length is unknown in advance and too large to store. Examples are sampling lines from a huge log or “Linked List Random Node”.

Algorithm R: keep the first $k$ items. For the $i$-th item after that (0-indexed), draw $j$ uniformly from $0..i$; if $j < k$, replace slot $j$ with the new item.

import random

def reservoir_sample(stream, k, rng=random):
    reservoir = []
    for i, item in enumerate(stream):
        if i < k:
            reservoir.append(item)
        else:
            j = rng.randrange(i + 1)
            if j < k:
                reservoir[j] = item
    return reservoir

Why it is uniform. Claim: after processing $n$ items, each is in the reservoir with probability $k/n$. Item number $i$ (1-indexed, $i > k$) enters with probability $k/i$. It survives step $m > i$ unless that step picks exactly its slot, which happens with probability $\frac{k}{m} \cdot \frac{1}{k} = \frac{1}{m}$. So

$P(\text{item } i \text{ present at the end}) = \frac{k}{i} \prod_{m=i+1}^{n} \frac{m-1}{m} = \frac{k}{i} \cdot \frac{i}{n} = \frac{k}{n}$

The product telescopes. The first $k$ items start with probability 1 and the same product with $i = k$ gives $k/n$.

Empirical check: 200,000 runs choosing 3 of 10 items, counting how often each item was chosen (expected frequency 0.3):

from collections import Counter

rng = random.Random(42)
trials, n, k = 200_000, 10, 3
counts = Counter()
for _ in range(trials):
    counts.update(reservoir_sample(range(n), k, rng))
expected = trials * k / n
chi2 = sum((counts[i] - expected) ** 2 / expected for i in range(n))
print([round(counts[i] / trials, 4) for i in range(n)])
print(f"chi-square = {chi2:.2f} (9 degrees of freedom, 5% critical value 16.92)")

[0.2986, 0.3023, 0.3001, 0.3011, 0.2998, 0.3005, 0.2983, 0.2989, 0.3002, 0.3002]
chi-square = 8.38 (9 degrees of freedom, 5% critical value 16.92)

(The inclusion counts are not independent, since each run picks exactly three items, so the chi-square figure is a rough sanity check rather than a formal test. Uniform marginals are also necessary but not sufficient; Algorithm R does in fact make every $k$-subset equally likely.)

$O(n)$ time with one random number per item, $O(k)$ space. When random number generation is the cost, Algorithm L computes how many items to skip before the next replacement, making the expected number of random draws $O(k(1 + \log(n/k)))$. For a weighted reservoir, the Efraimidis-Spirakis method keeps the $k$ items with the largest $u^{1/w}$ (with $u$ uniform in $(0,1)$ and $w$ the weight), using a heap. Distributed variants sample each partition and merge, which is how sampling over sharded data works.

If the length is known and the data is in memory, use random.sample(population, k) instead.

Random Pick with Weight

Given weights $w_0, \dots, w_{n-1}$, return index $i$ with probability $w_i / \sum w$, many times. Build prefix sums once; each weight then owns an interval of the number line whose length equals the weight. Draw a uniform point in $[0, \text{total})$ and binary search for the interval:

weights   1     3        0      6
prefix    1     4        4      10
interval [0,1) [1,4)    empty  [4,10)

import bisect
from itertools import accumulate

class WeightedPicker:
    def __init__(self, weights):
        self.prefix = list(accumulate(weights))
        self.total = self.prefix[-1]

    def pick(self, rng=random):
        target = rng.random() * self.total
        return bisect.bisect_right(self.prefix, target)

rng = random.Random(7)
weights = [1, 3, 0, 6]
picker = WeightedPicker(weights)
print(picker.prefix)
counts = Counter(picker.pick(rng) for _ in range(100_000))
print([round(counts[i] / 100_000, 4) for i in range(len(weights))])

[1, 4, 4, 10]
[0.0996, 0.3023, 0.0, 0.5981]

$O(n)$ build, $O(\log n)$ per pick, $O(n)$ space. The choice of bisect_right matters. A target in $[1, 4)$ must map to index 1, and an exact hit on a boundary such as $4.0$ belongs to the interval starting there. bisect_right returns the first index whose prefix is strictly greater than the target, which gets both cases right and never selects a zero-weight entry. bisect_left would send $4.0$ to index 1 and $1.0$ to index 0, both wrong. The common integer formulation, a target drawn from $1..\text{total}$ with bisect_left, is the same thing shifted by one and is equally correct.

The standard library does this already: random.choices(population, weights, k=...) builds cumulative weights and bisects. For $O(1)$ per pick after $O(n)$ setup, the alias method (Walker, with Vose’s numerically stable construction) splits the distribution into $n$ equal-probability columns of at most two outcomes each. It wins when the number of picks is enormous and the weights are fixed. If weights change between picks, a Fenwick tree over the weights supports both updates and sampling in $O(\log n)$.

Fisher-Yates shuffle

To shuffle an array uniformly, so that each of the $n!$ permutations is equally likely: for $i$ from $n - 1$ down to 1, swap position $i$ with a uniformly random position in $0..i$. Position $i$ is then fixed, chosen uniformly from the items that remain. The method was published by Fisher and Yates in 1938 for pencil-and-paper use; the in-place computer version is Durstenfeld’s (1964), and Knuth’s The Art of Computer Programming popularized it.

def fisher_yates(items, rng=random):
    a = list(items)
    for i in range(len(a) - 1, 0, -1):
        j = rng.randrange(i + 1)          # 0 <= j <= i, j == i allowed
        a[i], a[j] = a[j], a[i]
    return a

It makes $n(n-1)\cdots 2 = n!$ equally likely choice sequences, each giving a distinct permutation, so the distribution is exactly uniform. $O(n)$ time, in place. The two classic mistakes are excluding $j = i$ (which yields Sattolo’s algorithm, producing only single-cycle permutations, $(n-1)!$ of them) and drawing $j$ from the whole array at every step. The second looks plausible but makes $n^n$ equally likely sequences, and $n^n$ is not divisible by $n!$ for $n \ge 3$, so it cannot be uniform. Shuffling "abc" 600,000 times each way:

def naive_shuffle(items, rng):
    a = list(items)
    for i in range(len(a)):
        j = rng.randrange(len(a))      # wrong: draws from the whole array
        a[i], a[j] = a[j], a[i]
    return a

rng = random.Random(3)
trials = 600_000
good = Counter(tuple(fisher_yates("abc", rng)) for _ in range(trials))
bad = Counter(tuple(naive_shuffle("abc", rng)) for _ in range(trials))
for perm in sorted(good):
    print("".join(perm), f"{good[perm] / trials:.4f}", f"{bad[perm] / trials:.4f}")

abc 0.1664 0.1473
acb 0.1667 0.1847
bac 0.1662 0.1852
bca 0.1666 0.1863
cab 0.1670 0.1480
cba 0.1670 0.1485

Fisher-Yates gives $1/6 \approx 0.1667$ for every permutation. The naive version produces 27 equally likely outcomes, distributed 4, 5, 5, 5, 4, 4 over the six permutations: $4/27 \approx 0.148$ and $5/27 \approx 0.185$, exactly the measured skew.

random.shuffle implements Fisher-Yates. One limit applies to any shuffle driven by a pseudo-random generator: the generator has finitely many states, so it can produce at most that many permutations. Python’s Mersenne Twister has period $2^{19937} - 1$, and $2081!$ already exceeds $2^{19937}$, so for longer sequences most permutations can never be produced. That does not matter for simulations; it does matter for, say, card shuffling with security requirements, where random.SystemRandom().shuffle (backed by the operating system’s cryptographic source) is the right call.

Computational geometry

Geometry problems in practice reduce to a few primitives on points, and almost all of them are built on one: the cross product.

Orientation via the cross product

For points $O$, $A$, $B$, the 2D cross product of the vectors $OA$ and $OB$ is

$\text{cross}(O, A, B) = (A_x - O_x)(B_y - O_y) - (A_y - O_y)(B_x - O_x)$

It equals twice the signed area of triangle $OAB$. Its sign gives the orientation:

cross > 0          cross < 0          cross = 0
counterclockwise   clockwise          collinear
(left turn)        (right turn)

      B             A
     /               \
    /                 \               O ---- A ---- B
   O ---- A            O ---- B

def cross(o, a, b):
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

With integer coordinates this is exact, with no division, no square roots and no angles. That is the main reason to prefer it over computing slopes (which divide, and fail for vertical lines) or angles with atan2 (which round). Uses: the signed area of a polygon (the shoelace formula sums cross over consecutive vertices), whether a polygon is convex (all turns the same sign), which side of a line a point lies on, and the two algorithms below.

Segment intersection

Two segments $P_1P_2$ and $P_3P_4$ intersect in the general case exactly when each straddles the line through the other: $P_1$ and $P_2$ lie on opposite sides of line $P_3P_4$, and $P_3$ and $P_4$ lie on opposite sides of line $P_1P_2$. If some orientation is 0, a point is collinear with the other segment, and the segments touch only if that point also lies within the other segment’s bounding box.

def on_segment(p, q, r):
    """q is collinear with p-r: is it between them?"""
    return (min(p[0], r[0]) <= q[0] <= max(p[0], r[0]) and
            min(p[1], r[1]) <= q[1] <= max(p[1], r[1]))

def segments_intersect(p1, p2, p3, p4):
    d1, d2 = cross(p3, p4, p1), cross(p3, p4, p2)
    d3, d4 = cross(p1, p2, p3), cross(p1, p2, p4)
    if ((d1 > 0) != (d2 > 0) and d1 != 0 and d2 != 0 and
            (d3 > 0) != (d4 > 0) and d3 != 0 and d4 != 0):
        return True                                   # proper crossing
    if d1 == 0 and on_segment(p3, p1, p4): return True
    if d2 == 0 and on_segment(p3, p2, p4): return True
    if d3 == 0 and on_segment(p1, p3, p2): return True
    if d4 == 0 and on_segment(p1, p4, p2): return True
    return False

print(segments_intersect((0, 0), (4, 4), (0, 4), (4, 0)),    # X shape
      segments_intersect((0, 0), (2, 2), (3, 3), (5, 5)),    # collinear, disjoint
      segments_intersect((0, 0), (3, 3), (3, 3), (5, 1)),    # share an endpoint
      segments_intersect((0, 0), (1, 1), (1, 0), (2, -1)))   # lines cross, segments don't

True False True False

The second case is the one that naive implementations get wrong: all four orientations are 0, and without the bounding-box check two collinear but separate segments would be reported as intersecting. The function was tested against an exact rational-arithmetic reference on 30,000 random segment pairs on a small integer grid, where degenerate cases (shared endpoints, collinear overlaps, zero-length segments) are frequent.

To test many segments for any intersection, the Shamos-Hoey sweep line finds whether any pair intersects in $O(n \log n)$, and Bentley-Ottmann reports all $k$ intersections in $O((n + k) \log n)$.

Convex hull: Andrew’s monotone chain

The convex hull of a point set is the smallest convex polygon containing all of them, the shape a rubber band takes around nails. Andrew’s monotone chain (A. M. Andrew, 1979) sorts the points by $x$ (then $y$) and builds the lower hull left to right and the upper hull right to left. The invariant: the chain being built always makes strict left turns. Before appending a point, pop from the chain while the last two points and the new one do not turn left.

def convex_hull(points):
    pts = sorted(set(points))
    if len(pts) <= 2:
        return pts
    lower = []
    for p in pts:
        while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
            lower.pop()
        lower.append(p)
    upper = []
    for p in reversed(pts):
        while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
            upper.pop()
        upper.append(p)
    return lower[:-1] + upper[:-1]        # last point of each is the other's first

print(convex_hull([(0, 0), (2, 0), (1, 1), (2, 2), (0, 2), (1, 0), (1, 2)]))

[(0, 0), (2, 0), (2, 2), (0, 2)]

Walk-through on those points, sorted: $(0,0), (0,2), (1,0), (1,1), (1,2), (2,0), (2,2)$.

lower chain:
  add (0,0)                       [(0,0)]
  add (0,2)                       [(0,0),(0,2)]
  (1,0): (0,0)->(0,2)->(1,0) turns right, pop (0,2)    [(0,0),(1,0)]
  (1,1): left turn                [(0,0),(1,0),(1,1)]
  (1,2): (1,0)->(1,1)->(1,2) collinear, pop (1,1)      [(0,0),(1,0),(1,2)]
  (2,0): right turn, pop (1,2); (0,0)->(1,0)->(2,0) collinear, pop (1,0)
                                  [(0,0),(2,0)]
  (2,2): left turn                [(0,0),(2,0),(2,2)]
upper chain (right to left) ends as [(2,2),(0,2),(0,0)]
hull = lower[:-1] + upper[:-1] = (0,0) (2,0) (2,2) (0,2), counterclockwise

Complexity. $O(n \log n)$ for the sort; the two passes are $O(n)$ amortized, since each point is pushed and popped at most once per chain. $O(n)$ space. If the points arrive sorted, the hull is linear.

Collinear points. With <= 0 in the pop condition, points lying in the middle of a hull edge are dropped, and the output contains only true corners. Changing it to < 0 keeps them, but then fully collinear inputs need special handling, because the lower and upper chains both contain every point. Duplicates are removed by set up front; repeated points otherwise produce zero cross products in unexpected places. The function was checked on 3,000 random point sets against a brute-force hull (a pair of points is a hull edge exactly when every other point lies strictly to its left or on the segment between them), and each output was confirmed to make strict counterclockwise turns.

Monotone chain is usually preferred over Graham scan (1972), which sorts by polar angle around the lowest point and needs care with ties and collinear points at the start. Both are $O(n \log n)$, which is optimal for comparison-based algorithms since sorting reduces to computing a hull. Jarvis march (gift wrapping) is $O(nh)$ for $h$ hull vertices, better only when $h$ is very small. Classic problems: “Erect the Fence” (which requires keeping collinear boundary points, so it uses < 0), the farthest pair of points (rotating calipers over the hull, $O(n)$ after the hull), and the smallest enclosing rectangle.

Precision pitfalls with floats

Everything above is exact with integer coordinates, and Python’s integers make overflow a non-issue. With floating-point coordinates, orientation can come out with the wrong sign:

from fractions import Fraction

print(0.1 + 0.2 == 0.3, 0.1 + 0.2, math.isclose(0.1 + 0.2, 0.3))

pts = [(0.3, 0.9), (0.5, 1.5), (0.4, 1.2)]           # all on y = 3x
print(cross(*pts))                                   # float: says "clockwise"
exact = [tuple(Fraction(str(c)) for c in p) for p in pts]
print(cross(*exact))                                 # exact: collinear

EPS = 1e-9
def orientation(o, a, b):
    c = cross(o, a, b)
    return 0 if abs(c) < EPS else (1 if c > 0 else -1)
print(orientation(*pts))