The Art of the Algorithm

SECTION 0What Is an Algorithm?

An algorithm is a finite, well-defined sequence of computational steps that transforms an input into an output. The word honors the 9th-century mathematician al-Khwārizmī, but the idea is older than the name and far older than the computer.

When we analyze an algorithm we ask two separate questions, and we must never conflate them:

• Correctness — does it always produce the right output? (Section 3.)
• Efficiency — how do its time and memory requirements grow as the input grows? (Sections 1–2.)

The model of computation

To count "steps" meaningfully we adopt the RAM (Random-Access Machine) model. We assume each basic operation — an arithmetic op, a comparison, a memory read or write of a single word — costs one unit of time. We then count operations as a function of input size n. This abstraction lets us compare algorithms independently of hardware, compiler, or programming language: a fact that will outlive any particular CPU.

% Your first runnable example: sum 1..n two ways and compare cost.
% The closed form is O(1); the loop is O(n). Same answer, different work.

function s = sum_loop(n)        % O(n) time, O(1) space
  s = 0;
  for k = 1:n
    s = s + k;
  endfor
endfunction

function s = sum_formula(n)     % O(1) time -- Gauss's trick
  s = n * (n + 1) / 2;
endfunction

% Try it:  sum_loop(100), sum_formula(100)  -> both give 5050

SECTION 1Asymptotic Analysis

We don't care that one machine is twice as fast as another, nor about the cost of a single addition. We care about shape: how the running time scales as n → ∞. Asymptotic notation is the language for that shape.

Proof technique: produce the witnesses

To prove a Big-O claim you must exhibit a constant \(c\) and threshold \(n_0\) that work. Here is the canonical example.

The same function is also \(O(n^3)\) (true but loose) and \(\Omega(n^2)\) (since \(3n^2+5n+2 \ge 3n^2\) for all \(n\ge 1\), take \(c=3\)). Combining the two tight bounds gives \(\Theta(n^2)\) — the honest answer.

Proving a negative bound

The hierarchy you must memorize

The gap between \(n\log n\) and \(n^2\) is the difference between an algorithm finishing in a blink and one that never finishes at all. This is why analysis matters more than micro-optimization.

SECTION 2Recurrences & Divide-and-Conquer

Divide-and-conquer algorithms split a problem into smaller subproblems, solve those recursively, and combine the results. Their running time is naturally described by a recurrence relation — an equation that defines \(T(n)\) in terms of itself at smaller inputs.

A generic divide-and-conquer recurrence has the form

\( T(n) = a\,T(n/b) + f(n), \)

where \(a\ge 1\) subproblems are each of size \(n/b\) (with \(b>1\)), and \(f(n)\) is the cost to divide the problem and combine the answers. Three methods solve such recurrences.

Method 1 — Substitution (guess & verify by induction)

Method 2 — The recursion tree

Unroll the recurrence as a tree: the root does \(f(n)\) work, its \(a\) children each do \(f(n/b)\), and so on. Summing work level by level often reveals the answer by inspection. For \(T(n)=2T(n/2)+n\): each of the \(\log_2 n + 1\) levels contributes exactly \(n\) total work (level \(i\) has \(2^i\) nodes each costing \(n/2^i\)), giving \(n\cdot(\log_2 n + 1) = \Theta(n\log n)\). The tree picture is also the fastest way to guess the bound the substitution method then proves.

Method 3 — The Master Theorem

For the common shape \(T(n)=aT(n/b)+f(n)\), one theorem dispatches most cases mechanically. The pivot is comparing \(f(n)\) against the "watershed" function \(n^{\log_b a}\).

Worked example — merge sort. \(T(n)=2T(n/2)+\Theta(n)\). Here \(a=2,b=2\), so \(c^\star=\log_2 2 = 1\) and \(f(n)=\Theta(n)=\Theta(n^{1})=\Theta(n^{c^\star})\). That is Case 2, giving \(T(n)=\Theta(n\log n)\) — matching the substitution proof above.

Worked example — binary search. \(T(n)=T(n/2)+\Theta(1)\). Now \(a=1,b=2\), \(c^\star=\log_2 1 = 0\), and \(f(n)=\Theta(1)=\Theta(n^0)\): Case 2 again, \(T(n)=\Theta(\log n)\).

