76. Algorithms

Section 75 showed that the choice of algorithm gains a factor of thousands where rewriting the code gains percent. This section is about how to estimate an algorithm before you write it, and a few techniques that come up all the time in microcontroller programs: binary search and lookup in a calibration table, fast exponentiation, calculations without overflow, a stack for parsing expressions, working with bits.

Every answer in this section was obtained by running a program and recomputed independently.

Terms used in the documentation

Term	Meaning
time, space complexity	how running time, memory use grow with the input
big-O notation	the notation for the growth order
worst, average case	the least favourable input, a typical input
binary search	search by halving a sorted range
lookup table, interpolation	a table of precomputed points, values between them
exponentiation by squaring	raising to a power in O(log n) multiplications
binomial coefficient	the number of ways to choose k of n
reverse Polish notation	the operator written after its operands
shunting-yard algorithm	converting ordinary notation to RPN
count leading zeros (CLZ)	the number of zero bits above the highest one
population count	the number of one bits

76.1 How an algorithm is estimated

The exact running time depends on the processor, the compiler and the data. That is why algorithms are compared by how the work grows as the amount of data grows, and this is written in big-O notation, keeping only the leading term. An algorithm that takes 10n⁵ + 2n steps is O(n⁵).

For a thousand elements:

Two caveats for the chip.

For small n the constants matter more than the order. For ten elements, bubble sort with its tiny code can turn out better than a "proper" sort that drags in noticeably more code. The order starts to decide when the data runs to hundreds and thousands.

Space complexity matters as much as time complexity. An algorithm that needs a second array of the same size or deep recursion may simply not fit on a chip with 20 KB of RAM (section 66.7).

76.2 Binary search

In a sorted array, the value you look for is found by halving: compare with the middle and discard the half where it certainly is not.

int bsearch_int(const int *a, int n, int key)
{
    int lo = 0, hi = n - 1;
    while (lo <= hi) {
        int mid = lo + (hi - lo) / 2;
        if (a[mid] == key) return mid;
        if (a[mid] < key)  lo = mid + 1;
        else               hi = mid - 1;
    }
    return -1;                         /* нет */
}

In an array of a million elements, the last one is found in 20 steps, and the same number was needed to establish that a number is absent. Linear search would spend up to a million.

The classic mistake is bounds without ±1: lo = mid and hi = mid instead of lo = mid + 1 and hi = mid - 1. In the array {10, 20, 30, 40, 50} this variant finds 20 on the second step, but searching for 50, the last element, gets stuck forever: the middle between 3 and 4 is 3, the bound moves to 3, and everything repeats. A single mistake lo = mid is enough: then the search for 20 gets stuck instead. Always test a search on the end elements and on a missing one.

And a subtlety: (lo + hi) / 2 overflows on huge arrays, while lo + (hi - lo) / 2 never does. On a chip the arrays are small, but the habit is useful. The ready-made function is bsearch from the standard library (section 74.7).

76.3 🔧 Calibration table

A typical task on a chip: the sensor characteristic is nonlinear and given as a table of points. The table goes into Flash (const, section 68.3), the right segment is found by binary search, and inside the segment it is linear interpolation in integers (section 69.9):

static const int16_t adc_tab[]  = { 400,  800, 1200, 1600, 2000, 2400, 2800, 3200, 3600};
static const int16_t temp_tab[] = {-400, -200,    0,  150,  300,  450,  600,  800, 1000}; /* 0,1 °C */

int16_t adc_to_temp(int adc)
{
    int lo = 0, hi = 8;
    if (adc <= adc_tab[lo]) return temp_tab[lo];      /* ниже таблицы */
    if (adc >= adc_tab[hi]) return temp_tab[hi];      /* выше таблицы */
    while (hi - lo > 1) {
        int mid = (lo + hi) / 2;
        if (adc_tab[mid] <= adc) lo = mid; else hi = mid;
    }
    return temp_tab[lo] + (temp_tab[hi] - temp_tab[lo]) * (adc - adc_tab[lo])
                        / (adc_tab[hi] - adc_tab[lo]);
}

