Building a CPU Transformer Library

Why C is absolutely viable

C is the most emotionally honest language for this project. It gives you explicit file I/O, explicit memory ownership, explicit flat buffers, and easy access to SIMD intrinsics once the scalar version works. More importantly, the existence proof is overwhelming: llama.cpp and similar projects demonstrate that serious CPU transformer inference in C or C++ is not just possible but highly competitive.

The recommended file structure is small and unsurprising:

gguf.h / gguf.c        // file parsing and metadata
tokenizer.h / tokenizer.c
tensor.h / tensor.c    // matmul, rmsnorm, softmax, rope
model.h / model.c      // layer structs, KV cache, forward pass
main.c                 // CLI and generation loop

Start by defining plain structs with no hidden ownership rules. Resist the temptation to generalize early. A transformer checkpoint is regular enough that fixed-shape assumptions and a few dedicated structs will carry you surprisingly far.

typedef struct {
    void   *data;
    int     ndim;
    int64_t shape[4];
    int     gguf_type;
} Tensor;

typedef struct {
    Tensor wq, wk, wv, wo;
    Tensor w_gate, w_up, w_down;
    Tensor attn_norm, ffn_norm;
} Layer;

typedef struct {
    float *keys;
    float *values;
    int    seq_len;
    int    head_dim;
    int    n_kv_heads;
} KVCache;

typedef struct {
    int32_t vocab_size;
    int32_t n_layers;
    int32_t d_model;
    int32_t n_heads;
    int32_t n_kv_heads;
    int32_t head_dim;
    int32_t max_seq_len;
    float   rms_eps;
    Tensor  tok_embeddings;
    Tensor  final_norm;
    Tensor  lm_head;
    Layer  *layers;
    KVCache *cache;
} Model;

The GGUF loader is mostly a disciplined binary reader. Read the magic, verify the version, pull the metadata count and tensor count, then walk the metadata entries and tensor descriptors into temporary structs. Once you know all tensor names, shapes, types, and offsets, you can populate your Model object by name. Many implementations keep a name-to-descriptor table during loading, then resolve required tensors such as blk.0.attn_q.weight or the model-family equivalent.

static uint32_t read_u32(FILE *f) {
    uint32_t v;
    if (fread(&v, sizeof(v), 1, f) != 1) exit(1);
    return v;
}

static uint64_t read_u64(FILE *f) {
    uint64_t v;
    if (fread(&v, sizeof(v), 1, f) != 1) exit(1);
    return v;
}

static void read_exact(FILE *f, void *dst, size_t n) {
    if (fread(dst, 1, n, f) != n) exit(1);
}

int gguf_open(const char *path, GGUFFile *out) {
    FILE *f = fopen(path, "rb");
    if (!f) return -1;

    uint32_t magic = read_u32(f);
    uint32_t version = read_u32(f);
    if (magic != 0x46554747u) return -2; /* GGUF */
    if (version < 2 || version > 3) return -3;

    out->n_tensors = read_u64(f);
    out->n_kv      = read_u64(f);
    out->fp = f;
    return 0;
}

After the header comes the real parsing work: strings, metadata variants, tensor names, dimensions, and aligned offsets. The clean design is to parse descriptors first and only then decide how to expose the weight bytes. For a first implementation you can malloc a buffer per tensor and fread the raw bytes into it. For large models that is not ideal, but it is easy to reason about. Once correct, you can switch to memory mapping or direct file-backed views.

The biggest loader bug class is metadata mismatch. Validate every tensor shape against the config-derived expectation before you launch any compute. If the config says d_model = 2048 and your q-projection tensor claims a second dimension of 2304, stop immediately. C gives you no runtime safety net here beyond the checks you write yourself.

For tensor math, begin with a reference matmul so plain it hurts. A triple loop is enough for correctness validation. If your runtime only supports single-token decode at first, a matrix-vector multiply is often the actual hot path.

void matmul_f32(float *out, const float *w, const float *x,
                int rows, int cols) {
    for (int r = 0; r < rows; ++r) {
        float acc = 0.0f;
        const float *wr = w + (size_t)r * cols;
        for (int c = 0; c < cols; ++c) {
            acc += wr[c] * x[c];
        }
        out[r] = acc;
    }
}

That version is enough to verify embeddings, norms, projections, and the final output head. Once it is correct, tile it for cache and vectorize it. A typical next step is AVX2 or NEON intrinsics that load 8 or 16 floats at a time, multiply, and accumulate into vector registers before a horizontal reduction. The conceptual math does not change; only the way you feed the CPU pipelines changes.

