code/day71-mixed-precision/matmul_fp16.cuThis is the source used by the lesson and its recorded evidence. Compile commands and expected output live in the directory README.
// SPDX-License-Identifier: MIT
//
// Day 71: mixed precision. Day 16's tiled matmul three ways: FP32 all the
// way through, FP16 storage with FP32 accumulation, and FP16 storage with
// FP16 accumulation, each judged against a double-precision CPU reference
// computed from the original FP32 inputs. Charging the half paths for their
// own storage rounding is the point: the error being bounded is the whole
// pipeline's, not just the arithmetic's.
//
// The bound is the K-scaled rule from day 66: rtol = max(table, 4 * eps *
// sqrt(K)) with eps = 2^-mantissa-bits, and atol = rtol * max|reference|.
// The program prints the arithmetic next to every gate it applies.
//
// The FP16-accumulate kernel is shipped to be measured, not copied. Its
// error ordering against FP32 accumulation is reported as a prediction, not
// treated as correctness: cancellation can reverse a max-error ordering on
// a particular input even when both kernels are correct.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o matmul_fp16 matmul_fp16.cu
// Run: ./matmul_fp16
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cuda_bf16.h>
#include <cuda_fp16.h>
#include <cuda_runtime.h>
// The one error macro. This file is standalone, the way a Compiler Explorer
// embed is, so it carries its own verbatim copy. `err_` carries a trailing
// underscore so it cannot collide with a variable at the call site, and the
// do/while makes the macro one statement so it survives a braceless `if`.
#define CUDA_CHECK(call) \
do { \
cudaError_t err_ = (call); \
if (err_ != cudaSuccess) { \
std::fprintf(stderr, "CUDA error %s:%d: %s: %s\n", __FILE__, \
__LINE__, #call, cudaGetErrorString(err_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
// Day 16's tile: 16 x 16 is one 256-thread block, the course default. The
// FP16 tiles cost half the shared memory of the FP32 pair, 1.5 KiB total
// for all three kernels' worst case, nowhere near the T4's 48 KiB default.
constexpr int kTileDim = 16;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
// Square cases, all multiples of the tile, chosen for the two axes this day
// measures. K grows 8x across the ladder, so error growth with K is visible
// in one table. And the footprints straddle the T4's 4 MiB L2: three
// 256 x 256 buffers are 768 KiB and live in cache, while the 2048 case
// moves 48 MiB and has to face DRAM, where halving the bytes can pay.
constexpr int kNumCases = 3;
constexpr size_t kCaseSize[kNumCases] = {256, 1024, 2048};
constexpr size_t kMaxDim = 2048;
// Base rtol per dtype from the course tolerance table, and the explicit
// mantissa bit counts that set eps = 2^-p: 23 for FP32, 10 for FP16.
constexpr double kF32Rtol = 1e-5;
constexpr double kF16Rtol = 1e-2;
constexpr int kF32MantissaBits = 23;
constexpr int kF16MantissaBits = 10;
static_assert(kTileDim * kTileDim == 256,
"a 16 x 16 tile is one 256-thread block, the course default");
static_assert((kTileDim * kTileDim) % 32 == 0,
"block size must be a whole number of warps");
static_assert(kCaseSize[0] % kTileDim == 0 && kCaseSize[1] % kTileDim == 0 &&
kCaseSize[2] % kTileDim == 0,
"every case divides by the tile, so all three kernels time the "
"same clean geometry and the table compares formats, not "
"guards");
static_assert(kCaseSize[0] < kCaseSize[1] && kCaseSize[1] < kCaseSize[2] &&
kCaseSize[2] <= kMaxDim,
"cases ascend so the error-growth column reads as K grows, and "
"everything fits the one kMaxDim allocation");
static_assert(2 * kTileDim * kTileDim * sizeof(float) <= 48 * 1024,
"the FP32 tiles must fit the 48 KiB a block gets on a T4 "
"without the cudaFuncSetAttribute opt-in");
static constexpr size_t ceilDiv(size_t a, size_t b) {
return (a + b - 1) / b;
}
// Day 16's tiled kernel, unchanged: one 16 x 16 tile of C per block, both
// inputs staged through shared memory, FP32 everywhere. This is the
// baseline row of every table below.
//
// Memory: threadIdx.x is the fastest index, so a warp's half-row of the B
// tile fill reads 16 consecutive floats, 64 bytes. Nothing is laid out
// badly; the half kernels change the element width and nothing else.
//
// Launch assumption: exactly kTileDim x kTileDim threads per block and a
// grid that rounds up on both axes. Guards cover loads and the store, never
// the barrier.
__global__ void matmulTiledF32(const float* __restrict__ a,
const float* __restrict__ b,
float* __restrict__ c, size_t m, size_t n,
size_t k) {
__shared__ float tileA[kTileDim][kTileDim];
__shared__ float tileB[kTileDim][kTileDim];
const size_t col =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t row =
blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
const size_t tiles = (k + kTileDim - 1) / kTileDim;
float acc = 0.0f;
for (size_t tileIdx = 0; tileIdx < tiles; ++tileIdx) {
const size_t aCol = tileIdx * kTileDim + threadIdx.x;
const size_t bRow = tileIdx * kTileDim + threadIdx.y;
tileA[threadIdx.y][threadIdx.x] =
(row < m && aCol < k) ? a[row * k + aCol] : 0.0f;
tileB[threadIdx.y][threadIdx.x] =
(bRow < k && col < n) ? b[bRow * n + col] : 0.0f;
__syncthreads();
for (int p = 0; p < kTileDim; ++p) {
acc += tileA[threadIdx.y][p] * tileB[p][threadIdx.x];
}
__syncthreads();
}
if (row < m && col < n) {
c[row * n + col] = acc;
}
}
// The same product with __half storage and a float accumulator. Inputs and
// tiles are 16-bit; every value becomes float in a register at the moment
// it is used, so the arithmetic is the FP32 kernel's arithmetic exactly and
// the only new error is the storage rounding the inputs already carry.
//
// Memory: a warp's half-row of the B tile fill now reads 16 consecutive
// __half values, 32 bytes where the FP32 kernel read 64. Same addresses,
// same pattern, half the sectors: that halving is the entire speedup this
// kernel can offer, because it uses no tensor cores. Day 72 adds those.
//
// Launch assumption: identical to matmulTiledF32.
// snippet: half-kernel
__global__ void matmulF16StoreF32Acc(const __half* __restrict__ a,
const __half* __restrict__ b,
float* __restrict__ c, size_t m, size_t n,
size_t k) {
__shared__ __half tileA[kTileDim][kTileDim];
__shared__ __half tileB[kTileDim][kTileDim];
const size_t col =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t row =
blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
const size_t tiles = (k + kTileDim - 1) / kTileDim;
float acc = 0.0f;
for (size_t tileIdx = 0; tileIdx < tiles; ++tileIdx) {
const size_t aCol = tileIdx * kTileDim + threadIdx.x;
const size_t bRow = tileIdx * kTileDim + threadIdx.y;
tileA[threadIdx.y][threadIdx.x] =
(row < m && aCol < k) ? a[row * k + aCol] : __float2half(0.0f);
tileB[threadIdx.y][threadIdx.x] =
(bRow < k && col < n) ? b[bRow * n + col] : __float2half(0.0f);
__syncthreads();
for (int p = 0; p < kTileDim; ++p) {
acc += __half2float(tileA[threadIdx.y][p]) *
__half2float(tileB[p][threadIdx.x]);
}
__syncthreads();
}
if (row < m && col < n) {
c[row * n + col] = acc;
}
}
// end snippet
// The cautionary tale: __half storage and a __half accumulator. Every
// __hfma rounds the running sum to 10 mantissa bits, so once the partial
// sum is large the terms being added land between representable values and
// the error grows with K. The output is widened to float only at the store,
// after the damage is done.
//
// Memory: identical to matmulF16StoreF32Acc, which is what makes the error
// column of the table attributable to the accumulator alone.
//
// Launch assumption: identical to matmulTiledF32.
__global__ void matmulF16StoreF16Acc(const __half* __restrict__ a,
const __half* __restrict__ b,
float* __restrict__ c, size_t m, size_t n,
size_t k) {
__shared__ __half tileA[kTileDim][kTileDim];
__shared__ __half tileB[kTileDim][kTileDim];
const size_t col =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t row =
blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
const size_t tiles = (k + kTileDim - 1) / kTileDim;
__half acc = __float2half(0.0f);
for (size_t tileIdx = 0; tileIdx < tiles; ++tileIdx) {
const size_t aCol = tileIdx * kTileDim + threadIdx.x;
const size_t bRow = tileIdx * kTileDim + threadIdx.y;
tileA[threadIdx.y][threadIdx.x] =
(row < m && aCol < k) ? a[row * k + aCol] : __float2half(0.0f);
tileB[threadIdx.y][threadIdx.x] =
(bRow < k && col < n) ? b[bRow * n + col] : __float2half(0.0f);
__syncthreads();
for (int p = 0; p < kTileDim; ++p) {
acc = __hfma(tileA[threadIdx.y][p], tileB[p][threadIdx.x], acc);
}
__syncthreads();
}
if (row < m && col < n) {
c[row * n + col] = __half2float(acc);
}
}
// CPU reference in double over the ORIGINAL float inputs, never the
// half-rounded copies. A reference computed from the rounded inputs would
// forgive the storage error, and the storage error is half of what this day
// bounds. Written for obvious correctness: plain loops, no blocking.
static void matmulCpu(const float* a, const float* b, double* c, size_t m,
size_t n, size_t k) {
for (size_t row = 0; row < m; ++row) {
for (size_t col = 0; col < n; ++col) {
double total = 0.0;
for (size_t p = 0; p < k; ++p) {
total += static_cast<double>(a[row * k + p]) *
static_cast<double>(b[p * n + col]);
}
c[row * n + col] = total;
}
}
}
// Inputs that FP16 cannot represent exactly, on purpose. Day 16 filled its
// matrices with small integers so any mismatch had to be an indexing bug;
// those integers are also exact in __half, so reusing them here would show
// a storage error of exactly zero and prove nothing. 0.013f and 0.017f are
// not representable in any binary format, so every conversion to half
// rounds and the error budget is actually spent. The periods 97 and 89
// divide none of the case sizes, and both ranges stay under 1, so no
// partial sum can approach __half's 65,504 ceiling even when the
// accumulator is half.
static void makeInputs(std::vector<float>* h_a, std::vector<float>* h_b,
size_t aElems, size_t bElems) {
for (size_t e = 0; e < aElems; ++e) {
(*h_a)[e] = static_cast<float>(static_cast<int>(e % 97) - 48) * 0.013f;
}
for (size_t e = 0; e < bElems; ++e) {
(*h_b)[e] = static_cast<float>(static_cast<int>(e % 89) - 44) * 0.017f;
}
}
// The K-scaled tolerance. Rounding errors across a K-term dot product are
// uncorrelated, so they grow like sqrt(K), not K; the factor 4 is slack for
// a different-but-valid summation order. Day 66 established the rule; this
// day is the first one whose eps makes it bite.
// snippet: tolerance
static double kScaledRtol(double tableRtol, int mantissaBits, size_t kDim) {
const double eps = std::ldexp(1.0, -mantissaBits);
const double grown = 4.0 * eps * std::sqrt(static_cast<double>(kDim));
return grown > tableRtol ? grown : tableRtol;
}
// end snippet
// Worst error as a fraction of the gate |got - ref| <= atol + rtol * |ref|.
// A fraction above 1 anywhere is a failure, and the caller gets the first
// offending index so the report can name a row and a column. A non-finite
// output fails outright: inf here means an overflow bug, and NaN is never
// close enough.
static double worstGateFraction(const float* got, const double* want, size_t n,
double rtol, double atol, size_t* firstBad) {
double worst = 0.0;
*firstBad = n;
for (size_t i = 0; i < n; ++i) {
if (!std::isfinite(got[i])) {
*firstBad = i;
return HUGE_VAL;
}
const double err = std::fabs(static_cast<double>(got[i]) - want[i]);
const double frac = err / (atol + rtol * std::fabs(want[i]));
if (frac > worst) {
worst = frac;
}
if (frac > 1.0 && *firstBad == n) {
*firstBad = i;
}
}
return worst;
}
static double maxAbsError(const float* got, const double* want, size_t n) {
double worst = 0.0;
for (size_t i = 0; i < n; ++i) {
const double err = std::fabs(static_cast<double>(got[i]) - want[i]);
if (err > worst) {
worst = err;
}
}
return worst;
}
static double maxAbs(const double* x, size_t n) {
double worst = 0.0;
for (size_t i = 0; i < n; ++i) {
if (std::fabs(x[i]) > worst) {
worst = std::fabs(x[i]);
}
}
return worst;
}
// Times a launch with CUDA events and returns the mean milliseconds per run.
//
// This is the one template and the one lambda allowed in module 1 to 3 code.
// Copy it verbatim; the alternative is a second copy of the event
// boilerplate, which is how a warm-up goes missing from one of them.
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
// Warm up this kernel, not just the first kernel in the program. Lazy
// module loading has been the default since CUDA 12.2 on Linux, so the
// first launch of each kernel pays its own load.
for (int i = 0; i < kWarmupRuns; ++i) {
launch();
}
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < kTimedRuns; ++i) {
launch();
}
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError());
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / kTimedRuns;
}
// What this card's tensor cores can execute, from the compute capability
// read at run time. The floors: FP16 mma.sync.m16n8k8 needs sm_75; BF16 and
// TF32 mma need sm_80; E4M3 and E5M2 mma need sm_89; E2M1 (FP4) mma needs
// the architecture-specific sm_120a target (datacenter Blackwell reaches
// FP4 through tcgen05 instead, which no runtime CC test can promise, so the
// FP4 line names the target rather than claiming this card runs it). All
// from the PTX ISA mma and tcgen05 Target ISA Notes, checked 2026-09-01.
// This program's own kernels use no tensor cores at all; the report exists
// so the transcript states what this card could and could not run.
// snippet: format-report
static void printFormatSupport(int major, int minor) {
const int cc = major * 10 + minor;
struct FormatFloor {
const char* name;
int minCc;
};
const FormatFloor floors[] = {
{"FP16 (mma.sync.m16n8k8)", 75}, {"BF16 (mma, sm_80)", 80},
{"TF32 (mma, sm_80)", 80}, {"FP8 E4M3/E5M2 (mma, sm_89)", 89},
{"FP4 E2M1 (mma, sm_120a)", 120},
};
std::printf("tensor-core format support at compute capability %d.%d:\n",
major, minor);
for (const FormatFloor& f : floors) {
std::printf(" %-28s %s\n", f.name,
cc >= f.minCc ? "yes" : "no on this card");
}
}
// end snippet
// Four conversions the page quotes, printed so the transcript carries them:
// the largest finite __half, the first float that rounds past it to inf,
// and one value stored in both 16-bit formats to show the budgets differ.
static void printFormatDemos() {
std::printf("largest finite __half: %.1f\n",
__half2float(__float2half(65504.0f)));
std::printf("__float2half(65520.0f) isinf: %d\n",
std::isinf(__half2float(__float2half(65520.0f))) ? 1 : 0);
std::printf("300000.0f in __half: %g (isinf %d)\n",
__half2float(__float2half(300000.0f)),
std::isinf(__half2float(__float2half(300000.0f))) ? 1 : 0);
std::printf("300000.0f in __nv_bfloat16: %.1f\n",
__bfloat162float(__float2bfloat16(300000.0f)));
std::printf("1.0007f in __half: %.10f\n",
__half2float(__float2half(1.0007f)));
std::printf("1.0007f in __nv_bfloat16: %.10f\n\n",
__bfloat162float(__float2bfloat16(1.0007f)));
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
std::printf("GPU: %s (compute capability %d.%d)\n\n", prop.name, prop.major,
prop.minor);
printFormatSupport(prop.major, prop.minor);
std::printf("\n");
printFormatDemos();
const size_t maxElems = kMaxDim * kMaxDim;
std::vector<float> h_a(maxElems);
std::vector<float> h_b(maxElems);
std::vector<__half> h_ah(maxElems);
std::vector<__half> h_bh(maxElems);
std::vector<float> h_c(maxElems);
std::vector<double> h_want(maxElems);
float* d_a = nullptr;
float* d_b = nullptr;
__half* d_ah = nullptr;
__half* d_bh = nullptr;
float* d_c = nullptr;
CUDA_CHECK(cudaMalloc(&d_a, maxElems * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_b, maxElems * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_ah, maxElems * sizeof(__half)));
CUDA_CHECK(cudaMalloc(&d_bh, maxElems * sizeof(__half)));
CUDA_CHECK(cudaMalloc(&d_c, maxElems * sizeof(float)));
// Failure is recorded rather than returned, so every path falls through
// to the one cleanup block below and no cudaMalloc escapes its cudaFree.
int badCase = kNumCases;
const char* badKernel = nullptr;
size_t badIndex = 0;
for (int caseIdx = 0; caseIdx < kNumCases; ++caseIdx) {
const size_t size = kCaseSize[caseIdx];
const size_t elems = size * size;
const size_t outBytes = elems * sizeof(float);
makeInputs(&h_a, &h_b, elems, elems);
// The half copies are made on the host with the documented
// conversion helper, so the rounding the kernels inherit is
// round-to-nearest-even and reproducible anywhere.
for (size_t e = 0; e < elems; ++e) {
h_ah[e] = __float2half(h_a[e]);
h_bh[e] = __float2half(h_b[e]);
}
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), elems * sizeof(float),
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), elems * sizeof(float),
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_ah, h_ah.data(), elems * sizeof(__half),
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_bh, h_bh.data(), elems * sizeof(__half),
cudaMemcpyHostToDevice));
matmulCpu(h_a.data(), h_b.data(), h_want.data(), size, size, size);
const double refMax = maxAbs(h_want.data(), elems);
const double f32Rtol = kScaledRtol(kF32Rtol, kF32MantissaBits, size);
const double f16Rtol = kScaledRtol(kF16Rtol, kF16MantissaBits, size);
const double f32Atol = f32Rtol * refMax;
const double f16Atol = f16Rtol * refMax;
std::printf("case %zux%zux%zu, max |reference| = %.4f\n", size, size,
size, refMax);
std::printf(
" f32 gate: |got-ref| <= %.3e + %.3e*|ref|"
" (rtol = max(1e-5, 4*2^-23*sqrt(K)))\n",
f32Atol, f32Rtol);
std::printf(
" f16 gate: |got-ref| <= %.3e + %.3e*|ref|"
" (rtol = max(1e-2, 4*2^-10*sqrt(K)))\n",
f16Atol, f16Rtol);
const dim3 block(kTileDim, kTileDim);
const dim3 grid(static_cast<unsigned int>(
ceilDiv(size, static_cast<size_t>(kTileDim))),
static_cast<unsigned int>(
ceilDiv(size, static_cast<size_t>(kTileDim))));
// Three passes over one scaffold: launch, gate, time. The half
// kernels share a correctness gate; their error ordering is reported
// after the table as a prediction for the hardware run to judge.
struct Path {
const char* name;
int gateIsF16;
};
const Path paths[] = {
{"matmulTiledF32", 0},
{"matmulF16StoreF32Acc", 1},
{"matmulF16StoreF16Acc", 1},
};
double absErr[3] = {0.0, 0.0, 0.0};
float meanMs[3] = {0.0f, 0.0f, 0.0f};
std::printf(" %-22s %13s %10s %10s %8s\n", "kernel", "max abs err",
"gate frac", "time (ms)", "vs f32");
for (int pathIdx = 0; pathIdx < 3 && badKernel == nullptr; ++pathIdx) {
// The output is cleared before each correctness launch so a
// kernel that skips elements cannot inherit a right answer.
CUDA_CHECK(cudaMemset(d_c, 0, outBytes));
switch (pathIdx) {
case 0:
matmulTiledF32<<<grid, block>>>(d_a, d_b, d_c, size, size,
size);
break;
case 1:
matmulF16StoreF32Acc<<<grid, block>>>(d_ah, d_bh, d_c, size,
size, size);
break;
default:
matmulF16StoreF16Acc<<<grid, block>>>(d_ah, d_bh, d_c, size,
size, size);
break;
}
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_c.data(), d_c, outBytes, cudaMemcpyDeviceToHost));
const double rtol = paths[pathIdx].gateIsF16 ? f16Rtol : f32Rtol;
const double atol = paths[pathIdx].gateIsF16 ? f16Atol : f32Atol;
size_t firstBad = elems;
const double frac = worstGateFraction(h_c.data(), h_want.data(),
elems, rtol, atol, &firstBad);
absErr[pathIdx] = maxAbsError(h_c.data(), h_want.data(), elems);
if (firstBad != elems) {
badCase = caseIdx;
badKernel = paths[pathIdx].name;
badIndex = firstBad;
break;
}
switch (pathIdx) {
case 0:
meanMs[0] = timeKernel([&] {
matmulTiledF32<<<grid, block>>>(d_a, d_b, d_c, size,
size, size);
});
break;
case 1:
meanMs[1] = timeKernel([&] {
matmulF16StoreF32Acc<<<grid, block>>>(d_ah, d_bh, d_c,
size, size, size);
});
break;
default:
meanMs[2] = timeKernel([&] {
matmulF16StoreF16Acc<<<grid, block>>>(d_ah, d_bh, d_c,
size, size, size);
});
break;
}
std::printf(" %-22s %13.3e %10.4f %10.4f %8.2f\n",
paths[pathIdx].name, absErr[pathIdx], frac,
meanMs[pathIdx],
static_cast<double>(meanMs[0]) / meanMs[pathIdx]);
}
if (badKernel != nullptr) {
break;
}
// Cancellation can make either valid accumulation order's max error
// smaller on one input. Keep the ratio visible without turning the
// lesson's prediction into a program failure.
std::printf(
" observed f16-acc/f32-acc max-error ratio at K=%zu: %.3g\n\n",
size, absErr[2] / absErr[1]);
}
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_ah));
CUDA_CHECK(cudaFree(d_bh));
CUDA_CHECK(cudaFree(d_c));
if (badKernel != nullptr) {
const size_t n = kCaseSize[badCase];
std::fprintf(stderr,
"%s failed its gate at %zu (row %zu, col %zu) on case "
"%zux%zux%zu\n",
badKernel, badIndex, badIndex / n, badIndex % n, n, n, n);
return EXIT_FAILURE;
}
std::printf("all %d cases passed their gates\n", kNumCases);
return EXIT_SUCCESS;
}