The numbers in the table are made up, but the algorithm is checked: code 1000 gives −10.0 °C, 2000 gives 30.0 °C, 3300 gives 85.0 °C, and everything below the table is clamped to the lowest point, −40.0 °C. In a table of nine points the segment is found in exactly three steps, without a single floating-point operation.

76.4 Fast exponentiation

aⁿ can be computed with n multiplications, or by using the fact that aⁿ = (aⁿ/²)² for even n and aⁿ = aⁿ⁻¹ · a for odd n:

long binpow(long a, int n)
{
    if (n == 0) return 1;
    if (n & 1)  return binpow(a, n - 1) * a;
    long b = binpow(a, n >> 1);
    return b * b;
}

3¹³ = 1594323: 13 multiplications the direct way, 6 the fast way. For the power 1000 it is 15 multiplications instead of a thousand. The number of multiplications is the number of binary digits of the exponent minus one plus the number of ones in it: that is why the algorithm runs in O(log n).

The technique works for any operation where the order of parentheses does not matter. The most frequent use on a chip is modular exponentiation in cryptography (section 43).

76.5 Fibonacci in logarithmic time

Fibonacci numbers are obtained by raising the matrix [[0, 1], [1, 1]] to a power, which means by the same fast method:

Checked: F(10) = 55, F(90) = 2880067194370816120; the latter still fits in 64 bits. The lesson is wider than Fibonacci: the same problem is solved in O(2ⁿ), O(n) and O(log n), and the difference comes from the idea, not from the code.

76.6 Binomial coefficients and overflow

The number of combinations C(n, k) is how many ways there are to choose k items out of n. You cannot compute it with the factorial formula: the intermediate values are huge.

int cnk_naive(int n, int k)      /* 21 · 22 · … · 30, потом делить */
{
    int res = 1;
    for (int i = n - k + 1; i <= n; ++i) res *= i;
    for (int i = 2; i <= k; ++i)         res /= i;
    return res;
}

C(30, 10) computed this way gives −108 instead of 30045015: the product 21 · … · 30 does not fit in an int, and what gets divided after that is already garbage.

The right way is to multiply and divide in turn:

long cnk(int n, int k)
{
    long res = 1;
    for (int i = 1; i <= k; ++i)
        res = res * (n - k + i) / i;       /* деление всегда нацело */
    return res;
}

After the i-th step res holds C(n − k + i, i), an integer, and the division at each step is exact. The intermediate values stay close to the answer.

The long type gives headroom here only on the Linux computer where the example was checked: there it is 64-bit. On Cortex-M long is the same 32 bits as int (section 64.3). For C(30, 10) even those are enough: the largest intermediate value is 10015005 × 30 = 300450150, less than 2³¹. But C(34, 17) = 2333606220 does not fit in a 32-bit signed type even as the answer. On a chip, for headroom, you write uint64_t, both for res and for the result; binpow from 76.4 has the same issue.

A third way is Pascal's triangle, where each number is the sum of the two above it. It uses only additions, and the intermediate values never exceed the answer. Row 5: 1 5 10 10 5 1.

An example from chip life: how many bytes contain exactly three ones? C(8, 3) = 56, and a direct enumeration of all 256 bytes gives the same 56.

76.7 The stack at work: reverse Polish notation

In reverse Polish notation the operation is written after its operands: 1 2 + 4 * 3 + instead of (1 + 2) * 4 + 3. No parentheses are needed, and such an expression is evaluated in one pass with a stack (section 72.5): a number is pushed onto the stack, an operation pops two numbers, performs itself and pushes the result.

