Algorithms & Combinatorics — Sorting, Graph Colouring and Constraint Satisfaction

Algorithms are recipes for solving problems, and visualising them changes understanding irrevocably. Watching bubble sort grind through a nearly-sorted array while quicksort finishes in a flash makes O(n log n) visceral. Seeing a backtracker place and retract queens across a chessboard makes NP-hardness intuitive. Seven interactive simulations explore discrete algorithms from numerical integration to combinatorial search.

Why Visualise Algorithms?

Computer science education traditionally presents algorithms as pseudocode and complexity tables — O(n²) here, O(n log n) there. This works for theoretical analysis but fails to build intuition for when algorithms struggle. A sorting visualiser shows not just that merge sort is faster than insertion sort, but why: insertion sort wastes work moving elements one position at a time across large unsorted regions, while merge sort always divides the problem in half.

Similarly, seeing the N-queens backtracker retreat from dead-end positions makes vivid what "exponential worst case" means — and makes the contrast with the DSatur graph colouring heuristic that runs in polynomial time all the more striking. This spotlight covers seven simulations that span from classical sorting to open mathematical conjectures.

QuickSort · MergeSort · HeapSort Backtracking Graph Colouring Riemann Sums Collatz Conjecture Taylor Series Pascal's Triangle

Layer 1: The Art of Ordering — Sorting Algorithms

Sorting Algorithms Visualiser

Sorting is the canonical introductory algorithm problem, yet its depth is surprising. There are dozens of general-purpose sorting algorithms, each with different best-case, worst-case, and average-case complexities, memory requirements, stability properties, and practical performance profiles. A visualiser that plays them side-by-side makes these abstract properties concrete within seconds.

Bubble sort is O(n²) and painfully obvious to watch: it scans through the array repeatedly, swapping adjacent out-of-order elements, leaving each pass with the largest unsorted element "bubbled" to its final position. Quicksort chooses a pivot, partitions the array around it, and recursively sorts each half — on average O(n log n) but O(n²) on sorted or nearly-sorted input (addressed by randomised pivot selection). Merge sort is always O(n log n) but uses O(n) extra memory. Radix sort breaks the O(n log n) barrier but only for integer keys.

Sorting Algorithm Complexity Comparison

Algorithm       Best      Average   Worst     Space  Stable?
─────────────────────────────────────────────────────────────
Bubble Sort     O(n)      O(n²)     O(n²)     O(1)   ✓
Insertion Sort  O(n)      O(n²)     O(n²)     O(1)   ✓
Selection Sort  O(n²)     O(n²)     O(n²)     O(1)   ✗
Shell Sort      O(n log n) varies  O(n^1.5)   O(1)   ✗
Merge Sort      O(n lgn)  O(n lgn)  O(n lgn)  O(n)   ✓
Heap Sort       O(n lgn)  O(n lgn)  O(n lgn)  O(1)   ✗
QuickSort       O(n lgn)  O(n lgn)  O(n²)     O(lgn) ✗
Tim Sort        O(n)      O(n lgn)  O(n lgn)  O(n)   ✓
Radix Sort      O(nk)     O(nk)     O(nk)     O(n+k) ✓
Counting Sort   O(n+k)    O(n+k)    O(n+k)    O(k)   ✓

Information-theoretic lower bound for comparison sorts:
  Any comparison-based sort requires Ω(n log n) comparisons
  Proof: n! possible orderings → need ≥ log₂(n!) comparisons
         Stirling: log₂(n!) ≈ n·log₂(n) − n·log₂(e)

