Cache-Oblivious Divide-and-Conquer

Why Cache-Oblivious Divide-and-Conquer Matters

A matrix transpose performs almost no arithmetic, yet it can be slow because it moves data in an inconvenient order. Reading each source row consecutively writes a destination column with a large stride. Several writes may touch different cache lines before the program returns to the next word of any one destination line. If those lines have been evicted, the program fetches and writes them again.

Cache-oblivious divide-and-conquer organizes work into recursively smaller regions without consulting the cache capacity or line length. A depth-first traversal completes a small region and its descendants before moving elsewhere. This creates opportunities for locality at several scales, rather than selecting one hardware-specific tile size.

This lesson makes the idea concrete with an out-of-place matrix transpose. It is a representative algorithm, not a universal recipe for all recursive programs. You will implement the same traversal in Python, C++, and Java, trace every copy, and compare it with a row-major baseline using an explicit cache simulator. The matrix values stay the same; the order of accesses changes.

The distinction connects the earlier B/B+ tree lesson's page-transfer accounting with familiar divide-and-conquer reasoning from sorting. A cache line replaces a disk page as the transfer unit. Applications include reorganizing dense arrays, preparing data for numerical transforms, and designing recursive numerical kernels. The sample demonstrates transfer savings under stated assumptions; hardware speed still requires measurement.

Core Concepts and Terminology

The input is a flat row-major array $A$ of shape $r\times c$. The output is a separate flat row-major array $B$ of shape $c\times r$, satisfying $$ B_{ji}=A_{ij},\qquad 0\le i\lt r,\quad 0\le j\lt c. $$ Let $N=rc$. Source offset $ic+j$ becomes destination offset $jr+i$. The leading dimension, or row stride, is $c$ for A and $r$ for B. These strides remain fixed when recursion selects a smaller rectangle: a view does not repack the underlying array.

A recursive call describes a half-open source rectangle $[r_0,r_1)\times[c_0,c_1)$. Its height and width are $h=r_1-r_0$ and $w=c_1-c_0$. The corresponding output rectangle is $[c_0,c_1)\times[r_0,r_1)$. This interval representation handles odd dimensions and makes disjoint coverage easy to prove. The call's contract is to copy exactly this region to its transposed destination.

Two distinct roles for cache parameters. The transpose traversal knows only bounds, fixed strides, and how to load/store an element. Our observer has $\lambda$ words per cache line and $Q$ line frames, with data capacity $M=Q\lambda$ words. Changing $\lambda$ or $Q$ changes simulated traffic, but never changes the recursive decisions, leaf order, or result. The traversal never asks whether a rectangle fits in cache.

The classical ideal-cache model is fully associative and uses optimal offline replacement. Its tall-cache assumption is $M=\Omega(\lambda^2)$. For positive-size rectangular transpose, the classical bound is $\Theta(1+rc/\lambda)$ cache misses with $\Theta(rc)$ work; see Frigo, Leiserson, Prokop, and Ramachandran's Cache-Oblivious Algorithms, Section 3. Recursion alone does not establish such a bound: the layout and working-set analysis matter.

The executable laboratory uses LRU replacement. Every logical word is eight bytes for accounting, so a line transfers $8\lambda$ bytes and fast data capacity is $8Q\lambda$ bytes. These are modeled widths, not Python object sizes. Word address $a$ belongs to line $\lfloor a/\lambda\rfloor$ and slot $a\bmod\lambda$. Any line can occupy any cache frame; the least recently accessed resident line is evicted when all frames are occupied.

A missing load or store fetches one full line: this is write allocation. A store modifies its cached payload and marks the line dirty. Dirty eviction writes a full line to backing memory; clean eviction does not. At the end, a flush writes each remaining dirty line once and clears its dirty bit. We count $$ T=R+W_e+W_f, $$ where $R$ counts line reads, $W_e$ dirty eviction writes, and $W_f$ final flush writes. Hits are word references with no line fill. An eviction and a writeback are different events; a clean victim contributes only an eviction.

A starts at word zero. B starts at line-aligned address $$ p=\lambda\left\lceil\frac N\lambda\right\rceil. $$ The backing array has $2p$ words, including zero padding and initially zero B storage. Source setup, output initialization, and final inspection are outside the counted kernel; every kernel A load and B store goes through the simulated cache. Each comparison starts with a separate cold cache and finishes with a flush. Recursion metadata, instruction traffic, and trace storage are outside this data-cache model.

The diagram follows the first $2\times2$ source region of the sample. Its two source lines and two destination lines fit exactly in four frames. The program reaches this region by repeated splitting; it does not discover the fit with a capacity test.

How Cache-Oblivious Divide-and-Conquer Works, Step by Step