Ordinary notation is converted to Polish notation by the shunting-yard algorithm: numbers go straight to the output; an incoming operation first pops from the stack to the output all operations with a priority not lower than its own, and only then is pushed onto the stack. "(" is pushed onto the stack, ")" pops everything down to "(", at the end the stack is emptied to the output. Checked: (3+5)*12-17*3 → 3 5 + 12 * 17 3 * - → 45.

On a chip this means parsing commands and formulas received over a port, and simple built-in calculators for settings.

76.8 Bits as an algorithm

Counting zeros in binary notation. How many numbers from 1 to N have exactly K zeros? Enumeration with bit parsing: for N = 8, K = 1 it is 3 numbers (10, 101 and 110 in binary), for 13 and 2 it is 4, for 1000 and 5 it is 210.

Base-2 logarithm in one instruction. The index of the highest set bit is 31 − (number of leading zeros). On Cortex-M3 and M4 there is the clz instruction for this, and __builtin_clz turns into exactly that. On Cortex-M0 there is no such instruction, and the same calculation takes more than two dozen instructions. Population count has no dedicated instruction anywhere: on M3 __builtin_popcount is a dozen instructions.

The exclusive-OR cipher: each byte of the text is added modulo 2 to a byte of the key, and applying the same key again returns the text. But with the key HASSWORD and the text Hello world the first byte gives H ^ H = 0, and the encrypted "string" has length zero. Encrypted data is binary data, not a string: it is stored together with its length (section 68.8). And most importantly, this is not encryption: the key is recovered from a few bytes of known text. Real data protection is in section 43.

76.9 Choosing an algorithm on a chip

• The worst case matters more than the average. In a device with hard deadlines what matters is how long the most unlucky input takes: for binary search it is log n, for bubble sort n².
• In place, without a second array. Memory is more expensive than time.
• No deep recursion, or with a proven depth (section 66.7).
• Tables in Flash trade memory, of which there is a lot, for time, of which there is little: sines, calibrations, coefficients are computed in advance.
• Integers instead of floating point if there is no FPU (section 69.9).
• Simple while the data is small. A complex algorithm is justified when a measurement (section 73) has shown that the simple one does not keep up.

76.10 Pitfalls

Binary search hangs. The bounds move without ±1; test on the last element and on a missing one.

Binary search returns the wrong thing. The array is not sorted.

Combinations or factorials gave a negative number. Overflow of the intermediate product.

Exponentiation in a loop is slow. The fast method takes O(log n) multiplications.

Table interpolation jumps at the edges. Values outside the table are not handled.

Expression parsing breaks on parentheses. The operator stack is not bounded by the parentheses.

The encrypted string gets "cut off". A zero byte appeared in the result; store it with its length.

__builtin_clz(0) behaves unpredictably. For zero the result of the built-in function is undefined, and the compiler is entitled to assume that zero never comes. Zero is checked separately.

The logarithm is "slow" on F0. Cortex-M0 has no clz instruction.

The algorithm is fine on the computer and does not fit in the chip. The space complexity was not counted.

76.11 Summary

• Algorithms are compared by the growth of work: at n = 1000 it is 10 steps for O(log n), a thousand for O(n) and a million for O(n²).
• Binary search in a million elements takes 20 steps; the bounds must move by ±1, otherwise the search for the last element hangs.
• A calibration table in Flash, a binary search for the segment and integer interpolation are the standard way to linearize a sensor.
• Fast exponentiation: 3¹³ takes 6 multiplications instead of 13, the power 1000 takes 15.
• Fibonacci is solved in O(2ⁿ), O(n) and O(log n): the idea decides, not the code.
• C(30, 10) via a product in int gives −108; multiplying and dividing in turn gives 30045015.
• A stack evaluates reverse Polish notation, the shunting yard converts ordinary notation into it.
• clz gives the logarithm in one instruction on M3 and M4, but not on M0.
• XOR with a key is not encryption, and the result is binary data, not a string.

Next: section 77 "Compilation and compilers".