SECTION 3Correctness & Loop Invariants

A fast algorithm that returns wrong answers is worthless. The workhorse tool for proving an iterative algorithm correct is the loop invariant — a property that is true before the loop starts and stays true after every iteration, like an inductive hypothesis embedded in the code.

Case study: Insertion Sort

Insertion sort builds a sorted prefix one element at a time, exactly the way most people sort a hand of playing cards. Here it is in Octave.

function A = insertion_sort(A)
  % Sorts row vector A ascending, in place. Time: O(n^2) worst, O(n) best.
  for j = 2:numel(A)
    key = A(j);            % the element to insert into the sorted prefix
    i = j - 1;
    % shift everything larger than key one slot to the right
    while i >= 1 && A(i) > key
      A(i+1) = A(i);
      i = i - 1;
    endwhile
    A(i+1) = key;          % drop key into the gap
  endfor
endfunction

SECTION 4Sorting

Sorting is the most studied problem in computing — partly because it is everywhere, and partly because it is the perfect arena for learning every analysis technique at once. We meet the quadratic family, then the \(\Theta(n\log n)\) family, then prove that no comparison sort can ever beat \(n\log n\).

The quadratic family

Three classic \(O(n^2)\) sorts. They are slow on large inputs but simple, in-place, and instructive.

function A = bubble_sort(A)
  % Repeatedly swap adjacent out-of-order pairs; large items "bubble" up.
  n = numel(A);
  for i = 1:n-1
    swapped = false;
    for j = 1:n-i                 % last i items are already in place
      if A(j) > A(j+1)
        tmp = A(j); A(j) = A(j+1); A(j+1) = tmp;
        swapped = true;
      endif
    endfor
    if !swapped, break; endif      % early exit if a full pass made no swap
  endfor
endfunction

function A = selection_sort(A)
  % Each pass selects the minimum of the unsorted tail and places it.
  n = numel(A);
  for i = 1:n-1
    m = i;
    for j = i+1:n
      if A(j) < A(m), m = j; endif
    endfor
    tmp = A(i); A(i) = A(m); A(m) = tmp;
  endfor
endfunction

Merge sort — divide and conquer, provably \(\Theta(n\log n)\)

Split the array in half, sort each half recursively, then merge the two sorted halves in linear time. The merge step is the whole trick: comparing the fronts of two sorted lists lets you emit the global minimum in \(O(1)\).

function A = merge_sort(A)
  n = numel(A);
  if n <= 1, return; endif          % base case: 0 or 1 element is sorted
  mid = floor(n/2);
  L = merge_sort(A(1:mid));
  R = merge_sort(A(mid+1:end));
  A = merge(L, R);
endfunction

function M = merge(L, R)
  % Combine two sorted vectors into one sorted vector in O(|L|+|R|).
  M = zeros(1, numel(L) + numel(R));
  i = 1; j = 1; k = 1;
  while i <= numel(L) && j <= numel(R)
    if L(i) <= R(j)               % "<=" keeps the sort STABLE
      M(k) = L(i); i = i + 1;
    else
      M(k) = R(j); j = j + 1;
    endif
    k = k + 1;
  endwhile
  while i <= numel(L), M(k) = L(i); i++; k++; endwhile   % drain L
  while j <= numel(R), M(k) = R(j); j++; k++; endwhile   % drain R
endfunction

Quicksort — fast in practice, fragile in the worst case

Pick a pivot, partition the array so smaller elements go left and larger go right, then recurse on each side. No merge needed — the work is all in partitioning.

function A = quick_sort(A, lo, hi)
  if nargin < 2, lo = 1; hi = numel(A); endif   % allow quick_sort(A)
  if lo < hi
    [A, p] = partition(A, lo, hi);   % p is the pivot's final resting index
    A = quick_sort(A, lo, p - 1);
    A = quick_sort(A, p + 1, hi);
  endif
endfunction

function [A, i] = partition(A, lo, hi)
  pivot = A(hi);                   % Lomuto scheme: last element is the pivot
  i = lo - 1;
  for j = lo:hi-1
    if A(j) <= pivot
      i = i + 1;
      tmp = A(i); A(i) = A(j); A(j) = tmp;
    endif
  endfor
  i = i + 1;
  tmp = A(i); A(i) = A(hi); A(hi) = tmp;   % put pivot between the two sides