1. Validate the public input. Check dimensions, cache limits, source length, and values. Compute the aligned B base only after numeric limits are known. If either dimension is zero, return the empty transpose without entering recursion.
2. Start at the full region. Call visit(0,r,0,c,1). The final argument records depth, with the root at depth one. No data is moved merely by making the call.
3. Copy a scalar leaf. If $h=w=1$, load A at $r_0c+c_0$, then store the value at absolute B address $p+c_0r+r_0$. Record the source offset for the explanatory trace. The single scalar temporary is outside the modeled cache payload.
4. Otherwise split the longer dimension. If $h\ge w$, use midpoint $m=r_0+\lfloor h/2\rfloor$ and recurse on the top then bottom source rectangle. Ties go to this row split. If $w\gt h$, use $m=c_0+\lfloor w/2\rfloor$ and visit the left then right rectangle. Both children retain the original A and B strides.
5. Complete both children before returning. There is no numerical combine step: children have written disjoint output regions. The parent simply returns after both are finished. A child may reuse lines loaded by its sibling, depending on the cache and geometry.
6. Flush before observing output. A correct value can still reside only in a dirty cache payload. Flush the remaining dirty lines, then read the returned B values directly from backing memory outside the count.

A positive nonleaf always has a longer dimension at least two. Its floor midpoint lies strictly inside its interval, so both child rectangles are positive and strictly smaller. Guarding zero dimensions before the root prevents an empty interval from splitting forever.

For the diagram's $3\times5$ source, the first split is along columns at two: children have shapes $3\times2$ and $3\times3$. The left child's output is B rows zero and one, while the right child's output is B rows two through four. The left child next splits its three rows at one. The actual code uses a $1\times1$ base case throughout.

A fixed small leaf cutoff independent of cache geometry can preserve cache obliviousness while reducing recursive overhead. A cutoff chosen from the target cache's $M$ or $\lambda$ is an explicit cache-aware choice. Neither changes the requirement to prove coverage and preserve the original strides.

flowchart TD
    A["Validate input and align separate output"] --> B{"Either dimension zero?"}
    B -->|Yes| J["Flush and return output"]
    B -->|No| C["Call Visit on full rectangle; wait for it to return"]
    C --> J
    C -.-> D{"Visit: one row and one column?"}
    D -->|Yes| E["Load A element and store transposed B element"]
    D -->|No| F{"Height at least width?"}
    F -->|Yes| G["Split rows; recursively Visit top then bottom"]
    F -->|No| H["Split columns; recursively Visit left then right"]
    G --> I["Return from this Visit call to its caller"]
    H --> I
    E --> I

Worked Example with a Full Trace

Use $$ A=\begin{bmatrix}1&2&3&4\\5&6&7&8\\9&10&11&12\\13&14&15&16\end{bmatrix}, \qquad r=c=4,\quad\lambda=2,\quad Q=4. $$ A occupies addresses $0$ through $15$, lines $0$ through $7$. B starts at $p=16$ and occupies addresses $16$ through $31$, lines $8$ through $15$. A two-word line is sixteen logical bytes; four frames hold sixty-four bytes.

All recursive regions. The root's row split creates top $[0,2)\times[0,4)$ and bottom $[2,4)\times[0,4)$. Each splits its columns at two. Thus the four $2\times2$ regions are visited top-left, top-right, bottom-left, bottom-right. Each splits its two rows into single rows; each resulting $1\times2$ region splits its two columns into leaves. These are all internal calls: one root, two $2\times4$ calls, four $2\times2$ calls, and eight $1\times2$ calls. There are sixteen leaves, thirty-one calls total, and maximum depth five.

The source offsets are visited in this exact order: $$ 0,1,4,5,\;2,3,6,7,\;8,9,12,13,\;10,11,14,15. $$ The diagram records every copied value and destination offset. Each leaf performs its A load before its B store. M/H means that the source load misses and the destination store hits; M/M means both cause line fills. The reads and writes columns are cumulative before final flush.

Cache state after every leaf. The following independent replay gives line IDs from least to most recently accessed. A star marks a dirty B line. Each event lists its line ID; any victim is shown in parentheses, with a star when its eviction writes a dirty payload. Destination line IDs are absolute, including B's base. This makes all thirty-two references and all twelve evictions reconstructible.

Here dirty is the table's spelling of the star convention. For example, the first copy reads source line zero, fills destination line eight, and writes value one into its first slot. The second copy hits source line zero but fills destination line ten for value two. The third copy fills source line two, then hits destination line eight to add value five beside one. The fourth hits both lines and adds six beside two. Four line reads have now served eight word references.

At copy five, reading source line one evicts clean line zero; the store to line twelve evicts dirty line eight, writing the completed output pair $[1,5]$ to backing memory. Later copies follow the same ownership rule. At the end, resident lines are $[5,13,7,15]$ in LRU order; lines thirteen and fifteen are dirty. Final flush writes their pairs $[11,15]$ and $[12,16]$ and changes no recency order.

The complete transposed result is $$ B=\begin{bmatrix}1&5&9&13\\2&6&10&14\\3&7&11&15\\4&8&12&16\end{bmatrix}. $$ The recursive run has $32$ references, $16$ hits, $16$ reads, $6$ dirty eviction writes, and $2$ flush writes: $T=24$ transfers, or $384$ modeled bytes. The baseline reads A in flat order $0,1,\ldots,15$ and obtains exactly the same B. It has $8$ hits, $24$ reads, $13$ dirty eviction writes, and $3$ flush writes: $T=40$, or $640$ modeled bytes. Its twenty evictions exceed the recursive run's twelve. Both peak at four resident lines.