Tim Sort (Python's sorted(), Java Arrays.sort):
  Hybrid merge + insertion sort
  Identifies natural runs (sorted subsequences)
  O(n) on nearly-sorted data — real-world optimal

The Sorting Algorithms simulation runs 15+ algorithms simultaneously on the same shuffled array, displaying a colour-coded bar chart where height encodes value and colour encodes comparison/swap status. The swap counter and timer update live, making the O(n²) vs O(n log n) divide immediately visible. An audio mode maps bar height to pitch — you can literally hear quicksort finish before merge sort.

Practical wisdom: No single sorting algorithm dominates in practice. TimSort wins on nearly-sorted data (like database query results). Radix sort wins for fixed-width integers. Quicksort's cache-friendliness makes it faster than merge sort for random data despite identical O(n log n) complexity. Always profile, never assume.

Layer 2: Placing Non-Attacking Queens — N-Queens

N-Queens Backtracking Simulation

The N-queens problem is to place N queens on an N×N chessboard such that no two queens attack each other — no two may share a row, column, or diagonal. For N=8, there are 92 solutions; for N=12, there are 14 200; for N=20, over 39 billion. Counting them all requires an exhaustive search, but finding a single solution quickly is achievable with heuristics.

Backtracking — the systematic exploration of a search tree with early termination when a partial solution is detected to be infeasible — is the canonical algorithm. Place a queen in row 1, try each column, place one in row 2, try each valid column, and so on. When no valid placement exists in the current row, backtrack to the previous row and try the next column. The search tree has N! leaves in the worst case, but pruning makes it manageable.

N-Queens — Backtracking & Solution Counts

Backtracking pseudocode:
  solve(row):
    if row == N: solution found; record()
    for col in 0..N-1:
      if safe(row, col):
        place(row, col)
        solve(row + 1)
        remove(row, col)      ← backtrack

Safety check O(1) with bitmask tracking:
  cols    = set of occupied columns
  diag1   = set of occupied (row − col) diagonals
  diag2   = set of occupied (row + col) anti-diagonals
  safe iff col ∉ cols ∧ (row−col) ∉ diag1 ∧ (row+col) ∉ diag2

Solution counts (total, all rotations/reflections included):
  N   Solutions    N   Solutions
  4      2         12   14 200
  5      10        14   365 596
  6      4         16   14 772 512
  7      40        20   39 029 188 884
  8      92
  10     724

For large N: Q(N) ≈ (0.143 N)^N  (exponential growth)

NP-complete variant:
  Determining if a partial placement can be completed is NP-complete
  (a decision version studied in AI constraint satisfaction)

The N-Queens simulation visualises each step of the backtracker: the active queen being placed glows yellow, successfully placed queens glow green, and queens in conflict glow red. A solutions counter and backtrack counter update with each step. Set N from 4 to 12 and choose the animation speed to follow the algorithm's logic at your own pace, or hit "Auto" to race through all solutions.

Layer 3: Colouring Graphs — Graph Colouring

Graph Colouring Simulation

Graph colouring assigns colours to graph vertices such that no two adjacent vertices share a colour. The minimum number of colours required is the chromatic number χ(G), one of the most studied parameters in combinatorics. It has direct applications in scheduling (assign time slots to conflicting tasks), register allocation (assign CPU registers to variables), map colouring (the famous four-colour theorem), and frequency assignment in wireless networks.

Finding χ(G) exactly is NP-hard — but several polynomial-time heuristics produce near-optimal colourings quickly. The greedy algorithm assigns the smallest available colour to each vertex in some order. Welsh-Powell improves on this by ordering vertices by degree (highest first). DSatur (Degree of Saturation) is a dynamic-priority greedy — at each step it colours the vertex with the highest number of differently-coloured neighbours, which often achieves the optimal colouring.

Graph Colouring — Bounds & Algorithms

Chromatic number bounds:
  χ(G) ≥ ω(G)          (clique number — largest complete subgraph)
  χ(G) ≤ Δ(G) + 1      (Brook's bound; equality only for complete graphs or odd cycles)
  χ(G) ≤ Δ(G)          (Brook's theorem for connected non-complete non-odd-cycle graphs)

Four Colour Theorem (Appel & Haken, 1976):
  Any planar graph satisfies χ(G) ≤ 4
  (First major theorem proved with computer assistance — 1200+ configurations)

DSatur algorithm:
  saturation(v) = |{colours of neighbours of v}|  (number of distinct neighbour colours)
  1. Select uncoloured vertex v with max saturation (ties broken by degree)
  2. Assign smallest colour not used by any neighbour of v
  3. Repeat until all vertices coloured
  Time complexity: O(n²) naive, O((n+m) log n) with priority queue

Greedy upper bound:
  Greedy orders vertices arbitrarily → χ ≤ Δ + 1
  Welsh-Powell orders by degree ↓ → tighter bound in practice
  Worst case greedy: χ = Δ + 1 (e.g. bipartite graph with alternating high-degree vertices)

Petersen graph:
  10 vertices, 15 edges, Δ = 3, χ = 3
  3-regular, bridgeless — counterexample to many graph conjectures

The Graph Colouring simulation shows Greedy, Welsh-Powell, and DSatur running on preset and custom graphs. A step-by-step mode highlights which vertex is being coloured and why — you can see DSatur dynamically reorder its priority queue as colours are assigned. Drag vertices to adjust the layout; click empty space to add a new vertex; conflicts are highlighted in red when manual assignments create them.

Layer 4: Areas Under Curves — Riemann Integration

Riemann Integral Simulation

Numerical integration — approximating the area under a curve when an analytical antiderivative is unavailable — is one of the oldest and most practically important computational tasks. Riemann sums divide the integration interval into rectangles whose heights are sampled at the left endpoint, right endpoint, or midpoint of each subinterval. Trapezoidal and Simpson's rules use linear and quadratic interpolation respectively to achieve higher accuracy with fewer function evaluations.

The accuracy of each method depends on the smoothness of the integrand. For smooth functions, Simpson's rule achieves O(h⁴) error — the same result that requires millions of Riemann rectangles can be achieved with just 8 Simpson intervals. This is why numerical libraries use adaptive quadrature (bisecting intervals where the integrand varies rapidly), not uniform Riemann sums, for production calculations.

Numerical Integration Rules — Error Orders

∫ₐᵇ f(x) dx  with n equal subintervals of width h = (b−a)/n

Left Riemann sum:
  L_n = h · Σᵢ₌₀^{n−1} f(xᵢ)
  Error: O(h) = O(1/n)  — first-order method

Right Riemann sum:
  R_n = h · Σᵢ₌₁^n f(xᵢ)
  Error: O(h) = O(1/n)  — first-order method

Midpoint rule:
  M_n = h · Σᵢ f(xᵢ + h/2)
  Error: O(h²) — second-order; twice better than L/R for smooth f

Trapezoidal rule:
  T_n = h · [f(x₀)/2 + f(x₁) + f(x₂) + ... + f(xₙ₋₁) + f(xₙ)/2]
  Error: O(h²) — same order as midpoint but larger constant

Simpson's 1/3 rule (composite):
  S_n = (h/3) · [f(x₀) + 4f(x₁) + 2f(x₂) + 4f(x₃) + ... + 4f(xₙ₋₁) + f(xₙ)]
  (requires n even)
  Error: O(h⁴) — fourth-order; 100× fewer evaluations than Riemann for same accuracy

Convergence comparison for ∫₀¹ sin(x) dx = 1 − cos(1) ≈ 0.459698:
  n=10:  L error 0.0450,  Midpoint 0.00046,  Simpson 4.5×10⁻⁷
  n=100: L error 0.0045,  Midpoint 4.6×10⁻⁶, Simpson 4.5×10⁻¹¹

The Riemann Integration simulation draws each numerical method's rectangles or trapezoids in a different colour, showing the filled area and the uncovered gaps side by side. The live error display shows how each rule converges as you increase the number of subintervals with a slider — making the O(h⁴) advantage of Simpson's rule dramatic: increase n from 10 to 20 and the Simpson error drops by a factor of 16 while the Riemann errors drop by only 2.

Layer 5: An Unproven Conjecture — The Collatz Problem

Collatz Conjecture Simulation

Take any positive integer. If it is even, divide it by 2. If it is odd, multiply by 3 and add 1. Repeat. The Collatz conjecture states that every positive integer eventually reaches 1. It has been verified for all integers up to approximately 2.95 × 10²⁰, yet no general proof exists. Paul Erdős said "Mathematics is not yet ready for such problems."

The sequences exhibit seemingly chaotic behaviour. Starting from 27 requires 111 steps before reaching 1, climbing as high as 9232. Starting from 871 reaches 190 996 before descending. The stopping time (number of steps to reach 1) varies wildly and unpredictably — yet the conjecture asserts it is always finite. The Collatz problem sits at the intersection of number theory, dynamical systems, and computational complexity.

Collatz Conjecture — Formal Statement & Analysis

The 3n+1 function:
  T(n) = n/2       if n ≡ 0 (mod 2)
  T(n) = (3n+1)/2  if n ≡ 1 (mod 2)  [accelerated form]

Stopping time σ(n): smallest k such that T^k(n) = 1
Total stopping time τ(n): number of iterations until first reaching 1

Extremal examples:
  n=27:  σ=111, max value 9232
  n=871: σ=178, max value 190 996
  n=6171: σ=261, max value 975 400
  n=77031: σ=350, max value 21 933 016

Heuristic analysis (Terras, 1976):
  Average T(n) ≈ n · (3/4)  [under random-parity assumption]
  Expected stopping time: O(log n)
  But worst-case stopping time is not bounded by any proven function of log n

Verified range (as of 2024):
  Collatz conjecture holds for all n ≤ 2.95 × 10²⁰
  (Barina, 2020, distributed computation on GPU cluster)

Connection to open problems:
  Terry Tao (2019) proved "almost all" orbits reach a value below any given f(n)
  for f slowly tending to infinity — closest to a proof to date

The Collatz Conjecture simulation offers three visualisations: a log-scale sequence plot showing the trajectory of a single starting number, a stopping-time heat map showing σ(n) for n = 1 to N with colour proportional to stopping time, and a reverse Collatz tree that shows which numbers flow into any given value. The dramatic spikes in the stopping-time heat map — islands of slow sequences surrounded by fast neighbours — make the conjecture's unpredictability tangible.

Why is it so hard? The Collatz function mixes additive and multiplicative operations on integers, jumping between different residue classes in a way that defeats standard number- theoretic techniques. It cannot be analysed by the usual tools of divisibility proofs, generating functions, or algebraic identities — it truly is in a no-man's-land between solved mathematics and computational intractability.

Layer 6: Infinite Sums as Functions — Taylor Series

Taylor Series Simulation

Taylor's theorem says that any infinitely differentiable function can be approximated arbitrarily well near a point by a polynomial — and in many cases the approximation converges to the exact function everywhere. The polynomial's coefficients are determined entirely by the function's derivatives at a single point: f(a), f'(a), f''(a), and so on. This is a profound observation: the entire behaviour of sin(x) across all real numbers is encoded in the single number sin(0) = 0 and its derivatives 1, 0, −1, 0, 1, 0, …

Taylor series are the engine of scientific computation. They turn transcendental functions into polynomials, enabling fast evaluation on hardware that only natively performs addition and multiplication. They explain why sin(x) ≈ x for small x (a physics approximation used everywhere from optics to orbital mechanics) and why e^x ≈ 1+x for small perturbations. They provide the theoretical foundation for numerical differentiation, automatic differentiation, and neural network activation functions.

Taylor Series — Definition & Convergence

Taylor series of f(x) around x = a:
  f(x) = Σₙ₌₀^∞ f⁽ⁿ⁾(a) / n! · (x−a)ⁿ
       = f(a) + f'(a)(x−a) + f''(a)(x−a)²/2! + f'''(a)(x−a)³/3! + ...

Maclaurin series (a = 0):
  sin(x)   = x − x³/3! + x⁵/5! − x⁷/7! + ...        R = ∞
  cos(x)   = 1 − x²/2! + x⁴/4! − x⁶/6! + ...        R = ∞
  eˣ       = 1 + x + x²/2! + x³/3! + x⁴/4! + ...     R = ∞
  ln(1+x)  = x − x²/2 + x³/3 − x⁴/4 + ...            R = 1
  (1+x)^α  = 1 + αx + α(α−1)x²/2! + ...               R = 1  (binomial)
  arctan(x) = x − x³/3 + x⁵/5 − x⁷/7 + ...           R = 1
              (→ Leibniz formula: π/4 = 1 − 1/3 + 1/5 − ...)

Radius of convergence R:
  Determined by nearest singularity in the complex plane
  ln(1+x): singularity at x = −1 → R = 1
  1/(1−x) = Σ xⁿ: singularity at x = 1 → R = 1
  sin(x): no finite singularities → R = ∞

Remainder bound (Lagrange):
  |f(x) − Tₙ(x)| ≤ M · |x−a|^{n+1} / (n+1)!
  where M = max |f^{(n+1)}(t)| for t between a and x

The Taylor Series simulation overlays the partial sums T₁(x), T₂(x), ..., Tₙ(x) on top of the original function f(x). A slider adds one term at a time, with each new term animating in over 600 ms. A separate error plot shows |f(x) − Tₙ(x)| on a logarithmic scale, making the exponential convergence visible. Watch the approximation expand outward from the expansion point like a wave, converging on different functions at different rates depending on their radius of convergence.

Layer 7: Numbers in a Triangle — Pascal's Triangle

Pascal's Triangle Simulation

Pascal's triangle is constructed by the rule that each entry is the sum of the two entries above it, with edges padded with 1s. The resulting triangle contains the binomial coefficients C(n,k) = n!/(k!(n−k)!) — counting the number of ways to choose k items from n — which appear as the coefficients in the expansion of (a+b)ⁿ. But the triangle hides far more than binomial coefficients.

The diagonal sum of row n + row n−1 diagonals gives the Fibonacci numbers. Colouring cells by their parity (odd/even) produces Sierpiński's triangle fractal. Colouring by residue mod p for prime p produces a p-adic fractal with self-similar structure at all scales. The sums of rows are powers of 2; the hockey-stick identity gives cumulative sums of columns; the central column contains the central binomial coefficients which appear in random walk return probabilities.

Pascal's Triangle — Hidden Identities

Basic definition:
  C(n, k) = C(n−1, k−1) + C(n−1, k)    (Pascal's rule)
  C(n, 0) = C(n, n) = 1

Binomial theorem:
  (a + b)ⁿ = Σₖ₌₀^n C(n,k) · aⁿ⁻ᵏ · bᵏ

Row sums (set a=b=1):
  Σₖ C(n,k) = 2ⁿ   (each row sums to a power of 2)

Fibonacci diagonal sums:
  Σₖ C(n−k, k) = F(n+1)  (shallow diagonal sums = Fibonacci numbers)

Hockey-stick identity:
  Σⱼ₌ᵣ^n C(j, r) = C(n+1, r+1)

Parity pattern (Kummer's theorem):
  C(n, k) is odd ⟺ k AND n = k  in binary (bitwise AND)
  Equivalently: no carrying occurs in base-2 addition of k and (n−k)
  → Odd entries form Sierpiński triangle fractal

Central binomial coefficients:
  C(2n, n) = (2n)! / (n!)²  ≈ 4ⁿ / √(πn)  (Stirling approximation)
  Appears in: expected return time of symmetric random walk,
              Catalan numbers C_n = C(2n,n)/(n+1), lattice path counting

The Pascal's Triangle simulation renders up to 20 rows with three colouring modes: logarithmic scale (relative magnitudes), parity mode (which reveals Sierpiński's triangle), and modular colouring (which selects a prime modulus and colours by residue — producing different p-adic fractals for each prime). The Fibonacci and triangular-number diagonals can be separately highlighted. Hover on any cell to see the exact value and the formula C(n,k).

All Seven Simulations at a Glance

The Unity of Discrete Mathematics

These seven simulations might look like separate topics, but they are deeply connected through the lens of discrete mathematics. Sorting algorithms are specialisations of comparison networks whose optimal depth is related to information-theoretic lower bounds. Graph colouring and N-queens are both constraint satisfaction problems, sitting in the same complexity class (NP-complete for decision versions). Taylor series are the continuous analogue of Pascal's triangle row entries: the coefficients of e^x in a Taylor series are 1/n! = 1/C(n,n) · n!/n! — reciprocals of the diagonal entries of Pascal's triangle.

The Collatz problem defies connection to established mathematics — it stands alone as a reminder that simple rules can generate profound complexity that defeats the best tools we have. Numerical integration, meanwhile, might seem mundane, but it is the practical face of infinite series: Simpson's rule is simply evaluating a degree-3 Taylor polynomial at three points and integrating exactly.

Algorithms are the language in which we express computation. These simulations do not just show algorithms running — they reveal the mathematical structures that explain why they work, where they struggle, and what limits no algorithm can circumvent.