Loops are having a moment. The for and the while are the first control structures every programmer meets — and after decades of us hiding them behind map, reduce, vectorized array ops, and dataframe one-liners, they are suddenly the talk of the town again. Why? Because agentic AI is, underneath all the branding, a loop: perceive, decide, act, observe, repeat until done. while not solved: take_a_step(). The most autonomous software on earth in 2026 is a while True with good taste. The oldest primitive in computer science is cool again.

So let us talk about a loop. Specifically, the single most common way a for loop quietly destroys performance when you port code from a compiled language to an interpreted one — told through a real kernel in UXarray that computes the area of every face in an unstructured climate grid, ported from the C++ in TempestRemap’s GridElements.cpp.

The punchline up front: the port was faithful — line for line, it matched the C++ — and that faithfulness was exactly the bug. The same loop structure that is near-optimal in C++ is near-pessimal in Python, for reasons that are pure computer-science fundamentals.

The C++ that everyone copied

Here is the shape of the TempestRemap area routine, lightly trimmed:

// Mesh::CalculateFaceAreas — loop over every face
for (int i = 0; i < faces.size(); i++) {
    vecFaceArea[i] = CalculateFaceArea(faces[i], nodes);
}

Real CalculateFaceAreaQuadratureMethod(const Face & face, const NodeVector & nodes) {
    const int nOrder = 6;
    DataArray1D<double> dG;
    DataArray1D<double> dW;
    GaussQuadrature::GetPoints(nOrder, 0.0, 1.0, dG, dW);   // quadrature points/weights

    int nTriangles = face.edges.size() - 2;                // fan-triangulate the polygon
    for (int j = 0; j < nTriangles; j++) {
        Node node1 = nodes[face[0]];
        Node node2 = nodes[face[j+1]];
        Node node3 = nodes[face[j+2]];
        for (int p = 0; p < dW.GetRows(); p++)
        for (int q = 0; q < dW.GetRows(); q++) {
            dFaceArea += dW[p] * dW[q] * dJacobian;         // accumulate weighted area
        }
    }
    return dFaceArea;
}

This is good C++. It reads cleanly, and it runs fast. Look closely at three things it does, because each one is a trap in translation:

1. It calls GaussQuadrature::GetPoints inside the per-face routine — i.e. the quadrature weights are rebuilt for every single face.
2. It creates Node node1, node2, node3 inside the triangle loop — a fresh little struct per triangle.
3. It runs the whole thing as one serial for over faces.

In C++, none of these matter. And that is the entire lesson.

Why C++ shrugs

Cost of a loop iteration. In compiled C++, a for loop body is a handful of machine instructions. The loop counter lives in a register; the bounds check is one compare-and-branch the CPU’s branch predictor nails every time. There is no “cost of iterating” worth measuring. A billion iterations of trivial work is a fraction of a second.

Cost of Node node1 = nodes[...]. Node is a fixed-size value type — three doubles, 24 bytes, on the stack. “Allocating” it is subtracting from the stack pointer: one instruction, no bookkeeping, freed automatically when the scope ends. Creating three of them per triangle is free in any sense you can measure.

Cost of GetPoints per face. It fills two small stack/DataArray1D buffers with numbers. Cheap native code. Redundant across faces? Technically yes. Does anyone care at C++ speed? No — it is a few nanoseconds swallowed by a routine doing real arithmetic.

Cost of the serial loop. Even single-threaded, native throughput is high enough that a whole mesh finishes fast; and if you want more, OpenMP is one #pragma away.

So the C++ author made three “inefficiencies” that are not inefficiencies in C++. The compiled execution model makes per-iteration overhead, small stack allocations, and redundant cheap calls all effectively free. The code is clean and fast simultaneously. There was never any reason to optimize them away.

Why Python (and even Numba) punishes the exact same code

Now translate it faithfully. UXarray uses Numba — a just-in-time compiler that turns annotated Python into machine code — so a lot of the raw interpreter overhead is already gone. And yet the literal port was slow. Here is the faithful Python, and every line that was free in C++ now costs:

@njit(cache=True)
def calculate_face_area(x, y, z, quadrature_rule="triangular", order=4, ...):
    if quadrature_rule == "gaussian":
        dG, dW = get_gauss_quadrature_dg(order)   # (1) rebuilt PER FACE
    elif quadrature_rule == "triangular":
        dG, dW = get_tri_quadrature_dg(order)

    for j in range(num_triangles):
        node1 = np.array([x[0], y[0], z[0]])       # (2) HEAP allocation per triangle
        node2 = np.array([x[j+1], y[j+1], z[j+1]])
        node3 = np.array([x[j+2], y[j+2], z[j+2]])
        for p in range(len(dW)):
            if quadrature_rule == "gaussian":      # (3) STRING compare per point
                ...

# and the outer driver:
for face_idx in range(n_face):                     # (4) serial, one core
    face_area = calculate_face_area(...)