endfunction

The fundamental limit: why \(n\log n\) is a wall

Every sort above that works only by comparing elements is a comparison sort. Remarkably, we can prove that no comparison sort — present or future, however clever — can beat \(n\log n\) in the worst case.

Merge sort and heap sort match this bound, so they are asymptotically optimal among comparison sorts. (Counting sort and radix sort can beat it — but only by abandoning comparisons and exploiting structure in the keys, so the theorem does not apply to them.)

SECTION 5Searching

Finding an element in a collection. The lesson here is sharp: a tiny structural assumption — sortedness — collapses a linear scan into a logarithmic leap.

Linear search — \(O(n)\), no assumptions

function idx = linear_search(A, target)
  idx = -1;                       % sentinel: -1 means "not found"
  for i = 1:numel(A)
    if A(i) == target
      idx = i; return;            % first match wins
    endif
  endfor
endfunction

Binary search — \(O(\log n)\) on sorted data

Look at the middle. If it's the target, done. If the target is smaller, the answer can only be in the left half; if larger, the right half. Each step throws away half of what remains.

function idx = binary_search(A, target)
  % PRECONDITION: A is sorted ascending. Returns an index of target, or -1.
  lo = 1; hi = numel(A); idx = -1;
  while lo <= hi
    mid = floor((lo + hi) / 2);   % integer midpoint, no overflow in Octave
    if A(mid) == target
      idx = mid; return;
    elseif A(mid) < target
      lo = mid + 1;               % target must be to the right
    else
      hi = mid - 1;               % target must be to the left
    endif
  endwhile
endfunction

SECTION 6Graph Algorithms

A graph \(G=(V,E)\) is a set of vertices joined by edges — the right model for maps, networks, dependencies, and relationships of every kind. We store one as an adjacency matrix adj(u,v)=1 when an edge \(u\to v\) exists, which is wonderfully natural in Octave's matrix world.

Breadth-First Search — shortest paths in unweighted graphs

BFS explores in concentric rings: all vertices at distance 1, then distance 2, and so on, using a FIFO queue. It visits every vertex and edge once, so it runs in \(O(V+E)\).

function order = bfs(adj, start)
  % Returns vertices in the order BFS first reaches them from `start`.
  n = rows(adj);
  visited = false(1, n);
  order = [];
  queue = [start];
  visited(start) = true;
  while !isempty(queue)
    u = queue(1); queue(1) = [];   % dequeue front
    order(end+1) = u;
    for v = 1:n                    % enqueue unvisited neighbours
      if adj(u, v) && !visited(v)
        visited(v) = true;
        queue(end+1) = v;
      endif
    endfor
  endwhile
endfunction

Depth-First Search — go deep, then backtrack

DFS plunges as far as possible along each branch before backing up. It is the backbone of cycle detection, topological sorting, and connectivity. Also \(O(V+E)\).

function order = dfs(adj, start)
  n = rows(adj);
  visited = false(1, n);
  order = [];
  [order, ~] = dfs_visit(adj, start, visited, order);
endfunction

function [order, visited] = dfs_visit(adj, u, visited, order)
  visited(u) = true;
  order(end+1) = u;
  for v = 1:rows(adj)
    if adj(u, v) && !visited(v)
      [order, visited] = dfs_visit(adj, v, visited, order);  % recurse
    endif
  endfor
endfunction

Dijkstra's algorithm — shortest paths with weights

When edges carry non-negative weights, Dijkstra greedily finalizes the closest unfinished vertex, relaxing its outgoing edges. Here W(u,v) is the positive weight of edge \(u\to v\), and 0 means "no edge".

