Sometimes, we come across a youtube video showcasing some true programming black magic, something that leaves us astonished, and the only thing we can do is be left scratching our head. Happened to me, at least.

Today, we’ll dive into one of these videos, and I’ll try to explain just what is going on.
Here’s the stuff:

If you didn’t quite catch what is happening, the source code of the program is this:

#  include<stdio.h>//  .IOCCC                                         Fluid-  #
#  include <unistd.h>  //2012                                         _Sim!_  #
#  include<complex.h>  //||||                     ,____.              IOCCC-  #
#  define              h for(                     x=011;              2012/*  #
#  */-1>x              ++;)b[                     x]//-'              winner  #
#  define              f(p,e)                                         for(/*  #
#  */p=a;              e,p<r;                                        p+=5)//  #
#  define              z(e,i)                                        f(p,p/*  #
## */[i]=e)f(q,w=cabs  (d=*p-  *q)/2-     1)if(0  <(x=1-      w))p[i]+=w*/// ##
   double complex a [  97687]  ,*p,*q     ,*r=a,  w=0,d;    int x,y;char b/* ##
## */[6856]="\x1b[2J"  "\x1b"  "[1;1H     ", *o=  b, *t;   int main   (){/** ##
## */for(              ;0<(x=  getc (     stdin)  );)w=x  >10?32<     x?4[/* ##
## */*r++              =w,r]=  w+1,*r     =r[5]=  x==35,  r+=9:0      ,w-I/* ##
## */:(x=              w+2);;  for(;;     puts(o  ),o=b+  4){z(p      [1]*/* ##
## */9,2)              w;z(G,  3)(d*(     3-p[2]  -q[2])  *P+p[4      ]*V-/* ##
## */q[4]              *V)/p[  2];h=0     ;f(p,(  t=b+10  +(x=*p      *I)+/* ##
## */80*(              y=*p/2  ),*p+=p    [4]+=p  [3]/10  *!p[1])     )x=0/* ##
## */ <=x              &&x<79   &&0<=y&&y<23?1[1  [*t|=8   ,t]|=4,t+=80]=1/* ##
## */, *t              |=2:0;    h=" '`-.|//,\\"  "|\\_"    "\\/\x23\n"[x/** ##
## */%80-              9?x[b]      :16];;usleep(  12321)      ;}return 0;}/* ##
####                                                                       ####
###############################################################################
**###########################################################################*/

This can be compiled with

gcc -DG=1 -DP=4 -DV=8 endoh1.c -o endoh1

The source code is then fed as an input to the program, and we can see it… melting? What is going on??

What the program does is, it simulates the motion of a fluid from an initial state it is given as an input, using a discrete simulation model called smoothed particle hydrodynamics (SPH).
Hashtags are immobile barriers and all the other non-space characters are counted as part of the fluid.

That’s some serious black magic.

Let’s see how it works!

Making it readable

First of all, we should clean up the source a little, so it’s in a normal form. The author made it spell fluid with hashtags around it to make it a cool input, but it’s pretty impossible to read!

Whitespace and indentation are ignored by the C compiler, this is what makes it possible to have something like that be a valid program.
Fixing those gives us this pretty piece of code:

#include<stdio.h>
#include <unistd.h>
#include<complex.h>

#define h for(x=011;2012-1>x++;)b[x]
#define f(p,e) for(p=a; e,p<r; p+=5)
#define z(e,i) f(p,p[i]=e) f(q,w=cabs(d=*p-*q)/2-1) if(0<(x=1-w)) p[i]+=w*

double complex a[97687], *p, *q, *r=a, w=0, d;
int x,y;
char b[6856]="\x1b[2J\x1b[1;1H     ",
    *o=b,
    *t;
int main(){
    for(;0<(x=getc(stdin));)
        w = x > 10 ?
            (32 < x ? 4[*r++=w,r]=w+1, *r=r[5]=x==35, r+=9:0),w-I
            :(x=w+2);
    for(;;puts(o),o=b+4){
        z(p[1]*9,2) w;
        z(G,3)(d*(3-p[2]-q[2])*P + p[4]*V - q[4]*V) / p[2];
        h=0;
        f(p,(t=b+10+(x=*p*I)+80*(y=*p/2),*p+=p[4]+=p[3]/10*!p[1]))
            x = 0<=x && x<79 && 0<=y && y<23 ?
                1[1[*t|=8,t]|=4,t+=80]=1,*t|=2
                :0;
        h=" '`-.|//,\\|\\_\\/\x23\n"[x%80-9 ? x[b] : 16];
        usleep(12321);
    }
    return 0;
}

