I was born in 1982, a year that let me see, one after the other, the Commodore-64, the Amiga-600 and my first pc ever, an Olivetti-286 bought second or even third hand from a family friend. We are at the end of the ’80s, early ’90s at the latest, and I remember three things from that time:

  • computing was an “analogue” experience of sharing. We would all be together, now at my house, now at yours, trying to work out how these things worked.
  • a book on the C language, “Guida al Turbo C”, or something of the sort with the word Turbo in it. I remember buying it at a newsstand. I bought it because it came with floppy disks and because the cover had the words “Turbo” and “Programmazione” written on it. But above all “Turbo”.
  • after television ran one of the first news stories about “computer viruses”, my grandmother became very unhappy about leaving me in front of one of those machines, at the risk of catching some strange illness.

I can say that I have never abandoned C in thirty-odd years; even when I write in C++, everything stays very “C-ish”. If I write in Python I keep thinking about the chain reaction from the code to the interpreter, and I cannot help wondering which choice in the code is the more advantageous.

Meanwhile my way of programming built itself up over the years, punctuated by particular events. From the early days of slapping the keyboard trying, essentially, to make games, through to increasingly complicated programs.

The first evolution came when I read “C Interfaces and Implementations” (David Hanson) for the first time. I write “for the first time” because over the years I have “devoured” that book to the point where the pages have largely come away and the book has turned into more of a binder.

In those pages I learned what it means to model an idea onto code. The notion of creating “architectures” as a sum of atomic concepts.

The second great evolution I had when I began to study Algebra seriously.

Peano’s axioms are presented in many forms, but as time and reading went by, my favourite exposition became the one taken from an Italian textbook: “Aritmetica e Algebra” (Dikran Dikranjan, Maria Silvia Lucido). I give the axioms exactly as they appear there.

There exists a set $N$ and a map $s: N \rightarrow N$, called the successor, such that:

  1. $a \in N$
  2. if $n \in N$ then $s(n) \in N$ as well
  3. if $n \in N$, then $s(n) \neq a$
  4. $s$ is injective
  5. If the set $E$ contains $a$ and has the property that together with every $n \in E$ it also contains $s(n)$, then $N \subseteq E$.

What I like about this exposition is that the axioms are built without ever naming the natural numbers explicitly. Every set with these characteristics is in fact equivalent to $\mathbb{N}$.

The last proposition (the fifth) is nothing other than the principle of induction. An intuitive way of making it “visible” is to picture the game of toppling dominoes, the practice of lining up geometric tiles in sequence and then setting them off, creating a chain reaction. How can I say that the tiles, from the first to the last, will fall without failures? Thanks to the fact that the tiles are all of the same size and that the successor, in this case, is defined by a tile placed at a distance smaller than its own height.

We can therefore state that the set of domino tiles possesses the same properties as the set $N$, and that the definition of successor given above remains valid for a toppling chain made of infinitely many pieces.

The real turning point came when I read “Introduction to Algorithms” (Cormen, Leiserson et alii). Today for me, more than ever this book remains “The Book”. It tackles the subject formally and in pseudo code; it is a wonderful adventure through the recesses of code design.

In its opening pages it sets out the thing that changed my approach to programming for good: the loop invariant. The formal definition of what had until then been, for me, the plain for statement.

An iterative structure makes an algorithm correct if:

  • it is true at the first iteration of the loop (initialisation)
  • if it is true before an iteration of the loop, then it remains true before the next one (maintenance)
  • when the loop terminates, the invariant establishes the correctness of the algorithm

It is the text itself that declares the analogy with mathematical induction, while the revolution it set off in my approach is a continuous push towards the search for an “algebraisation” of programming. By then I was starting to think of a program as a system in which to build well defined operations on data.

Just as programming was taking a step towards abstraction and the object oriented paradigm was establishing itself, I decided at a certain point to set off on a journey in entirely the opposite direction. I began to study the hardware (dear old Tanenbaum, another book I consider fundamental) and assembly. Computers were evolving and I wanted to reason in terms of algebraic structures, not classes.