The saving is $16$ line transfers, or $40\%$ of the baseline's traffic. It is not a measured runtime improvement. Both versions still copy sixteen values and issue thirty-two modeled word references.

Implementation in Python, C++, and Java

Each panel is a complete standalone program with the same sample, traversal, write-back cache, and standard output. Run Python normally; compile C++ as C++17; save the Java panel as Main.java and run Main. The programs print separate line reads, dirty eviction writes, final flush writes, and their transfer sum. They also expose every recursive source offset.

Function contract. Supply A as a flat integer list/vector/array of exactly $rc$ elements. Dimensions must satisfy $0\le r,c\le64$; values lie in $[-10^6,10^6]$; $1\le\lambda\le64$ and $1\le Q\le128$. Python's final recursive argument defaults to True; C++ has the same default; Java requires the boolean argument explicitly. False selects the row-major baseline. The source is unchanged, and a separate result contains B and all counters.

Invalid dimension limits produce DIMENSION, invalid cache limits PARAMETER, wrong flat length SHAPE, and out-of-range values VALUE. Checks occur in that order before dimension multiplication or address arithmetic. Python first rejects non-exact integer parameter types and a non-exact boolean mode as PARAMETER; it rejects booleans, floats, strings, and other non-integer element types as VALUE. C++ uses typed int parameters, bool mode, and long long words; Java uses int, boolean, and long. The contract expects a non-null array/container; null or unrelated container objects are outside it.

Bounds keep all address/count arithmetic small, and values are copied without arithmetic. No matrix-size-dependent allocation occurs inside the recursive traversal. Python uses integer lists, C++ vectors, and Java primitive arrays; each cache line holds an actual copy of its backing payload. A store updates that payload, so forgetting dirty writeback would corrupt the observed output rather than merely change a counter.

Read the kernel separately from the observer. The visit method/function only reads interval lengths and divides them. The cache object receives $\lambda,Q$, and the initial storage setup uses $\lambda$ for alignment. Neither value participates in splitting or the scalar stopping rule. The source-order trace and call statistics remain identical for every allowed cache geometry. They are teaching instrumentation, not a requirement of a production transpose.

Zero dimensions produce an empty B with zero accesses, calls, evictions, writes, and peak resident lines. A singleton has one leaf call at depth one, one load, and one store; with only one frame, the source line is evicted before the distinct destination line can be fetched. Stores still allocate a line even though B was initially zero: that policy is part of the experiment.

from dataclasses import dataclass
from collections import OrderedDict


@dataclass
class Line:
    values: list
    dirty: bool = False


@dataclass
class Result:
    values: list
    accesses: int
    hits: int
    reads: int
    eviction_writes: int
    flush_writes: int
    evictions: int
    peak_lines: int
    calls: int
    max_depth: int
    order: list


class Cache:
    def __init__(self, memory, line_words, frames):
        self.memory, self.line_words, self.frames = memory, line_words, frames
        self.lines = OrderedDict()
        self.accesses = self.hits = self.reads = 0
        self.eviction_writes = self.flush_writes = self.evictions = self.peak_lines = 0

    def writeback(self, line_id, line):
        start = line_id * self.line_words
        self.memory[start:start + self.line_words] = line.values
        line.dirty = False

    def fetch(self, address):
        assert 0 <= address < len(self.memory)
        self.accesses += 1
        line_id = address // self.line_words
        if line_id in self.lines:
            self.hits += 1
            self.lines.move_to_end(line_id)
            return self.lines[line_id]
        self.reads += 1
        if len(self.lines) == self.frames:
            old_id, old = self.lines.popitem(last=False)
            self.evictions += 1
            if old.dirty:
                self.writeback(old_id, old)
                self.eviction_writes += 1
        start = line_id * self.line_words
        line = Line(self.memory[start:start + self.line_words])
        self.lines[line_id] = line
        self.peak_lines = max(self.peak_lines, len(self.lines))
        return line

    def load(self, address):
        return self.fetch(address).values[address % self.line_words]

    def store(self, address, value):
        line = self.fetch(address)
        line.values[address % self.line_words] = value
        line.dirty = True

    def flush(self):
        for line_id, line in self.lines.items():
            if line.dirty:
                self.writeback(line_id, line)
                self.flush_writes += 1