Quantized matmul is where GGUF-specific engineering enters. Instead of expanding a Q4 or Q6 block into a long temporary FP32 array, dequantize block-by-block inside the dot product. The kernel pattern is: read the block’s scale metadata, unpack the small integers, reconstruct approximate float values for a short span such as 32 weights, and accumulate against the matching activation fragment. This keeps memory traffic lower and matches the storage format’s design.

float dot_q4_0_block(const Q4_0_Block *blk, const float *x) {
    float acc = 0.0f;
    float d = blk->d;
    for (int i = 0; i < 16; ++i) {
        uint8_t q = blk->qs[i];
        int lo = (q & 0x0F) - 8;
        int hi = (q >> 4) - 8;
        acc += (lo * d) * x[i * 2 + 0];
        acc += (hi * d) * x[i * 2 + 1];
    }
    return acc;
}

The exact block structs differ by quantization type, but the idea repeats. A production kernel would vectorize the unpack and accumulation path, perhaps processing several blocks per iteration. The reference version, however, should be scalar and readable. Debuggability beats cleverness during the first pass.

RMSNorm and softmax are small but essential. Keep their reference forms around even after you optimize other parts, because they make excellent correctness anchors.

void rmsnorm(float *y, const float *x, const float *gamma,
             int n, float eps) {
    float sumsq = 0.0f;
    for (int i = 0; i < n; ++i) {
        sumsq += x[i] * x[i];
    }
    float inv = 1.0f / sqrtf((sumsq / (float)n) + eps);
    for (int i = 0; i < n; ++i) {
        y[i] = x[i] * inv * gamma[i];
    }
}

void softmax(float *x, int n) {
    float maxv = x[0];
    for (int i = 1; i < n; ++i) if (x[i] > maxv) maxv = x[i];
    float sum = 0.0f;
    for (int i = 0; i < n; ++i) {
        x[i] = expf(x[i] - maxv);
        sum += x[i];
    }
    for (int i = 0; i < n; ++i) x[i] /= sum;
}

Memory management deserves explicit planning because a transformer has both long-lived and short-lived state. Weights live for the life of the process. The KV cache grows with sequence length. Activations, scratch buffers, and temporary projection outputs are ephemeral and benefit from an arena allocator or a reusable workspace. A simple approach is one arena per inference context that is reset between decoding steps but sized large enough for the largest intermediate tensors used during one forward pass.

The final application loop in main.c can be pleasantly small once the pieces exist: load the model, tokenize the prompt, prefill the cache over the prompt tokens, then repeatedly run one decode step, sample a token, print it, and stop at EOS. If you want streaming output, decode each sampled token immediately with the tokenizer and printf it as it arrives.

gcc -O3 -march=native -o inference main.c tensor.c gguf.c tokenizer.c -lm

That build line is intentionally minimal. You can add platform-specific SIMD flags, profiling hooks, and threading later. The important thing is that C lets you start from zero dependencies and grow only what the profiler justifies. Once your engine answers "What is the capital of France?" with a plausible token stream, you can iteratively harden it into something much closer to production quality.

Why C# is a serious systems option here

C# on .NET 8+ is far more suitable for this project than old stereotypes suggest. You get strong binary I/O primitives, Span<T> and Memory<T> for zero-copy slicing, MemoryMappedFile for large checkpoints, System.Numerics.Vector<float> for portable SIMD, and hardware intrinsics when you want to get closer to the metal. You also get a pleasant type system for representing model metadata without fighting raw pointers everywhere.

A clean file layout looks like this:

GgufReader.cs
Tokenizer.cs
Tensor.cs
Model.cs
Sampling.cs
Program.cs

The central type can stay very small. A tensor is just a flat memory region plus a shape. Whether the backing store comes from a managed array, a memory-mapped view, or a custom quantized buffer is an implementation detail that higher layers should not need to care about immediately.

public enum TensorKind
{
    F32,
    F16,
    Q4_0,
    Q4_K_M,
    Q6_K,
    Q8_0
}

public sealed record Tensor(Memory<float> Data, int[] Shape);

public sealed record QuantizedTensor(
    Memory<byte> Raw,
    int[] Shape,
    TensorKind Kind,
    int BlockSize);

The GGUF reader can use BinaryReader for the structured header and then switch to memory mapping for the heavy tensor region. That split keeps the parser simple while still letting the OS page in weight blocks on demand. A common pattern is: read metadata into ordinary managed objects, store tensor descriptors with offsets into the file, then create a view accessor or mapped stream that serves slices to the math layer.

public sealed record GgufTensorInfo(
    string Name,
    TensorKind Kind,
    long[] Shape,
    long Offset);

