Numerical algorithms

Summary

The previous lessons of the unit taught how to write loops and how to reason about them with an invariant and a variant; this lesson applies that reasoning to four classic numerical algorithms, older than computers and still in use. The primality test by trial division, which stops at $\sqrt{n}$ for a reason that is proved, and whose stopping condition overflows if it is written carelessly. The Euclidean algorithm for the greatest common divisor, with its invariant, its termination, its trap with negative numbers and its comparison with std::gcd. Change of base by successive divisions and Horner’s rule, which rebuilds a number from its digits and evaluates a polynomial with the same scheme. And bisection, which approximates a root by splitting an interval in two, guaranteed by the intermediate value theorem, with a number of iterations that is computed in advance. Each algorithm is presented as a mathematical object: with a proof that it does what it says, a proof that it terminates and a first count of its steps.

Learning objectives

Implement the primality test by divisions up to $\sqrt{n}$ and justify why it suffices to reach $\sqrt{n}$ .
Implement the Euclidean algorithm, state its invariant (the greatest common divisor does not change) and compare it with std::gcd.
Convert a number to another base by successive divisions and evaluate a positional representation, or a polynomial, with Horner’s rule.
Approximate a root with bisection up to a tolerance, and compute how many iterations it requires.

Divisors and primes

Trial division

An integer $n \geq 2$ is prime if its only positive divisors are 1 and $n$ ; if it has another, it is composite; 0 and 1 are neither. The definition immediately suggests an algorithm: try each candidate $d$ between 2 and $n - 1$ and declare the number composite as soon as one divides it. It is the smallest-divisor loop of the lesson Loops, and it is correct. Now, for a ten-digit prime it performs almost ten billion divisions, and almost all of them are useless, for the following reason.

PropositionIt suffices to test up to

\sqrt{n}

If $n \geq 2$ is composite, it has a divisor $d$ with $2 \leq d \leq \sqrt{n}$ .

Proof

$n = a \cdot b, \quad 2 \leq a \leq b$
If $n$ is composite, it is written as a product of two factors, neither equal to 1 nor to $n$ . Call $a$ the smaller of the two.
$a^2 \leq a \cdot b = n$
Multiplying $a \leq b$ by the positive number $a$ gives $a^2 \leq ab = n$ , whence $a \leq \sqrt{n}$ : the smaller factor of every factorization lies to the left of $\sqrt{n}$ .
$d \leq \sqrt{n} \iff \frac{n}{d} \geq \sqrt{n}$
The divisors therefore come in pairs $(d, n/d)$ mirrored around $\sqrt{n}$ : with $n = 36$ , the pairs are $(2, 18)$ , $(3, 12)$ , $(4, 9)$ and $(6, 6)$ . If no candidate up to $\sqrt{n}$ divides $n$ , none beyond it divides it either, because its partner would have appeared earlier.

It follows that the number of divisions drops from about $n$ to about $\sqrt{n}$ : for $n$ close to two billion, from two billion to forty-six thousand. The program uses the flag pattern of the lesson Designing correct loops: isPrime starts out true and the loop turns it off if it finds a divisor.

C++ program

#include <iostream>#include <print>int main(){    int n{};    std::cin >> n;    bool isPrime{n >= 2};    for (int d{2}; isPrime && d <= n / d; ++d) {        if (n % d == 0) {            isPrime = false;        }    }    if (isPrime) {        std::println("{} is prime", n);    } else {        std::println("{} is not prime", n);    }}

Program input

2147483647

Output

2147483647 is prime

Three decisions deserve comment. The flag is initialized with n >= 2, so that 0, 1 and the negative numbers are declared not prime without entering the loop: forgetting that case is the typical error, because with $n = 1$ the loop performs no iteration and the program would declare 1 prime. The condition of the for has two parts: the flag, which stops the loop at the first divisor, and the bound. And the bound is written d <= n / d, which, for positive integers, is equivalent to $d^2 \leq n$ (integer division truncates, but $d$ does not exceed the truncated quotient exactly when it does not exceed $n/d$ , because $d$ is an integer). The number in the input, $2^{31} - 1$ , is the largest int and is prime; the program decides it with $46\,339$ remainders, those of the candidates from 2 to $46\,340$ .

The loop invariant is: no integer from 2 to $d - 1$ divides $n$ , as long as the flag remains true. It holds on entry, when the interval from 2 to 1 is empty; it is maintained, because if $d$ does not divide $n$ the interval grows by a number that does not divide it; and on exit with the flag true, the invariant and the false condition, $d^2 > n$ , say together that no candidate up to $\sqrt{n}$ divides $n$ , which is what the proposition requires. The variant is the number of candidates still to be tested, that is, of integers between $d$ and $\sqrt{n}$ : it decreases by one per iteration and does not go below zero.

The condition that overflows

The natural way to write the bound is d * d <= n, and it is an error. With the previous input, the loop reaches $d = 46\,341$ , and then the product $46\,341^2 = 2\,147\,488\,281$ exceeds the maximum of int.

C++ program

#include <iostream>#include <print>int main(){    int n{};    std::cin >> n;    bool isPrime{n >= 2};    for (int d{2}; isPrime && d * d <= n; ++d) {        if (n % d == 0) {            isPrime = false;        }    }    std::println("Prime? {}", isPrime);}

Program input

2147483647

Output

main.cpp:9:33: runtime error: signed integer overflow: 46341 * 46341 cannot be represented in type 'int'

The program exited with code 1.

The undefined behavior sanitizer stops the program at the product. Without it, signed overflow is undefined behavior, and the standard does not say what happens. What is observed with GCC on Linux, with and without optimization, is instructive: the product wraps around and becomes negative or smaller than $n$ , and since no int exceeds the maximum, the condition d * d <= n never stops holding; the loop goes on testing candidates for about two billion iterations, until it reaches $d = n$ , which divides $n$ , and the program declares a prime number composite. An overflow error does not show up as an overflow, but as a false answer. The form d <= n / d cannot overflow, because it divides instead of multiplying; the lesson Bits and the math library offered another way out, computing the bound once with std::sqrt, but then the exactness of the bound depends on floating-point rounding, and integer division does not have that problem.

The lesson “Prime numbers and the fundamental theorem of arithmetic”, in the course Algebra and Geometry I, proves that every integer greater than 1 is a product of primes in exactly one way; trial division is the most direct way of finding that product, dividing $n$ by each divisor found and continuing with the quotient.