def transpose(a, rows, cols, line_words, frames, recursive=True):
    if (any(type(x) is not int for x in (rows, cols, line_words, frames))
            or type(recursive) is not bool):
        raise ValueError("PARAMETER")
    if not (0 <= rows <= 64 and 0 <= cols <= 64):
        raise ValueError("DIMENSION")
    if not (1 <= line_words <= 64 and 1 <= frames <= 128):
        raise ValueError("PARAMETER")
    if len(a) != rows * cols:
        raise ValueError("SHAPE")
    if any(type(x) is not int or not -1000000 <= x <= 1000000 for x in a):
        raise ValueError("VALUE")
    n = rows * cols
    output_base = (n + line_words - 1) // line_words * line_words
    memory = a.copy() + [0] * (2 * output_base - n)
    cache = Cache(memory, line_words, frames)
    order = []
    calls = max_depth = 0

    def copy_cell(i, j):
        source = i * cols + j
        value = cache.load(source)
        cache.store(output_base + j * rows + i, value)
        order.append(source)

    def visit(r0, r1, c0, c1, depth):
        nonlocal calls, max_depth
        calls += 1
        max_depth = max(max_depth, depth)
        height, width = r1 - r0, c1 - c0
        if height == 1 and width == 1:
            copy_cell(r0, c0)
        elif height >= width:
            mid = r0 + height // 2
            visit(r0, mid, c0, c1, depth + 1)
            visit(mid, r1, c0, c1, depth + 1)
        else:
            mid = c0 + width // 2
            visit(r0, r1, c0, mid, depth + 1)
            visit(r0, r1, mid, c1, depth + 1)
    if rows and cols:
        if recursive:
            visit(0, rows, 0, cols, 1)
        else:
            for i in range(rows):
                for j in range(cols):
                    copy_cell(i, j)
    cache.flush()
    return Result(memory[output_base:output_base + n], cache.accesses, cache.hits,
                  cache.reads, cache.eviction_writes, cache.flush_writes,
                  cache.evictions, cache.peak_lines, calls, max_depth, order)


def show(label, r):
    print(f"{label} accesses={r.accesses} hits={r.hits} reads={r.reads} "
          f"eviction_writes={r.eviction_writes} flush_writes={r.flush_writes} "
          f"transfers={r.reads + r.eviction_writes + r.flush_writes} "
          f"evictions={r.evictions} peak_lines={r.peak_lines} "
          f"calls={r.calls} max_depth={r.max_depth}")


def main():
    a = list(range(1, 17))
    plain = transpose(a, 4, 4, 2, 4, False)
    rec = transpose(a, 4, 4, 2, 4, True)
    assert plain.values == rec.values
    print("shape=4x4 -> 4x4 line_words=2 cache_lines=4")
    print("line_bytes=16 cache_bytes=64 equal=true")
    print("B=[" + ",".join(map(str, rec.values)) + "]")
    print("recursive_source_order=[" + ",".join(map(str, rec.order)) + "]")
    show("row_major", plain)
    show("recursive", rec)


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cassert>
#include <iostream>
#include <list>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <vector>
using namespace std;
using Word = long long;

struct Result {
    vector<Word> values;
    int accesses, hits, reads, eviction_writes, flush_writes, evictions, peak_lines;
    int calls, max_depth;
    vector<int> order;
};
struct Line {
    vector<Word> values;
    bool dirty;
    list<int>::iterator position;
};
struct Cache {
    vector<Word>& memory;
    int line_words, frames;
    list<int> lru;
    unordered_map<int,Line> lines;
    int accesses = 0, hits = 0, reads = 0;
    int eviction_writes = 0, flush_writes = 0, evictions = 0, peak_lines = 0;
    Cache(vector<Word>& m, int l, int q) : memory(m), line_words(l), frames(q) {}
    Cache(const Cache&) = delete;
    Cache& operator=(const Cache&) = delete;
    void writeback(int id, Line& line) {
        copy(line.values.begin(), line.values.end(), memory.begin() + id * line_words);
        line.dirty = false;
    }
    Line& fetch(int address) {
        assert(address >= 0 && address < static_cast<int>(memory.size()));
        ++accesses;
        int id = address / line_words;
        auto found = lines.find(id);
        if (found != lines.end()) {
            ++hits; lru.splice(lru.end(), lru, found->second.position); return found->second;
        }
        ++reads;
        if (static_cast<int>(lines.size()) == frames) {
            int old_id = lru.front(); auto old = lines.find(old_id);
            ++evictions;
            if (old->second.dirty) { writeback(old_id, old->second); ++eviction_writes; }
            lines.erase(old); lru.pop_front();
        }
        int start = id * line_words;
        vector<Word> values(memory.begin() + start, memory.begin() + start + line_words);
        lru.push_back(id);
        auto inserted = lines.emplace(id, Line{move(values), false, prev(lru.end())});
        peak_lines = max(peak_lines, static_cast<int>(lines.size()));
        return inserted.first->second;
    }
    Word load(int address) { return fetch(address).values[address % line_words]; }
    void store(int address, Word value) {
        Line& line = fetch(address); line.values[address % line_words] = value; line.dirty = true;
    }
    void flush() {
        for (auto& entry : lines)
            if (entry.second.dirty) { writeback(entry.first, entry.second); ++flush_writes; }
    }
};