public sealed class GgufReader
{
    public static GgufModel Open(string path)
    {
        using var fs = File.OpenRead(path);
        using var br = new BinaryReader(fs, Encoding.UTF8, leaveOpen: true);

        uint magic = br.ReadUInt32();
        uint version = br.ReadUInt32();
        if (magic != 0x46554747) throw new InvalidDataException("Not GGUF");
        if (version < 2 || version > 3) throw new NotSupportedException();

        ulong tensorCount = br.ReadUInt64();
        ulong kvCount = br.ReadUInt64();

        var metadata = ReadMetadata(br, kvCount);
        var tensors  = ReadTensorInfos(br, tensorCount);
        var mmf = MemoryMappedFile.CreateFromFile(path, FileMode.Open);
        return new GgufModel(path, metadata, tensors, mmf);
    }
}

Once the tensor map exists, model construction is mostly name resolution. Pull the expected metadata values such as vocabulary size, hidden size, layer count, head count, and RMSNorm epsilon. Then resolve tensor descriptors into per-layer objects. If the model stores some tensors in float form and others in quantized block form, keep that distinction explicit instead of prematurely normalizing everything to float arrays. Preserving the quantized representation is what lets you avoid blowing up RAM usage.

The hot path lives in math helpers that accept ReadOnlySpan<float> or raw Span<float> buffers. This keeps call overhead low and avoids unnecessary allocations. Vector<float> gives you an easy first layer of SIMD without committing to architecture-specific code immediately.

public static void MatMul(
    Span<float> output,
    ReadOnlySpan<float> matrix,
    ReadOnlySpan<float> vector,
    int rows,
    int cols)
{
    int width = Vector<float>.Count;
    for (int r = 0; r < rows; r++)
    {
        var acc = Vector<float>.Zero;
        int rowBase = r * cols;
        int c = 0;
        for (; c <= cols - width; c += width)
        {
            var w = new Vector<float>(matrix.Slice(rowBase + c, width));
            var x = new Vector<float>(vector.Slice(c, width));
            acc += w * x;
        }

        float sum = 0f;
        for (int i = 0; i < width; i++) sum += acc[i];
        for (; c < cols; c++) sum += matrix[rowBase + c] * vector[c];
        output[r] = sum;
    }
}

When you want more control, System.Runtime.Intrinsics opens the door to Avx, Avx2, Fma, and AdvSimd. The usual strategy is to keep a scalar reference path, a portable Vector<T> path, and then one or two specialized intrinsic paths guarded by IsSupported. That gives you reach across x64 and ARM64 while still letting fast machines accelerate the inner loop.

Quantized tensors fit naturally into C# if you treat them as raw bytes plus typed decode routines. The dequantizer should not materialize a full floating-point matrix unless absolutely necessary. Instead, decode just enough of a block to feed the current dot product. Here is a simplified pattern:

public static float DotQ40(
    ReadOnlySpan<byte> block,
    ReadOnlySpan<float> x)
{
    float scale = BitConverter.Int16BitsToHalf(
        BinaryPrimitives.ReadInt16LittleEndian(block)).ToSingle();

    float acc = 0f;
    for (int i = 0; i < 16; i++)
    {
        byte q = block[2 + i];
        int lo = (q & 0x0F) - 8;
        int hi = (q >> 4) - 8;
        acc += (lo * scale) * x[i * 2 + 0];
        acc += (hi * scale) * x[i * 2 + 1];
    }
    return acc;
}

RMSNorm, RoPE, and softmax all translate directly into idiomatic C#. MathF covers the scalar math, Span<float> handles slicing, and the code remains readable enough to cross-check against a Python reference. A surprisingly effective debugging strategy is to implement a “dump layer 0 token 3” mode that prints the first eight values after each operator so you can compare them line-by-line against a known-good implementation.

One pleasant extra in C# is asynchronous I/O. You generally do not need async in the inner inference loop, but it can be useful while staging auxiliary artifacts or pre-reading non-tensor metadata. A loader that reads the GGUF header and tokenizer data using Memory<byte> buffers can overlap some disk operations without making the math code itself asynchronous.

dotnet publish -c Release -r linux-x64 --self-contained

The generation loop then becomes straightforward. Encode the prompt, prefill the cache token by token, run single-token decode for each new step, apply temperature and top-k, append the sampled token, and stream detokenized text to the console. The result feels high-level compared with C, but the core numerical work remains explicit enough that performance tuning is still meaningful. For many engineers, C# hits a useful balance: vastly more ergonomic than C, yet still close enough to the metal to build a legitimate CPU inference engine.