function dist = dijkstra(W, src)
  % W: n-by-n weight matrix, W(u,v) > 0 is an edge, 0 = no edge. Weights >= 0.
  n = rows(W);
  dist = inf(1, n);
  dist(src) = 0;
  visited = false(1, n);
  for iter = 1:n
    % pick the unvisited vertex with the smallest tentative distance
    u = -1; best = inf;
    for v = 1:n
      if !visited(v) && dist(v) < best
        best = dist(v); u = v;
      endif
    endfor
    if u == -1, break; endif        % remaining vertices unreachable
    visited(u) = true;
    for v = 1:n                     % relax edges out of u
      if W(u, v) > 0 && !visited(v) && dist(u) + W(u, v) < dist(v)
        dist(v) = dist(u) + W(u, v);
      endif
    endfor
  endfor
endfunction

SECTION 7Dynamic Programming

Dynamic programming (DP) solves a problem by combining solutions to overlapping subproblems, each solved once and cached. It applies precisely when a problem has two features.

The motivating example: Fibonacci

The naive recursion recomputes the same values exponentially many times — \(T(n)=T(n-1)+T(n-2)+\Theta(1)\), which is \(\Theta(\varphi^{\,n})\). Caching collapses it to \(\Theta(n)\).

function f = fib_naive(n)         % EXPONENTIAL: O(phi^n) -- do not use for big n
  if n <= 1
    f = n;
  else
    f = fib_naive(n-1) + fib_naive(n-2);
  endif
endfunction

function f = fib_dp(n)            % BOTTOM-UP: O(n) time, O(n) space
  memo = zeros(1, n+1);
  memo(1) = 0;                   % fib(0); Octave is 1-indexed so memo(k+1)=fib(k)
  if n >= 1, memo(2) = 1; endif  % fib(1)
  for i = 3:n+1
    memo(i) = memo(i-1) + memo(i-2);
  endfor
  f = memo(n+1);
endfunction

0/1 Knapsack — the canonical DP table

Given items with integer weights and values and a capacity \(W\), choose a subset of maximum value that fits. The recurrence considers each item either left out or taken:

\(\;dp[i][w] = \max\big(\,dp[i\!-\!1][w],\; v_i + dp[i\!-\!1][w-w_i]\,\big)\;\) (the second term only if \(w_i \le w\)).

function best = knapsack(weights, values, W)
  % 0/1 knapsack. Time O(n*W), space O(n*W). Column j stores capacity j-1.
  n = numel(weights);
  dp = zeros(n+1, W+1);          % dp(i+1, w+1) = best value, first i items, cap w
  for i = 1:n
    for w = 0:W
      if weights(i) <= w
        take  = dp(i, w - weights(i) + 1) + values(i);
        leave = dp(i, w + 1);
        dp(i+1, w+1) = max(leave, take);
      else
        dp(i+1, w+1) = dp(i, w+1);     % item i too heavy; cannot take it
      endif
    endfor
  endfor
  best = dp(n+1, W+1);
endfunction
% Example: knapsack([1 3 4 5], [1 4 5 7], 7)  ->  9

Longest Common Subsequence

The edit-distance cousin behind diff and DNA alignment. Match extends the diagonal; a mismatch takes the better of dropping one character from either string.

function len = lcs(X, Y)
  % Length of the longest common subsequence of vectors/strings X and Y. O(mn).
  m = numel(X); n = numel(Y);
  dp = zeros(m+1, n+1);
  for i = 1:m
    for j = 1:n
      if X(i) == Y(j)
        dp(i+1, j+1) = dp(i, j) + 1;                  % extend the match
      else
        dp(i+1, j+1) = max(dp(i, j+1), dp(i+1, j));   % drop from X or from Y
      endif
    endfor
  endfor
  len = dp(m+1, n+1);
endfunction
% Example: lcs('ABCBDAB', 'BDCAB')  ->  4   ("BCAB")

SECTION 8Greedy Algorithms

A greedy algorithm makes the locally optimal choice at each step and never reconsiders. When it works, it is simpler and faster than DP — but it works only when the problem has a special structure, and proving it requires care.

Activity selection

Given activities with start and finish times, schedule the most that don't overlap. The greedy rule — always take the activity that finishes earliest among those still compatible — is optimal.

function chosen = activity_select(starts, finishes)
  % Returns indices of a maximum set of mutually compatible activities.
  [finishes, idx] = sort(finishes);     % sort by finish time, keep original idx
  starts = starts(idx);
  chosen = [idx(1)];                     % always take the earliest-finishing one
  last_finish = finishes(1);
  for i = 2:numel(starts)
    if starts(i) >= last_finish          % compatible: starts after last finishes
      chosen(end+1) = idx(i);
      last_finish = finishes(i);
    endif
  endfor