struct Traversal {
    Cache& cache;
    int rows, cols, output_base, calls = 0, max_depth = 0;
    vector<int> order;
    Traversal(Cache& c, int r, int s, int base) : cache(c), rows(r), cols(s), output_base(base) {}
    void copy_cell(int i, int j) {
        int source = i * cols + j;
        Word value = cache.load(source);
        cache.store(output_base + j * rows + i, value);
        order.push_back(source);
    }
    void visit(int r0, int r1, int c0, int c1, int depth) {
        ++calls; max_depth = max(max_depth, depth);
        int height = r1 - r0, width = c1 - c0;
        if (height == 1 && width == 1) copy_cell(r0, c0);
        else if (height >= width) {
            int mid = r0 + height / 2;
            visit(r0, mid, c0, c1, depth + 1);
            visit(mid, r1, c0, c1, depth + 1);
        } else {
            int mid = c0 + width / 2;
            visit(r0, r1, c0, mid, depth + 1);
            visit(r0, r1, mid, c1, depth + 1);
        }
    }
};
Result transpose(const vector<Word>& a, int rows, int cols, int line_words,
                 int frames, bool recursive = true) {
    if (rows < 0 || rows > 64 || cols < 0 || cols > 64)
        throw invalid_argument("DIMENSION");
    if (line_words < 1 || line_words > 64 || frames < 1 || frames > 128)
        throw invalid_argument("PARAMETER");
    int n = rows * cols;
    if (a.size() != static_cast<size_t>(n)) throw invalid_argument("SHAPE");
    for (Word x : a) if (x < -1000000 || x > 1000000) throw invalid_argument("VALUE");
    int base = (n + line_words - 1) / line_words * line_words;
    vector<Word> memory(2 * base, 0);
    copy(a.begin(), a.end(), memory.begin());
    Cache cache(memory, line_words, frames);
    Traversal traversal(cache, rows, cols, base);
    if (rows && cols) {
        if (recursive) traversal.visit(0, rows, 0, cols, 1);
        else for (int i = 0; i < rows; ++i)
            for (int j = 0; j < cols; ++j) traversal.copy_cell(i, j);
    }
    cache.flush();
    return {vector<Word>(memory.begin()+base, memory.begin()+base+n),
            cache.accesses, cache.hits, cache.reads, cache.eviction_writes,
            cache.flush_writes, cache.evictions, cache.peak_lines,
            traversal.calls, traversal.max_depth, move(traversal.order)};
}
void show(const string& label, const Result& r) {
    cout << label << " accesses=" << r.accesses << " hits=" << r.hits << " reads=" << r.reads
         << " eviction_writes=" << r.eviction_writes << " flush_writes=" << r.flush_writes
         << " transfers=" << r.reads+r.eviction_writes+r.flush_writes
         << " evictions=" << r.evictions << " peak_lines=" << r.peak_lines
         << " calls=" << r.calls << " max_depth=" << r.max_depth << '\n';
}
int main() {
    vector<Word> a; for (int i = 1; i <= 16; ++i) a.push_back(i);
    Result plain = transpose(a,4,4,2,4,false), rec = transpose(a,4,4,2,4,true);
    assert(plain.values == rec.values);
    cout << "shape=4x4 -> 4x4 line_words=2 cache_lines=4\n"
         << "line_bytes=16 cache_bytes=64 equal=true\nB=[";
    for (size_t i=0;i<rec.values.size();++i) { if(i)cout<<','; cout<<rec.values[i]; }
    cout << "]\nrecursive_source_order=[";
    for (size_t i=0;i<rec.order.size();++i) { if(i)cout<<','; cout<<rec.order[i]; }
    cout << "]\n"; show("row_major",plain); show("recursive",rec);
    return 0;
}

import java.util.Arrays;
import java.util.LinkedHashMap;
import java.util.Map;