Still pretty much unintelligible.

A frequently used trick here is comma-separated expressions. In C, when several expressions are chained together by commas, they behave as a single expression, in the sense that the whole is evaluated to only one value. Namely, the result of the last chained expression.
This saves space and also makes it harder to understand if used abundantly, which is exactly the point of this program.

Another obscure C feature is the fact that you can refer to array elements not only with array[index], but also with index[array].
Now that we understand this, let’s look at the expression

32 < x ? 4[*r++=w,r]=w+1, *r=r[5]=x==35, r+=9:0

In C, you can inline if statements like this: condition ? value_if_true : value_if_false.
So our condition is x > 32, if this is true, the expression evaluates to 4[*r++=w,r]=w+1, *r=r[5]=x==35, r+=9. If false, it gives 0 as a result.
But what is 4[*r++=w,r]=w+1, *r=r[5]=x==35, r+=9 ?
Let’s look at the single expressions:

4[*r++=w, r]=w+1 does 3 separate things. Following the order of operations, r[0] is set to w, r is incremented by 1, and then r[4] is set to w+1, since *r++=w, r evaluates to the last expression r.
*r=r[5]=x==35 sets both r[0] and r[5] to the boolean value of x==35. Since in our program x is a char, it checks if x is equal to #.
Finally, r+=9 increments r by 9. The new value of r is what the whole thing evaluates to.

Now, the resulting value never gets assigned to anything, since, as you can see if you read the code carefully, this whole thing is chained to another expression with a comma, so the result of the expression on that line is simply w-I.

Now that we got the concept, let’s see what the code looks like when written in a more standard fashion:

#include <stdio.h>
#include <unistd.h>
#include <complex.h>
#include <string.h>

#define CLEAR_SCREEN "\x1b[2J"
#define RESET_CURSOR "\x1b[1;1H"
#define LIQUID_CHARS " '`-.|//,\\|\\_\\/#\n"

/**
    Particles have five fields:
        0: position: complex, (row, -col)
        1: is_wall: int/bool, 1 if wall
        2: density
        3: force
        4: velocity

        Hypotesis: particles have radius 2, mass 1
**/
#define POS 0
#define IS_WALL 1
#define DENSITY 2
#define ACCEL 3
#define VEL 4

#define PARTICLE_LEN 5

#define P_RADIUS 2
#define P_MASS 1
#define P0 1.5

double complex array[97687], *p, *q, *arr_end=array, w=0, diff;
int x,y;
// clears the screen and sets cursor at upper left corner
char buffer[6856]=CLEAR_SCREEN RESET_CURSOR,
    *output=buffer,
    *t;

double dist(double complex p1, double complex p2){
    return cabs(p1-p2);
}

double complex W(double complex r){
    return (r/P_RADIUS-1)*(r/P_RADIUS-1);
}