(1) The quadrature rebuild. In C++ this was a cheap native call you could ignore. In Python it constructs NumPy arrays and runs interpreted branching to pick the table — and it happens once per face. On a million-face mesh that is a million redundant table constructions. This is a textbook loop-invariant computation: the weights depend only on (rule, order), never on which face, so the result is identical every iteration. The compiler-optimization name for the fix is loop-invariant code motion — hoist it out, compute once, reuse.

(2) np.array([...]) per triangle. This is the one that really hurts, and it is the deepest C++/Python impedance mismatch. In C++, Node was 24 bytes on the stack — free. In Python/NumPy, np.array([...]) is a heap allocation: the runtime finds free memory, builds a full ndarray object (header, shape, strides, refcount, a separate data buffer), hands back a pointer, and later the garbage collector must reclaim it. There is no such thing as a cheap 3-element NumPy array. The faithful port turned a free stack value into an expensive heap object — and put it inside the innermost loop, so it pays that tax on every triangle of every face. (I wrote a whole post on this heap-vs-stack distinction — it is the same villain wearing a different hat.)

(3) if quadrature_rule == "gaussian" per quadrature point. A string comparison is not free: it compares characters. In C++ this branch was on a cheap value inside optimized code; here it runs on a Python string, deep inside the hottest loop, evaluated for every quadrature point of every triangle of every face. Replace the repeated string compare with a boolean computed once (strength reduction — swap an expensive operation for a cheap equivalent).

(4) The serial outer loop. Every face’s area is independent — no face reads another face’s result. This is the definition of an embarrassingly parallel problem. C++ left it serial because native serial was already fast enough; Python left it serial because the port was faithful. But an interpreted-then-JITed language needs the parallelism more, not less, to claw back the per-iteration cost.

The fix: change the loop, not the math

Not one line of the geometry — the spherical-triangle Jacobian, the quadrature accumulation — changed. Only the loop’s structure changed, applying four classic optimizations:

# (1)+(3): compute quadrature ONCE, outside the face loop; pick rule with a bool
if quadrature_rule == "gaussian":
    dG, dW = get_gauss_quadrature_dg(order); is_gaussian = True
else:
    dG, dW = get_tri_quadrature_dg(order);   is_gaussian = False

# (4): the faces are independent -> prange runs them across all CPU cores.
#      each iteration writes a distinct area[face_idx], so there is no data race.
for face_idx in prange(n_face):
    ...
    area[face_idx] = _face_area_from_quadrature(fx, fy, fz, dG, dW, is_gaussian, ...)

# (2): inside the kernel, hoist the shared apex node out of the triangle loop
node1 = np.array([x[0], y[0], z[0]])   # built once, not per triangle
for j in range(num_triangles):
    node2 = ...; node3 = ...

Four fundamentals, one loop:

• Loop-invariant code motion — quadrature weights out of the per-face path.
• Strength reduction — string compare → boolean.
• Hoisting — the shared fan-apex node built once, not per sub-triangle.
• Data parallelism — prange + parallel=True, safe because each iteration writes its own index (no shared mutation, no race).

Measured on a 5,400-face cubed-sphere grid, best-of-seven, same machine: 17.9 ms → 2.8 ms, a 6.3× speedup, with byte-identical areas (total still 4π on the unit sphere to 1e-10). On bigger meshes and more cores the prange term scales further.

The lesson for the CS grad

The instinct when porting is fidelity: match the source line by line, and trust that correct structure yields correct behavior. For correctness, fidelity is exactly right. For performance, fidelity is a trap, because performance is not a property of the code’s structure — it is a property of the execution model the structure runs on.

C++ and Python have opposite execution models in three ways that all converge inside a loop:

So the rules of thumb, learned the hard way:

1. A faithful port preserves the algorithm, not the performance. Re-derive the performance separately, on the target language’s terms.
2. Hunt for loop-invariants at the port boundary. Anything a C++ author left inside a loop “because it was cheap” may be expensive in Python. Hoist it.
3. Small fixed-size allocations are the silent killer. A C++ struct on the stack becomes a Python heap object. If it lives in a loop, it is a fire. Keep values as scalars in registers, or allocate once outside.
4. Independent loops want to be parallel. prange/vectorization is not premature optimization in interpreted-language ports — it is recovering ground the execution model took away.
5. Profile the loop, not the math. The temptation is to suspect the hard trigonometry. The cost was in the plumbing around it. Measure where time goes before touching what looks complicated.

The math was never the bottleneck. The loop was. Which feels right for 2026 — the year the loop, that first and oldest primitive, quietly runs both our agents and, if we are not careful, our performance regressions. Respect the loop. It is having its moment for a reason.

The area-kernel cleanup lives in UXarray (issue #1571); the original C++ is TempestRemap. For the heap-vs-stack half of this story, see the previous post.