public class Main {
    static final class Result {
        final long[] values;
        final int accesses, hits, reads, evictionWrites, flushWrites, evictions, peakLines;
        final int calls, maxDepth;
        final int[] order;
        Result(long[] v, Cache c, int calls, int maxDepth, int[] order) {
            this.calls = calls; this.maxDepth = maxDepth; this.order = order;
            values = v; accesses = c.accesses; hits = c.hits; reads = c.reads;
            evictionWrites = c.evictionWrites; flushWrites = c.flushWrites;
            evictions = c.evictions; peakLines = c.peakLines;
        }
    }
    static final class Line {
        final long[] values;
        boolean dirty;
        Line(long[] values) { this.values = values; }
    }
    static final class Cache {
        final long[] memory;
        final int lineWords, frames;
        final LinkedHashMap<Integer,Line> lines = new LinkedHashMap<>(16,0.75f,true);
        int accesses, hits, reads, evictionWrites, flushWrites, evictions, peakLines;
        Cache(long[] memory, int lineWords, int frames) {
            this.memory = memory; this.lineWords = lineWords; this.frames = frames;
        }
        void writeback(int id, Line line) {
            System.arraycopy(line.values,0,memory,id*lineWords,lineWords); line.dirty = false;
        }
        Line fetch(int address) {
            if (address < 0 || address >= memory.length) throw new AssertionError("address");
            ++accesses;
            int id = address / lineWords; Line line = lines.get(id);
            if (line != null) { ++hits; return line; }
            ++reads;
            if (lines.size() == frames) {
                Map.Entry<Integer,Line> old = lines.entrySet().iterator().next();
                ++evictions;
                if (old.getValue().dirty) { writeback(old.getKey(),old.getValue()); ++evictionWrites; }
                lines.remove(old.getKey());
            }
            int start = id*lineWords;
            line = new Line(Arrays.copyOfRange(memory,start,start+lineWords));
            lines.put(id,line); peakLines = Math.max(peakLines,lines.size()); return line;
        }
        long load(int address) { return fetch(address).values[address%lineWords]; }
        void store(int address, long value) {
            Line line = fetch(address); line.values[address%lineWords] = value; line.dirty = true;
        }
        void flush() {
            for (Map.Entry<Integer,Line> e : lines.entrySet())
                if (e.getValue().dirty) { writeback(e.getKey(),e.getValue()); ++flushWrites; }
        }
    }
    static int padded(int n, int lineWords) { return (n+lineWords-1)/lineWords*lineWords; }
    static final class Traversal {
        final Cache cache;
        final int rows, cols, outputBase;
        final int[] order;
        int copied, calls, maxDepth;
        Traversal(Cache cache, int rows, int cols, int base) {
            this.cache=cache; this.rows=rows; this.cols=cols; outputBase=base;
            order=new int[rows*cols];
        }
        void copyCell(int i, int j) {
            int source=i*cols+j;
            long value=cache.load(source);
            cache.store(outputBase+j*rows+i,value);
            order[copied++]=source;
        }
        void visit(int r0, int r1, int c0, int c1, int depth) {
            ++calls; maxDepth=Math.max(maxDepth,depth);
            int height=r1-r0,width=c1-c0;
            if(height==1 && width==1) copyCell(r0,c0);
            else if(height>=width) {
                int mid=r0+height/2;
                visit(r0,mid,c0,c1,depth+1); visit(mid,r1,c0,c1,depth+1);
            } else {
                int mid=c0+width/2;
                visit(r0,r1,c0,mid,depth+1); visit(r0,r1,mid,c1,depth+1);
            }
        }
    }
    static Result transpose(long[] a, int rows, int cols, int lineWords,
                            int frames, boolean recursive) {
        if(rows<0 || rows>64 || cols<0 || cols>64) throw new IllegalArgumentException("DIMENSION");
        if(lineWords<1 || lineWords>64 || frames<1 || frames>128) throw new IllegalArgumentException("PARAMETER");
        int n=rows*cols;
        if(a.length!=n) throw new IllegalArgumentException("SHAPE");
        for(long x:a) if(x < -1000000 || x > 1000000) throw new IllegalArgumentException("VALUE");
        int base=padded(n,lineWords);
        long[] memory=new long[2*base]; System.arraycopy(a,0,memory,0,n);
        Cache cache=new Cache(memory,lineWords,frames);
        Traversal traversal=new Traversal(cache,rows,cols,base);
        if(rows!=0 && cols!=0) {
            if(recursive) traversal.visit(0,rows,0,cols,1);
            else for(int i=0;i<rows;++i) for(int j=0;j<cols;++j) traversal.copyCell(i,j);
        }
        cache.flush();
        return new Result(Arrays.copyOfRange(memory,base,base+n),cache,
                          traversal.calls,traversal.maxDepth,traversal.order);
    }
    static void show(String label, Result r) {
        System.out.println(label+" accesses="+r.accesses+" hits="+r.hits+" reads="+r.reads
            +" eviction_writes="+r.evictionWrites+" flush_writes="+r.flushWrites
            +" transfers="+(r.reads+r.evictionWrites+r.flushWrites)
            +" evictions="+r.evictions+" peak_lines="+r.peakLines
            +" calls="+r.calls+" max_depth="+r.maxDepth);
    }
    static String compact(long[] a) { return Arrays.toString(a).replace(" ",""); }
    static String compact(int[] a) { return Arrays.toString(a).replace(" ",""); }
    public static void main(String[] args) {
        long[] a=new long[16]; for(int i=0;i<16;++i)a[i]=i+1;
        Result plain=transpose(a,4,4,2,4,false),rec=transpose(a,4,4,2,4,true);
        if(!Arrays.equals(plain.values,rec.values)) throw new AssertionError("transpose");
        System.out.println("shape=4x4 -> 4x4 line_words=2 cache_lines=4");
        System.out.println("line_bytes=16 cache_bytes=64 equal=true");
        System.out.println("B="+compact(rec.values));
        System.out.println("recursive_source_order="+compact(rec.order));
        show("row_major",plain); show("recursive",rec);
    }
}

Expected standard output, identical in all three languages:

• shape=4x4 -> 4x4 line_words=2 cache_lines=4
• line_bytes=16 cache_bytes=64 equal=true
• B=[1,5,9,13,2,6,10,14,3,7,11,15,4,8,12,16]
• recursive_source_order=[0,1,4,5,2,3,6,7,8,9,12,13,10,11,14,15]
• row_major accesses=32 hits=8 reads=24 eviction_writes=13 flush_writes=3 transfers=40 evictions=20 peak_lines=4 calls=0 max_depth=0
• recursive accesses=32 hits=16 reads=16 eviction_writes=6 flush_writes=2 transfers=24 evictions=12 peak_lines=4 calls=31 max_depth=5

Complexity, Correctness, and Invariants