int main(){
    // parse input
    for(;0<(x=getc(stdin));){
        // here, w is position: real is line num, negative img is col num
        if(x <= '\n')
            w = (int)w+2;
        else{
            if(x > ' '){
                // each char is 2 particles, one on top of the other
                double complex *part1 = arr_end;
                double complex *part2 = arr_end + PARTICLE_LEN;
                //set position
                part1[POS]=w;
                part2[POS]=w+1;
                // set is_wall
                part1[IS_WALL] = part2[IS_WALL] = (x=='#');
                // next pos in array
                arr_end += 2 * PARTICLE_LEN;
            }
            w-=I;
        }
    }

    while(1){
        // calculate density
        for(p=array; p<arr_end; p+=PARTICLE_LEN){
            // extra density if wall?
            p[DENSITY]=p[IS_WALL]*9;
            // integrate within domain, weighting with kernel funC
            for(q=array; q<arr_end; q+=PARTICLE_LEN){
                double rij = dist(p[POS], q[POS]);
                if(rij < P_RADIUS)
                    p[DENSITY] += P_MASS*W(rij);
            }
        }
        // calculate interaction force
        for(p=array; p<arr_end; p+=PARTICLE_LEN){
            // gravity component
            p[ACCEL]=G;
            for(q=array; q<arr_end; q+=PARTICLE_LEN){
                diff = p[POS] - q[POS];
                double rij = cabs(diff);
                if(rij < P_RADIUS){
                    double complex press = (p[DENSITY]+q[DENSITY]-2*P0)*P*diff;
                    double complex visc = (p[VEL] - q[VEL])*V;
                    p[ACCEL] +=
                        (1-rij/P_RADIUS)/ p[DENSITY] *
                        (press - visc);
                }
            }
        }
        // clear buffer
        for(x=10;x<2011;++x)
            buffer[x]=0;
        // fill buffer using marching cubes
        for(p=array; p<arr_end; p+=PARTICLE_LEN){
            // get coords
            x = -cimag(p[POS]);
            y = creal(p[POS])/2;
            // get pointer to buffer element
            t=buffer+strlen(CLEAR_SCREEN RESET_CURSOR)+x+80*y;
            // update position and velocity
            p[POS]+=p[VEL]+=p[ACCEL]/10*!p[IS_WALL];
            // render with marching cubes
            if(x>=0 && x<79 && y>=0 && y<23){
                t[0] |= 8;
                t[1] |= 4;
                t += 80;
                t[0] |= 2;
                t[1] = 1;
            }
        }
        for(x=10;x<2011;++x)
            buffer[x] = LIQUID_CHARS[x%80-9 ? buffer[x] : 16];
        usleep(12321);
        puts(output);
        output=buffer+strlen(CLEAR_SCREEN);
    }
    return 0;
}

There is clearly a lot going on here, which I am going to explain in the next sections.

The algorithm

So how does it work?
In short, it converts the input text into discrete particles, simulates their movement, and then draws them in some fancy way.

Input parsing

First, let’s look at how the particles are created from the input text. Here’s the code for that:

// parse input
for(;0<(x=getc(stdin));){
    // here, w is position: real is line num, negative img is col num
    if(x <= '\n')
        w = (int)w+2;
    else{
        // whitespace characters are skipped
        if(x > ' '){
            // each char is 2 particles, one on top of the other
            double complex *part1 = arr_end;
            double complex *part2 = arr_end + PARTICLE_LEN;
            //set position
            part1[POS]=w;
            part2[POS]=w+1;
            // set is_wall
            part1[IS_WALL] = part2[IS_WALL] = (x=='#');
            // next pos in array
            arr_end += 2 * PARTICLE_LEN;
        }
        w-=I;
    }
}

This program relies pretty heavily on the C complex number library. In general, any vector here is represented as a complex number. This makes it easier to manage the two x,y components together, instead of having two variables.
Here, the program scans the input line by line, updating a w complex variable which denotes the current position in the 2D space of the simulation.
Each character is 1 unit wide and 2 units high. So, each non-whitespace character gets converted into 2 particles.

The logic is pretty simple. The particles are stored in a complex number array, each with 5 components:

POS: the 2D position vector.
IS_WALL: either 0 or 1, the particle is immobile if it is a wall.
DENSITY: the density of the liquid where the particle is (will explain this later).
ACCEL: the 2D acceleration vector of the particle.
VEL: the 2D velocity vector.

Only POS and IS_WALL are initialized at this stage.
As you can see in the code, the position of a particle is taken from the w variable. IS_WALL is set to 1 if the corresponding character is a #.
A thing to note is that the x coordinate (left-right) is negative. That’s because of some operations the original code did later to extract the column number, not really necessary anymore, but I thought I’d leave it for reflect the original.

The simulation algorithm

This program uses Smoothed Particle Hydrodynamics (SPH) to simulate the behaviour of the fluid. The basic idea behind it is to have lots of discrete particles, with special rules governing their interaction, which on the bigger scale simulate a fluid.
There are three parameters governing this interaction:

G, the gravitational constant, here set to 1
P, the pressure, here set to 4
V, the viscosity, here set to 8 These are passed in at compile time as -DG=1 -DP=4 -DV=8.