Why PowerShell is slow, strange, and still educationally excellent

A pure PowerShell transformer runtime is not a practical serving stack, but it is a powerful teaching tool because it removes all mystique. PowerShell can open binary files via .NET, store arrays, call [Math] for floating-point operations, and print tokens as they stream out. It will be slow. There is no hidden SIMD miracle. Most loops will run on scalar [double] or [float] values interpreted by the PowerShell engine. But the algorithm still works, which is the important lesson.

The basic approach is to read the GGUF file with [System.IO.BinaryReader], parse enough metadata to locate tensors and tokenizer entries, and store model weights as flat arrays. For a tiny educational model you may even choose to convert quantized tensors to [float[]] once during load. That wastes memory but keeps the math code legible. Once the reference version works, you can selectively move hot loops into inline C# via Add-Type or into precompiled helper assemblies while keeping the orchestration in PowerShell.

That educational value is why this variant matters. You are not trying to beat C. You are trying to demonstrate that transformer inference is ultimately file parsing, array math, caching, and sampling. PowerShell makes every step painfully visible.

function Read-UInt32([System.IO.BinaryReader]$br) {
    return $br.ReadUInt32()
}

function Read-String([System.IO.BinaryReader]$br) {
    $len = $br.ReadUInt64()
    $bytes = $br.ReadBytes([int]$len)
    return [System.Text.Encoding]::UTF8.GetString($bytes)
}

function Invoke-MatMul([float[]]$matrix, [float[]]$vector, [int]$rows, [int]$cols) {
    $out = [float[]]::new($rows)
    for ($r = 0; $r -lt $rows; $r++) {
        [double]$sum = 0.0
        $base = $r * $cols
        for ($c = 0; $c -lt $cols; $c++) {
            $sum += $matrix[$base + $c] * $vector[$c]
        }
        $out[$r] = [float]$sum
    }
    return $out
}

function Invoke-RmsNorm([float[]]$x, [float[]]$gain, [double]$eps) {
    [double]$sumsq = 0.0
    for ($i = 0; $i -lt $x.Length; $i++) { $sumsq += $x[$i] * $x[$i] }
    $inv = 1.0 / [Math]::Sqrt(($sumsq / $x.Length) + $eps)
    $y = [float[]]::new($x.Length)
    for ($i = 0; $i -lt $x.Length; $i++) {
        $y[$i] = [float]($x[$i] * $inv * $gain[$i])
    }
    return $y
}

That is already enough to build the core of a tiny decoder-only runtime. You still need GGUF metadata resolution, RoPE, attention over the cache, and sampling, but the math structure mirrors the C and C# variants almost line for line. The main difference is throughput. Nested loops that are merely “boring” in C become glacial in PowerShell once hidden widths and layer counts grow.

The generation loop can still be pleasantly direct. The following example assumes you already have a loaded model object with helper methods for tokenization, forward stepping, and detokenization. The point is not that PowerShell magically makes those helpers fast; the point is that the control flow is exactly the same control flow every serious engine uses.

$model = Load-GgufModel ".\model.gguf"
$tokens = [System.Collections.Generic.List[int]]::new()
$tokens.AddRange((Tokenize $model "What is the capital of France?"))
$eos = $model.EosTokenId

for ($i = 0; $i -lt $tokens.Count; $i++) {
    [void](Forward-Step $model $tokens[$i] $i)
}

$position = $tokens.Count
while ($true) {
    $logits = Forward-Step $model $tokens[$position - 1] ($position - 1) -DecodeOnly
    $next = Sample-TopK -Logits $logits -Temperature 0.7 -TopK 40
    $tokens.Add($next)

    $piece = Detokenize-Token $model $next
    Write-Host -NoNewline $piece

    if ($next -eq $eos) { break }
    $position++
}

In a real educational script you would expand Forward-Step into explicit functions: lookup embedding row, apply attention norm, project to Q/K/V, rotate Q and K, append to the cache, compute attention scores, apply softmax, mix values, add the residual, run the FFN, then finally produce logits. That sounds large, but each function is small. The bulk of the source code is array indexing, not new mathematics.

Where PowerShell becomes surprisingly useful is in experimentation. You can quickly dump tensor metadata, inspect tokenizer maps, compare the first eight logits against a Python reference, or prototype sampling behavior without recompiling anything. If you are about to write a C version, a PowerShell script can serve as an executable notebook: ugly, slow, but extremely transparent.

The biggest optimization trick is to push the hot loops into .NET while keeping PowerShell as the shell. Add-Type lets you embed small C# helper classes for matmul or quantized dot products. Alternatively, you can compile a tiny assembly once and call into it. That hybrid preserves the inspectability of PowerShell but removes the catastrophic overhead from the deepest loops.