endfunction

SECTION 9Hashing & Hash Tables

Binary search bought us \(O(\log n)\) by demanding sorted data. Hashing strikes a bolder bargain: give up order entirely, and lookups, inserts, and deletes all become \(O(1)\) — on average. It is the structure quietly powering nearly every dictionary, set, and database index you have ever touched.

Because we are squeezing a huge key space into \(m\) slots, two keys can map to the same bucket — a collision. The pigeonhole principle guarantees this once \(n > m\), so every real implementation needs a resolution strategy.

A chaining hash table in Octave

Octave ships a real map as containers.Map, but building one by hand makes the mechanics concrete. We use a cell array of buckets, each bucket a matrix of [key, value] rows.

function T = ht_new(m)
  % create a hash table with m empty buckets
  T.m = m;
  T.buckets = cell(1, m);          % each cell holds a k-by-2 matrix of [key value] rows
endfunction

function h = ht_hash(T, key)
  % division method, shifted into Octave's 1-based index range
  h = mod(key, T.m) + 1;
endfunction

function T = ht_put(T, key, value)
  h = ht_hash(T, key);
  B = T.buckets{h};
  if ~isempty(B)
    row = find(B(:,1) == key, 1);  % key already stored in this bucket?
    if ~isempty(row)
      B(row, 2) = value;           % update in place
      T.buckets{h} = B;
      return;
    endif
  endif
  T.buckets{h} = [B; key, value];  % otherwise chain the new pair on
endfunction

function [found, value] = ht_get(T, key)
  h = ht_hash(T, key);
  B = T.buckets{h};
  found = false; value = nan;
  if ~isempty(B)
    row = find(B(:,1) == key, 1);  % scan only this one chain
    if ~isempty(row)
      found = true; value = B(row, 2);
    endif
  endif
endfunction

SECTION 10Amortized Analysis

Some operations are usually cheap but occasionally expensive. Amortized analysis asks the honest question: averaged over a worst-case sequence of operations, what does each one cost? No probability is involved — the averaging is over time, not over random inputs.

The motivating example: a growing array

Octave's convenient A(end+1) = x hides a trap. If the array is reallocated and fully copied on every append, \(n\) appends cost \(1+2+\dots+n=\Theta(n^2)\). The cure is to double the capacity whenever the array fills, so copying happens only rarely.

function S = darray_new()
  S.data = zeros(1, 1);            % backing store with spare capacity
  S.cap  = 1;                      % current capacity
  S.n    = 0;                      % number of elements actually used
endfunction

function S = darray_push(S, x)
  if S.n == S.cap                  % full: grow by doubling
    S.cap = 2 * S.cap;
    bigger = zeros(1, S.cap);
    bigger(1:S.n) = S.data(1:S.n); % copy old contents -- the costly step
    S.data = bigger;
  endif
  S.n = S.n + 1;
  S.data(S.n) = x;                 % cheap: write into existing capacity
endfunction

SECTION 11P, NP & the Edge of the Map

Some problems we can solve quickly; some we can only check quickly; and for a vast, important class, nobody knows which it is. This is the deepest open question in computer science.

Clearly \(P \subseteq NP\) (if you can solve it fast, you can verify it fast). Whether \(P = NP\) — whether every quickly-checkable problem is also quickly-solvable — is unsolved, and carries a \$1,000,000 Clay Millennium Prize. The hardest problems in NP are called NP-complete: the traveling salesman, Boolean satisfiability, graph coloring, and thousands more. A polynomial algorithm for any one of them would solve them all.

SECTION 12Quick Reference

A pocket card for the road. Octave idioms on the left, complexity facts on the right.

Octave survival kit

The complexity cheat sheet

SECTION 13Practice Exercises

Theory sticks when you wrestle with it. Try each problem with pencil first; the worked solution is one click away. They span the whole guide — asymptotics, recurrences, invariants, hashing, amortization, and graphs.

m = A(1);
for i = 2:numel(A)
  if A(i) > m
    m = A(i);
  endif
endfor