The principle underlying SPH is that we can define the average of some function A in the point r as $\langle A\rangle = \int_{\Omega_\mathbb{r}} A(\mathbb{r'})W(\lvert\mathbb{r-r'}\rvert)d\mathbb{r'}$
$\Omega_\mathbb{r}$ is usually a circular area around the position $\mathbb{r}$ , and W is the so-called kernel function. It is used as a weight for the contributions in the integration area and has to satisfy several properties, but we’re not going to go into that here. In this implementation, it’s defined as $W(d) = (\frac{d}{R}-1)^2$ , where R is the radius of the integration area, specifically 2.

In practice, we are going to use a discrete version of this. Instead of integrating for each point in a radius of 2, we are going to integrate on each particle present in that area. So for each particle i we have $\langle A_i \rangle = \sum_{j\in \Omega_{\mathbb{r}_i}}\frac{m_j}{\rho_j}A(\mathbb{r}_j)W(\lvert \mathbb{r}_i-\mathbb{r}_j\rvert)$ , where $m_j$ is the mass of particle j (here always 1), and $\rho_j$ is its density.

In particular, it’s how we calculate the density of the fluid in various points (the location of each particle).
We have then that $\langle \rho_i \rangle = \sum_{j\in \Omega_{\mathbb{r}_i}}\frac{m_j}{\rho_j}\rho_j W(\lvert \mathbb{r}_i-\mathbb{r}_j\rvert) = \sum_{j\in \Omega_{\mathbb{r}_i}}m_j W(\lvert \mathbb{r}_i-\mathbb{r}_j\rvert)$ .
And that’s exactly what this part does:

// calculate density
for(p=array; p<arr_end; p+=PARTICLE_LEN){
    // extra density if wall
    p[DENSITY]=p[IS_WALL]*9;
    // integrate within domain, weighting with kernel func
    for(q=array; q<arr_end; q+=PARTICLE_LEN){
        double rij = dist(p[POS], q[POS]);
        if(rij < P_RADIUS)
            p[DENSITY] += P_MASS*W(rij);
    }
}

After that, we need to calculate the force each particle is subject to. There are two components:

gravitational force, which is a constant downward force of intensity $Gm_i$
interaction force, which is created by the, well, interactions with other particles.

To calculate the force on particle i exerted by particle j, SPH uses this formula:
$\mathbb{f}_{ij}^{interact}=\frac{m_j}{\pi R^4\rho_j}(1-q_{ij})[15k(\rho_i + \rho_j - 2\rho_0)\frac{(1-q_{ij})}{q_{ij}}\mathbb{r}_{ij}-40\mu \mathbb{v}_{ij}]$
Here, $\mathbb{r}_{ij}$ is the vector from i to j, $\rho_i, \rho_j$ are the densities of the particles, $\mathbb{v}_{ij}$ their relative velocity, $q_{ij} = \frac{\lVert \mathbb{r}_{ij} \rVert}{R}$ and $\mu$ is the viscosity.
$k$ and $\rho_0$ are constants of the fluid.

Now, we use a more simplified version, this one: $\mathbb{a}_{ij}^{interact} = \frac{1-q_{ij}}{\rho_j}[P(\rho_i+\rho_j-2\rho_0)\mathbb{r}_{ij} - \mu \mathbb{v}_{ij}]$ , which is scaled later by a constant factor to make the animation realistic.
$P$ and $\mu$ are the pressure and viscosity parameters, and $\rho_0$ is set to $1.5$ .

In code:

// calculate interaction force
for(p=array; p<arr_end; p+=PARTICLE_LEN){
    // gravity component
    p[ACCEL]=G;
    for(q=array; q<arr_end; q+=PARTICLE_LEN){
        diff = p[POS] - q[POS];
        double rij = cabs(diff);
        if(rij < P_RADIUS){
            double complex press = (p[DENSITY]+q[DENSITY]-2*P0)*P*diff;
            double complex visc = (p[VEL] - q[VEL])*V;
            p[ACCEL] +=
                (1-rij/P_RADIUS)/ p[DENSITY] *
                (press - visc);
        }
    }
}

The acceleration, scaled accordingly, is then added to the velocity, which is added to the position.

All of these steps are repeated over and over, simulating the evolution of the fluid one step at a time.

ASCII Fluid Simulations

Analysis of Yusuke Endoh's ASCII fluid dynamics simulator

Making it readable

The algorithm

Input parsing

The simulation algorithm

Output