Work, best/average/worst cases. For positive $N=rc$, every leaf copies one element and every internal call has exactly two children. There are $N$ leaves and $N-1$ internal calls, hence $$ \text{calls}=2N-1,\qquad U=2N,\qquad H+R=2N. $$ Here $U$ counts word references and $H$ hits. Work is $\Theta(N)$ in best, average, and worst cases for fixed dimensions; values do not affect traversal. There is no probabilistic average-case assumption for cache performance. The row-major baseline has the same $2N$ references, no recursive calls, and constant kernel control storage.

Each row split halves a positive row interval, and each column split halves a positive column interval. Root depth is one, so $$ \text{maximum depth}\le1+\lceil\log_2 r\rceil+\lceil\log_2 c\rceil. $$ The recursive control stack is $O(\log r+\log c)$ for positive dimensions. Calls execute sequentially: sibling subtrees do not both remain expanded on the stack. For an empty shape, kernel work/calls/traffic are zero, while validation/setup still has a fixed cost.

Exact traffic bounds for this simulator. For $N\gt0$, let $D=\lceil N/\lambda\rceil$. Distinct aligned A and B regions each touch $D$ lines. Every A line must be read, every B line must be allocated at its first store, and every B line must eventually be written back. Thus $$ R\ge2D,\qquad W_e+W_f\ge D,\qquad T\ge3D. $$ A writeback needs at least one store since that line's previous writeback or initial load. There are $N$ stores and $2N$ accesses, so $$ R\le2N,\qquad W_e+W_f\le N,\qquad 3D\le T\le3N. $$ These bounds include final flush and partial padded lines. A sufficient capacity for both orders to attain $3D$ is $Q\ge2D$, retaining all touched lines, but it is not necessary. In the sample, $D=8$ and the lower bound is $24$ transfers; recursion attains it with four frames even though full residency would require sixteen. With $\lambda=1$, the bounds meet at $3N$, and this model offers neither traversal a line-traffic advantage.

The sample's capacity sweep, with $r=c=4$ and $\lambda=2$, is:

This shows both improvement and ties. The exact finite traffic depends on shape, alignment, replacement, write policy, and capacity; a single asymptotic expression is not a counter prediction.

Why recursive locality can be asymptotically optimal. For near-square regions, conceptually stop the analysis at side length $s$ on the order of $\sqrt M$, with a sufficiently small constant for headroom. Source plus destination occupy $O(s^2/\lambda+s)$ lines, including row fragments. Tall cache gives $s=\Omega(\lambda)$, making this $O(s^2/\lambda)$ and allowing the footprint to fit. Descendants continue executing, but finish inside that footprint. Summing over $O(N/s^2)$ such groups gives $O(N/\lambda)$ movement when both dimensions are sufficiently large; add the constant for a small whole input. Thin rectangles need a separate boundary argument, not an area-only fit claim. Prokop's thesis gives the full rectangular analysis.

The ideal-cache theorem counts misses under optimal replacement. Our finite numbers use LRU and explicitly include writes. LRU can be related to optimal replacement with constant-factor capacity comparisons, and regularity conditions support asymptotic transfer of suitable bounds; this does not make their same-capacity miss counts equal. In particular, arbitrary simulator settings need not meet the theorem's tall-cache assumptions. A write-allocating write-back execution with $O(1+N/\lambda)$ fills also has only that order of dirty writebacks, because each loaded line can be written back at most once before its next load, including final flush.

Space and simulation cost. Out-of-place output requires $\Theta(N)$ words. The demonstration additionally retains the input copy, padded backing memory, returned result, and $N$ trace indices: $O(N+\lambda)$ host storage excluding cache payload. Resident cache payload is at most $Q\lambda$ words, with $O(Q)$ replacement metadata, plus the logarithmic control stack. Trace arrays and metadata are outside the modeled data-cache capacity. With expected constant-time hash lookup, the software simulator costs $O(N+\lambda T)$ plus setup, because it copies a whole payload on each read/write transfer. Its timing is not the hardware timing of an uninstrumented transpose.

Correctness by region induction. A call on $[r_0,r_1)\times[c_0,c_1)$ must copy exactly the cells in that source rectangle to $[c_0,c_1)\times[r_0,r_1)$ in B. At a leaf, the global offset formulas write exactly the required value. A row split partitions the source into disjoint top/bottom rectangles and partitions the output into corresponding disjoint column regions. A column split partitions source columns and output rows. By induction, both children fulfill their contracts, and their union fulfills the parent's contract without overwriting another child's cell. The decreasing positive area proves termination. At the root, every output cell is written once.

Cache-coherence invariant. For any logical word, its newest value is in its dirty resident payload if present, otherwise in backing memory. Fetch preserves the current value; store modifies the authoritative resident payload; dirty eviction writes it before removal. Flush therefore materializes every B value before result inspection. Disjoint A/B regions ensure no output write destroys an unread source value. This memory invariant and the region proof establish correct values independently of the cache's hit rate.

Variants, Pitfalls, and Optimizations

Preserve strides when splitting views. A $2\times2$ view inside a $3\times5$ A does not have source stride two. Using its width as the leading dimension reads unrelated cells. Keep original $c$ and $r$ in the global offset formulas; change only interval bounds. The diagram contrasts a correct source value six with a wrong value three caused by this mistake.