Years later I found part of my way of thinking again in the ideas set out in the architectural paradigm known as Data-Oriented Design (DOD) (for anyone wishing to go deeper, I refer you to the near-manifesto of one of its bards, Richard Fabian, https://www.dataorienteddesign.com/dodbook/).

In Data-Oriented Design one does not start from the single element and from what it knows how to do, but from the set of elements that undergo the same transformation. Data structures and algorithms are designed to treat them as a single block. This is where I found again the way of reasoning of group theory. In a group an element exists only through the operation that binds it to the others, and what is studied is the pair set plus operation, closed upon itself.

Let me be clear, this is a methodological analogy, not a structural one. There is no mathematical formalisation here; I am talking about what shaped my approach to programming.

The last evolution was brought about by the arrival of AI programming assistants. My relationship with these tools (that is what we are talking about, and we should keep it in mind) settled over time, and it had its turning point the moment I realised that my greatest adversary all these years had always been them: my hands. However hard I tried, my hands could not keep up with the ideas.

And from this followed the question: can programming and writing code be separated?

That they are two distinct things is almost obvious. Programming, for me, has always been an intellectual act; it is culture, art if you like. Writing code is a physical act, an executive one. Less obvious is that they can be pulled apart, because we are inclined to hold that it is by writing that one discovers when an idea does not stand up, and in this sense writing is not only execution, it is also verification.

My answer is that what gets delegated is neither the programming nor the writing, it is the typing. The judgement on what is written stays where it was.

Adopting these new tools has meant more than one change for me, and has brought an unexpected variety of approaches. I have not stopped writing code; on the contrary, I have more time to study new solutions and new languages. What is more, speeding up the analysis of code in search of weak points has led to my ideas being called into question immediately.

These are no small matters. The path research $\rightarrow$ experimentation $\rightarrow$ implementation $\rightarrow$ failure $\rightarrow$ research $\rightarrow$ implementation $\rightarrow$ … meant a great deal of time before seeing results.

And it is not true that the whole process amounted to a form of learning anyway, because the parts of that process with real added value were diluted among many other attendant practices, often repetitive and mechanical. As for the question of whether I let the AI assistant run free or step in, that depends on many aspects; the notion of delegation is never a trivial one.

Let us take an example. For narrative purposes I will give only the essential blocks of the code generated by the AI, without the surrounding material.

We need to simulate the motion of a million particles.

A solution produced without any particular guidance (other than the choice of C++ as the language) led the model to sketch out a decidedly Object Oriented design.

class Particle {
public:
    Particle(double x, double y, double z,
             double vx, double vy, double vz)
        : x(x), y(y), z(z), vx(vx), vy(vy), vz(vz) {}

    void update(double dt, double t) {
        x += vx * dt;
        y += vy * dt;
        z += vz * dt;
    }

private:
    double x, y, z;
    double vx, vy, vz;
    double temperature = 0.0;
};

And here is the computational part.

...
Particle* particelle = new Particle[N];
...
for (int i = 0; i < N; i++) 
    particelle[i].update(dt, t);

It makes sense, but what will the hardware have to do in order to execute this code?

The Particle class occupies 56 bytes in memory (six double values for position and velocity, one for the temperature). In the dynamically allocated memory Particle particelle[N], particle i begins at offset i × 56. The problem is that 56 is not a divisor of 64 (the size of a cache line). The remainder (i × 56) mod 64 changes with every particle and comes back into alignment only after 8 particles, because the least common multiple of 56 and 64 is 448 bytes, that is exactly 7 cache lines.

Working out the offsets of the first eight:

$$ \begin{aligned} \text{particle 0: offset 000 } &\rightarrow \text{ fits within one cache line} \\ \text{particle 1: offset 056 } &\rightarrow \text{ straddles two lines}\\ \text{particle 2: offset 112 } &\rightarrow \text{ straddles}\\ \text{ particle 3: offset 168 } &\rightarrow \text{ straddles}\\ \text{ particle 4: offset 224 } &\rightarrow \text{ straddles}\\ \text{ particle 5: offset 280 } &\rightarrow \text{ straddles}\\ \text{ particle 6: offset 336 } &\rightarrow \text{ straddles}\\ \text{ particle 7: offset 392 } &\rightarrow \text{ fits entirely within one line}\\ \end{aligned} $$

Six particles out of eight have their fields split across two different cache lines. It follows that update() has to fetch two 64 byte lines from memory instead of one. For every representation, moreover, it holds that 8 bytes out of 56 are not used in the computation (temperature).

The reasoning above is already optimistic, because it assumes that the array begins exactly at the start of a cache line. In practice, however, this assumption is not guaranteed. The statement new Particle[N] allocates memory through operator new[], which guarantees only the platform’s default alignment, typically 16 bytes on 64 bit systems. It does not guarantee a 64 byte alignment, that is, alignment to the size of a cache line.

A simple experimental check confirms this behaviour. In one of the runs, the address returned was indeed aligned to 16 bytes, but its remainder with respect to 64 was 16 (address mod 64 = 16). This value depends on the allocator and on the platform, not on the individual run, and in any case it is not under the program’s control.

Consequently, the offsets calculated earlier (0, 56, 112, …) are not referred to the start of a cache line, but to the initial address actually returned by the allocator. The conclusion, however, does not change: 6 particles out of 8 still occupy two cache lines. All that changes is the exact position at which this pattern begins, which depends on the alignment assigned to the array.

One factor contributing to the slowdown of the loop is represented by cache misses. When the data requested by the processor are not present in the cache (L1, L2 or L3), the CPU is forced to retrieve them from main memory (RAM), whose access times are far higher. During this wait the processor remains largely idle, introducing numerous stall cycles that significantly reduce overall performance. If the loop accesses data with poor spatial or temporal locality, the number of cache misses grows and execution time can increase appreciably, even when the number of computational operations remains unchanged.

And so we have to explain to a model how the class necessarily has to be decomposed into something completely different from the notion of “encapsulation”.

struct alignas(32) vec3 {
    double x, y, z;
    double _pad;   
};

struct scalar{
  double t;  
};

Here alignas(32) tells the compiler that every vec3 must begin at a memory address that is a multiple of 32. The exact 32 bytes I already get from this, because the size of a type is always a multiple of its alignment. The _pad member is declared in order to make that padding explicit, padding which the compiler would otherwise add anyway, but without a name.

This little trick will allow us to have sets characterised by well defined operations that are also well projected onto the hardware.

Both from the point of view of memory allocation:

vec3* position  = static_cast<vec3*>(std::aligned_alloc(32, static_cast<size_t>(N) * sizeof(vec3)));

vec3* velocity = static_cast<vec3*>(std::aligned_alloc(32, N * sizeof(vec3)));

scalar* temperature = static_cast<scalar*>(std::malloc(N * sizeof(scalar)));

There is one detail, though, that should not be overlooked. aligned_alloc returns uninitialised memory, and if at initialisation time I write only x, y and z, the _pad field keeps whatever bits happened to be in that region of memory. In the scalar loop nothing happens, because the three components are updated one by one. In the vectorised version, on the other hand, the load reads four double values in one go and the operation updates them all together, hence those bits too.

The results for x, y and z remain correct, because the four values are processed independently of one another. The cost is hidden elsewhere: if those random bits happen to correspond to denormal values, the floating point unit pays a penalty that on many microarchitectures is a heavy one, and it pays it on a million particles at every step. It is also the worst kind of thing to diagnose, because it does not always show up: the pages the operating system hands over are zeroed, so at the first allocation _pad is zero and everything is fine. The residual bits appear only when the allocator returns memory that has already been used and freed, that is, in long runs and not in short tests. Code that gives the right numbers and every so often runs far slower than it should.

It is enough to zero the memory right after allocating it.

std::memset(position, 0, static_cast<size_t>(N) * sizeof(vec3));
std::memset(velocity, 0, static_cast<size_t>(N) * sizeof(vec3));

This arrangement brings me a further advantage, in that I can reuse the same code while exploiting the processor’s SIMD operations.

SIMD (Single Instruction, Multiple Data) is a processing paradigm that allows the processor to perform the same operation on several pieces of data at once. On modern x86 architectures, the AVX (Advanced Vector Extensions) instructions operate on vector registers 256 bits wide, enough to hold four double values (4 × 64 bits) in a single variable.

In the code, the variable used for this purpose is of type __m256d, a type defined in <immintrin.h> that represents an AVX register containing four double precision numbers. The data are loaded from the position or velocity array by means of the _mm256_load_pd instruction, processed simultaneously, and then stored back to memory with _mm256_store_pd. Precisely because _mm256_load_pd requires the memory address to be aligned to 32 bytes, the corresponding arrays are allocated with aligned_alloc(32, ...).

temperature is not subject to the same constraint, since in this loop it is never read in blocks of 4 with a vector load; a plain malloc is enough.

So I can now go from this:

void update(vec3* position, vec3* velocity, int n, double dt) {
    for (int i = 0; i < n; i++) {
        position[i].x += velocity[i].x * dt;
        position[i].y += velocity[i].y * dt;
        position[i].z += velocity[i].z * dt;
    }
}

To this:

void update_simd(vec3* position, vec3* velocity, int n, double dt) {
    const __m256d vdt = _mm256_set1_pd(dt);
    for (int i = 0; i < n; i++) {
        __m256d v = _mm256_load_pd(reinterpret_cast<double*>(&velocity[i])); 
        __m256d p = _mm256_load_pd(reinterpret_cast<double*>(&position[i]));
         p = _mm256_fmadd_pd(v, vdt, p);     
        _mm256_store_pd(reinterpret_cast<double*>(&position[i]), p);
    }
}

Or else I can decide to move the computation onto the GPU, and since I want to try something other than CUDA, I try OpenCL.

OpenCL (Open Computing Language) is an open standard, an alternative to Nvidia’s CUDA and, like the latter, designed to run parallel computations by exploiting different kinds of hardware devices, such as GPUs, CPUs and dedicated accelerators. The idea behind OpenCL is to separate the program into two parts: the kernel, that is the code executed in parallel on the compute device, and the host code, which runs on the CPU and takes care of preparing the data, transferring them to the GPU and launching the execution.

In the OpenCL execution model, the work is divided into many independent units called work-items. In our example, each work-item would therefore take care of updating a single particle, and hence of reading a block of 32 consecutive bytes corresponding to the particle’s entire vector (double4). This lets the GPU access the data efficiently, because each element already contains all the necessary information in a single structure.

It should be said that this is not the only route, nor even the most orthodox one. Data-Oriented Design normally answers the same problem with three separate arrays, one for x, one for y and one for z, where the padding disappears altogether whereas in vec3 there remain 8 bytes out of 32. If I chose vec3 all the same, it is because those 32 contiguous bytes serve me twice over, as the alignment required by _mm256_load_pd and as the double4 read by an OpenCL work-item. A single layout that I can reuse along both routes without rethinking the access pattern. Not the best solution in absolute terms, the compromise best suited to this case.

This is exactly where what I was saying at the beginning about my reconsideration of object oriented programming comes in. The DOD approach fitted perfectly into the way I had already developed my approach to programming. Quite apart from the fact that a 3D vector is in itself an algebraic structure, in this example it has a twofold nature in the way I understand it as an object projected onto the hardware: a set of blocks of bytes defined by a characteristic, and well defined operations that act upon them.

This is an example of what I try today to convey to a tool capable of speeding up my idea of “programming”. Sometimes I find myself writing detailed specifications, sometimes I find myself “discussing” step by step, starting from an idea written in pseudo-code or prototyped.

The fact remains that, whereas in an early phase I was constantly trying to review and to step in, I have now realised that what I am standing in front of is a sort of Mirror of Galadriel.

In “The Lord of the Rings”, this elvish artefact acts as a kind of psychic amplifier that responds directly to the nature, the desires and the destiny of whoever looks into it. It does not project fixed images, but tunes itself to the soul of the observer, showing unique visions based on their sensibility and on their role in Middle-earth.

Now, I do not mean to say that an LLM tunes itself to the soul of whoever uses it; far more simply, it reacts directly to the tokens we have decided to pass to it.

So in the end, when a coding assistant works badly, I believe the great Machine Learning classic always applies: garbage in, garbage out.