So yes: PowerShell can run transformer inference on CPU. The honest performance expectation is terrible, but the intellectual payoff is excellent. If you can make the prompt "What is the capital of France?" produce plausible output in PowerShell, you have internalized the algorithm at a much deeper level than most framework users ever need to.

Why Bash is absurd, and why that absurdity teaches something real

A Bash implementation is not a recommendation. It is a stress test for your own understanding. If you can express the forward loop, token stepping, and a toy quantized matmul using shell tools, then you truly understand which parts of transformer inference are essential and which parts are ergonomics or performance engineering. Bash strips the stack down to raw orchestration, pipes, and small text-processing helpers.

The price is obvious. Bash has no native tensors, no native floating-point arithmetic, and no pleasant binary parsing API. You end up leaning on od, xxd, dd, awk, sed, and maybe bc. A “vector” is often a whitespace-separated line. A “matrix” is an array of lines or a text file on disk. The result is wildly slow, but slowness is not the point here.

For GGUF parsing, the shell can read headers using byte offsets. This is clumsy but real: od or xxd can extract integers, and small helper functions can advance an offset through the file as fields are decoded. Tokenizer handling is hardest in pure Bash; for a serious experiment you would usually preprocess the tokenizer into a simpler TSV form once and then let Bash orchestrate a greedy match loop.

read_u32() {
  local file="$1" offset="$2"
  od -An -t u4 -N 4 -j "$offset" "$file" | awk '{print $1}'
}

read_f32() {
  local file="$1" offset="$2"
  od -An -t f4 -N 4 -j "$offset" "$file" | awk '{print $1}'
}

matmul_awk() {
  local matrix_file="$1" vector_file="$2" rows="$3" cols="$4"
  awk -v rows="$rows" -v cols="$cols" '
    NR==FNR { x[FNR-1]=$1; next }
    {
      r = int((NR-1)/cols)
      c = (NR-1)%cols
      acc[r] += $1 * x[c]
    }
    END {
      for (i=0; i<rows; i++) print acc[i]+0.0
    }' "$vector_file" "$matrix_file"
}

That snippet is already enough to communicate the idea. A matrix stored as one scalar per line can be multiplied by a vector stored the same way, and the result can flow into the next stage. It is computationally ridiculous, but mathematically faithful. Add an RMSNorm function in awk, a softmax in awk, and a few helper scripts for cache management, and you have a real if comically slow inference pipeline.

The generation loop itself can stay short because Bash is good at orchestration. Again, imagine the tokenizer and weight-extraction helpers already exist. The shell’s job is to call them in the right order and append the sampled token each time.

MODEL="model.gguf"
TOKENS=$(tokenize "$MODEL" "What is the capital of France?")
EOS=$(get_eos_id "$MODEL")
POS=0

for t in $TOKENS; do
  forward_prefill "$MODEL" "$t" "$POS" >/dev/null
  POS=$((POS + 1))
  LAST="$t"
done

while true; do
  LOGITS=$(forward_decode "$MODEL" "$LAST" $((POS - 1)))
  NEXT=$(printf '%s\n' "$LOGITS" | sample_topk_awk 0.7 40)
  detokenize_piece "$MODEL" "$NEXT"
  [ "$NEXT" = "$EOS" ] && break
  LAST="$NEXT"
  POS=$((POS + 1))
done

In a fully pure Bash version, forward_prefill, forward_decode, and sample_topk_awk would themselves be shell-plus-awk programs manipulating text files or pipes. That is possible, but this is also the place where practicality re-enters the conversation. A very sensible compromise is to keep Bash as the conductor and delegate the numerical hot loop to a tiny C helper executable. Then Bash still exposes the sequence logic, file paths, and control flow, while the part that actually multiplies numbers runs at a saner speed.

The reason this variant is worth discussing in a technical book is not that anyone should deploy it. It is that Bash annihilates excuses. Once you see the algorithm expressed in shell terms, it becomes obvious that transformer inference is a deterministic dataflow problem. All the intimidating machinery in mature runtimes is there for performance, safety, and model breadth, not because the core idea is fundamentally inaccessible.

There is also one practical niche: Bash can be a thin glue layer around a family of tiny purpose-built helpers. One helper might extract GGUF metadata, another might perform quantized dot products, another might tokenize. The shell script then becomes a reproducible experimental harness for validating the whole loop end-to-end. Even if you never ship the Bash version, building it teaches you where every dependency boundary truly sits.