Handle odd and empty dimensions. Use half-open intervals and floor midpoints. The longer dimension of a positive nonleaf is at least two, so neither child is empty. Do not recurse at all when the original shape has a zero dimension. Tests like $3\times5$, $1\times64$, $64\times1$, and $0\times64$ expose assumptions hidden by square powers of two.

A recursive split does not copy a submatrix. Python slicing that constructs new matrices, or packing at every recursion level, adds allocation and memory movement beyond this model. Carry bounds and strides or a true noncopying view. An explicit depth-first stack can replace language recursion while preserving visit order; push the second child before the first when using a LIFO stack.

Use a suitable base case. Scalar leaves make the proof and trace simple but can be expensive in practice. A small fixed loop kernel reduces call overhead while keeping the traversal independent of $M,\lambda$. Selecting the cutoff by cache-size tuning makes that part cache aware. Vectorized or register-level transpose kernels can be useful leaves, but their reference trace and measured speed differ from the supplied scalar laboratory.

Out of place and in place are different contracts. The program writes into a separate aligned B region. Reusing the same storage overwrites unread values: for a $2\times2$ matrix, writing $B_{10}=A_{01}$ destroys the old $A_{10}$ before it is read. Square in-place transpose can swap pairs and treat diagonal blocks specially. Rectangular in-place transpose generally needs more careful permutation/cycle handling. Reusing the presented copy operation unchanged is incorrect.

Do not infer cache optimality from recursion alone. Balanced recursion can still access poor layouts or revisit too much data between uses. Long-stride row fragments affect the true line footprint; sufficient area in words does not by itself prove sufficient line capacity. Recursive FFT, matrix multiplication, and static recursive search layouts require their own dependency, layout, and cache analyses. Ordinary two-way mergesort is recursive, but that fact alone does not establish the optimal external-memory sorting bound; cache-oblivious sorting structures such as funnels need additional machinery.

Multilevel opportunity has assumptions. One recursive schedule contains regions at many sizes, which motivates reuse at more than one memory level. The classic multilevel results impose hierarchy and replacement conditions. Real caches have limited associativity, alignment effects, prefetch, TLB behavior, write buffers, and coherence traffic. Our one-level fully associative model omits all of those, including stack and instruction references. The theorem is not a universal guarantee for every machine.

Compare equivalent boundaries. Both runs begin cold, use the same logical layout, perform one load then one store per value, and flush at the same point. Omitting flush hides dirty output traffic. Changing to non-temporal stores, assuming an entire overwritten line needs no allocation, or counting initialization in only one run changes the experiment. Such policies can be valuable, but must be described and measured separately.

Measure useful outcomes. The $40\to24$ transfer change is a verified model result. A small scalar recursive Python implementation may still run slower than a simple loop because of call and object overhead. For production numerical work, compare uninstrumented kernels or established array-library operations with the actual layout, repeated measurements, and hardware counters where available. Check output first; cache efficiency cannot repair a wrong permutation.

Practice Problems and Real-World Use

Build an independent value oracle by enumerating B in row-major order: for each output row $j$, append $A_{0j},A_{1j},\ldots,A_{r-1,j}$. Do not use the implementation's recursive partition logic to compute expected values. Separately enumerate source offsets with an explicit stack and replay their load/store addresses using a timestamp-based recency oracle. This checks permutation correctness, split order, and cache accounting independently.

Useful properties include: the source stays unchanged; B has $rc$ values; the source-order trace is a permutation of $0,\ldots,N-1$ for positive N; recursive calls equal $2N-1$; every access is a hit or a line read; peak resident lines never exceed $Q$; and final backing output equals the direct transpose. Transposing twice with reversed shapes restores the original values, though its address trace and cache counters need not be the same. Changing values or cache geometry must not change this algorithm's recursive source order.

The supplied programs passed $3669$ shared inputs per language: $3654$ accepted cases and $15$ rejected cases. Tests compare every output value, leaf offset, call/depth statistic, and cache counter against independent flat-transpose, explicit-stack, and timestamp-cache oracles. They include exhaustive shapes from $0$ through $8$ in each dimension across several cache settings, randomized rectangles, odd and skinny shapes, $64\times64$ bounds, extreme allowed values, invalid numeric limits, bad lengths, and out-of-range values. Python additionally passed $24$ exact-type rejection checks, $3200$ generic load/store checks with backing-memory verification and idempotent flush, and $400$ comparisons of traversal independence across cache geometries. C++ ran with address and undefined-behavior sanitizers. All three standalone samples print byte-identical output.

In systems work, a materialized transpose can make a later traversal contiguous, but its cost must be amortized by that later work. A strided view alone can avoid the copy when consumers support it; a consumer requiring contiguous rows may justify materialization. Recursive tiling also informs image layout conversion and numerical kernels, but boundaries, dependencies, and the actual storage format still determine legality and benefit.

For deeper study, the original paper provides separate treatments of transpose, multiplication, FFT, and cache-oblivious sorting, while Prokop's thesis develops the model and analyses at greater length. Use the worked laboratory to understand exact traffic, then read those proofs to distinguish a portable traversal from a theorem about its memory behavior.