Forem: Serhii Herasymov

Fractal Monsters Evolution

Serhii Herasymov — Sat, 20 Dec 2025 14:55:50 +0000

Today we will draw geometric fractals, which are given undeservedly little attention. Meanwhile, each fractal here is a small masterpiece that amazes the imagination!

Algorithm

Let's consider one interesting algorithm using the example of a well-known fractal - the Levy curve. We have two points A1 and A2. Let's find a point A3 for which the angle A1 = 45°, and the angle A3 = 90°.

We do the same operation for points A1, A3 and A3, A2.

And then recursively for each pair of points.

Third iteration:

At iteration 14, this curve looks like this:

More clearly in JavaScript:

  <canvas id="myCanvas"></canvas>
  <script>
    window.onload=function(){
      var canvas=document.getElementById('myCanvas');
      var context=canvas.getContext('2d');
      canvas.width=800, canvas.height=800;
      context.beginPath();
      recurs(context, 17, 400, 250, 400, 550); //рекурсивная функция, рисующая фрактальную кривую
      context.stroke();
    };
    var angle=45*Math.PI/180; //переводим углы в радианы
    function recurs(context, n, x0, y0, x1, y1){ //n итераций, Точки A1(x0,y0) и A2(x1,y1)
      if (n==0){
        context.moveTo(x0,y0);
        context.lineTo(x1,y1);
      }else{
        var xx=Math.cos(angle)*((x1-x0)*Math.cos(angle)-(y1-y0)*Math.sin(angle))+x0;
        var yy=Math.cos(angle)*((x1-x0)*Math.sin(angle)+(y1-y0)*Math.cos(angle))+y0; //находим точку A3
        recurs(context, n-1, x0, y0, xx, yy); //A1, A3
        recurs(context, n-1, xx, yy, x1, y1); //A3, A2
      }
    }
  </script>

By slightly modifying the above code, we can draw another well-known fractal curve - the Harter-Haytheway Dragon. To do this, starting from the second iteration, we will alternate angles of 45° and -45°

Third iteration:

At iteration 14, this curve looks like this:

Modified fragment of the source code

      }else{
        angle2=angle*k; // <-
        var xx=Math.cos(angle2)*((x1-x0)*Math.cos(angle2)-(y1-y0)*Math.sin(angle2))+x0;
        var yy=Math.cos(angle2)*((x1-x0)*Math.sin(angle2)+(y1-y0)*Math.cos(angle2))+y0;
        recurs(context, n-1, x0, y0, xx, yy, 1); //последний аргумент функции - k
        recurs(context, n-1, xx, yy, x1, y1, -1);
      }

We have changed our function just a little bit, adding one more angle to it, and the fractal has changed beyond recognition! In dynamics (we change the second angle in 5° increments):

And here our inquisitive mind starts asking questions. What if we try using other angle values, instead of [45,-45]? Say, [30,-60]. Or [30,-30,-60,60]. Let's rewrite our function so that it accepts an array of angles as one of the arguments.

The final version of the function

    function recurs(context, arr, position, n, x0, y0, x1, y1){
      if (n==0){
        context.fillRect(x1,y1, 1,1);
        return position;
      }else{
        if (!position[n]) position[n]=0;
        var xx=Math.cos(arr[position[n]])*((x1-x0)*Math.cos(arr[position[n]])-(y1-y0)*Math.sin(arr[position[n]]))+x0;
        var yy=Math.cos(arr[position[n]])*((x1-x0)*Math.sin(arr[position[n]])+(y1-y0)*Math.cos(arr[position[n]]))+y0;
        position[n]++;
        if (position[n]==arr.length) position[n]=0;
        position=recurs(context, arr, position, n-1, x0, y0, xx, yy);
        position=recurs(context, arr, position, n-1, xx, yy, x1, y1);
        return position;
      }
    }
// The array position stores the angle number from the array arr for the current recursion depth (n).

From now on we will draw only dots, without connecting them with lines... purely for aesthetic reasons. Fractal [30,-30,-60,60] drawn with lines and dots:

Special attention should be paid to the fact that all the dots A3 lie on a circle with a diameter of A1A2:

If the dots are inside (or outside) the circle, the fractal is too compressed (or vice versa - goes beyond the drawing area). Visually:

With this cosine we align the dots along the circle:

var xx=Math.cos(arr[position[n]])* ...
var yy=Math.cos(arr[position[n]])* ...

It is easy to see that using an angle of 120° we get the same point as with an angle of -60°. Therefore, all angles are in the range (-90°,90°).

Let's try to enter different combinations of angles into our script.

"Little Masterpieces"

Several random fractals with angles of -45°/45° (clickable):

With angles of -30°/30°/-60°/60°:

With angles of -15°/15°/-30°/30°/-45°/45°/-60°/60°/-75/75°:

… We can go on like this for a very long time. If we use angles of 15°/30°/45°/60°/75° (positive and negative) and, say, draw “octagonal” fractals — how many fractals can we draw? About 10^8 = 100,000,000 pieces. In addition, the set of “octagonal” fractals does not intersect with the set of, say, “nonagonal” fractals (except for fractals with all the same angles). In general, a lot of fractals can be drawn with this algorithm.

Manually searching for something in this set, without knowing how this “something” should look, is a bit difficult. What prevents us from adding a genetic algorithm to this? Nothing. Let’s add it!

Evolution

A genetic algorithm is a heuristic search algorithm that uses the mechanisms of biological evolution. The algorithm is simple, intuitive and at the same time quite effective.

It is divided into several stages.

An initial population is created, consisting of a small number of individuals. Each individual is an array of fixed length with a set of genes (genotype). The genotypes of the initial population are filled randomly. In relation to our problem, genes are the angles with which we will draw fractals.

The next stage is selection. For each individual, we calculate the fitness coefficient. The better the individual, the higher the coefficient. We sort the population by this coefficient. We divide the population in half - we will form the next population from the most adapted. (How to divide and what to choose for the next population is a topic that can be debated for a long time).

The next stage is crossing. From the most adapted, we select (randomly) two individuals - parents. We make two offspring: we randomly take a part of the genes of one parent and a part of the genes of the other, fill them into the genotype of one of the offspring. The remaining genes go to the second offspring. We take the next two parents. We do the same operation. So until the parents run out. We form a new population from the parents and offspring.

The final stage is mutations. We take a certain percentage of individuals from the new population, for each of them (individuals) we randomly select a gene and replace it with a random value. In this way, we add genes to our gene pool that were not in the initial population and which, in the long run, can form the best results. In addition, by increasing the percentage of mutations, we can to some extent solve the problem of convergence to a local optimum.

Each stage is more clearly in JavaScript.

Initial population

Generating angles

function randomangl(min, max, step){
  if(step==0) step=1;
  var rand=Math.floor(Math.random() * (max/step - min/step + 1)) + min/step;
  if(rand==0) rand=1;
  rand*=step;
  if(rand>max) rand=max;
  return Math.round(rand);
}

The randomangl function is called with three arguments: the minimum angle, the maximum angle, and the step. This is useful if we want to populate the initial population with specific angles. For example, we can call the function with the following arguments:

randomangl(-45, 45, 90) — generates angles of -45°/45°
randomangl(-60, 60, 30) — generates angles of -60°/-30°/30°/60°
randomangl(-75, 75, 15) — generates angles of -75°/-60°/-45°/-30°/-15°/15°/30°/45°/60°/75°

We do not use angles of 0°. The fractal [45°,0°] would look like this:

We get the same fractal [45°], only in the place where we try to draw 0° — the points merge (all further recursion in this place becomes meaningless). For the same reason, there is no point in generating fractals using 90°. For the initial population, it is best to call the randomangl(-75, 75, 15) function. For mutations — randomangl(-89, 89, 1).

Create the initial population

population=[];
fitness=[];
for(var n=anglemin; n<=anglemax; n++){
  population[n]=[];
  fitness[n]=[];
  for(var i=0; i<PopulationSize; i++){
    population[n][i]=[];
    fitness[n][i]=0;
    for(var j=0; j<n; j++){
      population[n][i][j]=randomangl(min, max, step);
      if(j==0) population[n][i][j]=Math.abs(population[n][i][j]);
    }
  }
}

Create several populations at once with different numbers of angles (anglemin, anglemax). Initialize a separate array with fitness coefficients. The first gene of each individual is forced to be positive — the fractals [45°,-45°] and [-45°,45°] are symmetrical.

Selection

In general, the fitness of each individual in classical genetic algorithms is determined by a special fitness function. The function evaluates how well each individual solves the task. In our case, it is not possible to clearly formulate the task (we do not know how this “something” should look), so we will attach it to the user algorithm as a fitness function.

There are two options. We could draw all the fractals from the populations and ask the user to select the best ones, but some difficulties may arise - how to choose the best fractals when they all look impressive? It is much easier to compare two random fractals.

Interface:

There are two canvases on the page: myCanvas1 and myCanvas2. And two buttons: onclick="selecter(1);" and onclick="selecter(2);"

We call the function that draws two random fractals:

Draw two random fractals

function get2fractals(){
  PopulationNumber=Math.floor(Math.random() * (anglemax - anglemin + 1))+anglemin;
  drawcanvas(1);
  drawcanvas(2);
}

function drawcanvas(n){
  var canvas=document.getElementById('myCanvas'+n);
  var context=canvas.getContext('2d');
  canvas.width=canvasSize, canvas.height=canvasSize;
  context.fillStyle = 'rgb(255,255,255)';
  context.fillRect (0, 0, canvas.width, canvas.height);
  var position=[];

  if(n==1){
    FractalNumber[1]=Math.floor(Math.random() * (PopulationSize));
  }else{
    do{
      var rnd = Math.floor(Math.random() * (PopulationSize));
    }while(rnd==FractalNumber[1])
    FractalNumber[2]=rnd;
  }

  var array1fractal=population[PopulationNumber][FractalNumber[n]];

  context.fillStyle = 'rgb(0,0,0)';
  draw(context, array1fractal, position, iteration, xA, yA, xB, yB);
}

In the get2fractals function, we select a random population.
In the drawcanvas function, we check if(n==1) — to avoid drawing two identical fractals. We store the population number and the numbers of the two fractals in global variables (PopulationNumber and FractalNumber[]).

When the user clicks the Select button next to the fractal they like, we call the following function:

Select

function selecter(n){
    fitness[PopulationNumber][FractalNumber[n]]++;
    get2fractals();
}

Which increases the fitness coefficient for the selected fractal and draws two new fractals.

We will assume that the user was conscientious and made no fewer choices than the number of individuals in the population before launching the evolution mechanism.

Crossing and mutations

function sortf(a, b) {
  if (a[1] < b[1]) return 1;
  else if (a[1] > b[1]) return -1;
  else return 0;
}

function evolution(){
  var mutation=document.getElementById("mutatepercent").value;
  var exchange=document.getElementById("exchangepercent").checked;

  var sizehalf=PopulationSize/2;
  var sizequarter=sizehalf/2;
  for(var n=anglemin; n<=anglemax; n++){

    var arrayt=[];
    for(var i=0; i<PopulationSize; i++){
      arrayt[i]=[];
      arrayt[i][0]=population[n][i];
      arrayt[i][1]=fitness[n][i];
    }
    arrayt.sort(sortf);
    arrayt.length=sizehalf;
    population[n].length=0;
    fitness[n].length=0;

    for(var i=0; i<sizequarter; i++){
      var i0=i*4;
      var i1=i*4+1;
      var i2=i*4+2;
      var i3=i*4+3;

      var removed1=Math.floor(Math.random()*(arrayt.length));
      var parent1f = arrayt.splice(removed1,1);
      var parent1=parent1f[0][0];
      var removed2=Math.floor(Math.random()*(arrayt.length));
      var parent2f = arrayt.splice(removed2,1);
      var parent2=parent2f[0][0];

      var child1=[];
      var child2=[];

      for(var j=0; j<n; j++){
        var gen=Math.round(Math.random());
        if(gen==1){
          child1[j]=parent1[j];
          child2[j]=parent2[j];
        }else{
          child1[j]=parent2[j];
          child2[j]=parent1[j];
        }
      }

      population[n][i0]=parent1;
      population[n][i1]=parent2;
      population[n][i2]=child1;
      population[n][i3]=child2;

      fitness[n][i0]=0;
      fitness[n][i1]=0;
      fitness[n][i2]=0;
      fitness[n][i3]=0;
    }

    if(mutation>0){
      var m=100/mutation;
      for(var i=0; i<PopulationSize; i++){
        var rnd=Math.floor(Math.random()*(m))+1;
        if(rnd==1){
          var mutatedgene=Math.floor(Math.random()*(n));
          population[n][i][mutatedgene]=randomangl(minm, maxm, stepm);
        }
      }
    }
  }

  if(exchange){
    ///...
  }
  get2fractals();
}

We take each population in turn. We write all individuals to a temporary array, and we also write their fitness coefficients there. We sort the temporary array by these coefficients. We cut off half (the fittest). We reset the old two arrays so that nothing extra is left there. We pull two ancestors from the temporary array, form two descendants. Fill the population with ancestors and descendants.

We make mutations. A small note. Instead of the percentage of mutations in the entire population, we use the probability of mutation of a single individual.

The exchange option is available. This option implements a model of migration of individuals between populations.

Migration model

if(exchange){
  var n=anglemin;
  while(n<population.length){
    var rnd=Math.round(Math.random());
    var np=n+1;
    if(rnd==1 && (typeof population[np]!="undefined")){

      var i1=Math.floor(Math.random() * (PopulationSize));
      var i2=Math.floor(Math.random() * (PopulationSize));
      var tempa1=population[np][i1];
      var tempa2=population[n][i2];
      var indexofexcessgene=Math.floor(Math.random()*(tempa1.length));
      var excessgene=tempa1.splice(indexofexcessgene,1);
      tempa2.splice(indexofexcessgene,0,excessgene[0]);
      population[np][i1]=tempa2;
      population[n][i2]=tempa1;

      n+=2;
    }else{
      n++;
    }
  }
}

If the Exchanges checkbox is set, one random individual with a 50% probability can exchange genes with a random individual from the next population (population n + 1) in turn for each population. There are more genes in the genotypes of the next population, we perform the exchange as follows. We select a random gene from an individual in population n + 1, insert it into the same position in the genotype of an individual from population n. We swap individuals:

It remains to add bells and whistles (css buttons, saving the entire process in localStorage, displaying the Parents Tree ...).

You can run it in a browser by going to the project website

The result of my selections is in the picture at the beginning of the publication.

And finally, I will add that 2D is good, but 3D is even better!

You can rotate the fractal in three dimensions with your mouse here. After a slash, enter the angle pairs /alpha/beta/alpa/beta/… in the address bar.

Alpha is the rotation angle we are already familiar with. The new angle Beta is the rotation of point A3 around the axis A1A2.

Bug in 3D algorithm

There is a funny bug (feature) in the 3D algorithm, which was not fixed. To draw the Harter-Hayteway Dragon, you need to use 4 pairs of angles: 45/180/45/0/45/0/45/180.
The 3D algorithm can draw two-dimensional fractals that cannot be drawn with the 2D algorithm. For example:

45/0/45/180/45/180/45/180/45/0

45/180/45/180/45/0/45/180/45/180/

A small gallery of two-dimensional fractals drawn by this algorithm.

Repository

Author

Serhii Herasymov

sergeygerasimofff@gmail.com

https://github.com/xcontcom

Perfect Shuffle: The Strange Fractals Hidden in a Simple Algorithm

Serhii Herasymov — Wed, 10 Dec 2025 18:57:08 +0000

The Perfect Shuffle

Repository

I have always been fascinated by elementary algorithms that create complex patterns. There is something fundamentally intriguing about such algorithms. One such algorithm is the Perfect Shuffle. Let's explore its unusual properties and create some striking fractals using this algorithm.

This article includes numerous images, GIF animations, and a bit of bullshit (12%).

The Perfect Shuffle, known among card magicians as the Faro Shuffle, captivated me when I began learning card tricks. Mastery starts with sleight of hand. After a week of practicing with a deck of cards, I became proficient in basic shuffling techniques: Riffle Shuffle, Faro Shuffle, and Charlier Pass. While practicing the Faro Shuffle, I arranged the deck with all black suits at the beginning and red suits at the end. After several shuffles, I noticed an intriguing order in the cards. I set the deck aside and turned to programming.

The Perfect Shuffle algorithm reorders elements in an array. It is remarkably simple:

Divide the array into two equal parts.
Place elements from the first half into even positions and elements from the second half into odd positions.

Visually:

Here’s a JavaScript implementation:

function shuffle(array) {
    var half = array.length / 2;
    var temparray = [];
    for (var i = 0; i < half; i++) {
        temparray[i * 2] = array[i + half];
        temparray[i * 2 + 1] = array[i];
    }
    return temparray;
}

var array = [1, 2, 3, 4, 5, 6, 7, 8];
array = shuffle(array);
console.log(array); // [5, 1, 6, 2, 7, 3, 8, 4]

A Brief Digression

Card enthusiasts call this version the "in-shuffle." There is also an "out-shuffle," where elements from the first half are placed in odd positions:

In the out-shuffle, outer elements remain in place, while the middle is shuffled like the in-shuffle. For simplicity, we will focus on the in-shuffle and assume arrays have an even number of elements, as odd-length arrays leave the last element static.

The Perfect Shuffle has fascinating properties. After a certain number of iterations, the array returns to its original order. This is because the algorithm is deterministic: each state has a unique predecessor, and the number of possible states is finite. The number of iterations required to restore the original order is less than or equal to the array length.

For example:

An array of 12 elements requires 12 iterations:

An array of 14 elements needs only 4 iterations:

The following image (clickable) shows arrays with lengths from 2 to 20 elements:

(unshuffle)

These arrays return to their original state after the following number of iterations:
2, 4, 3, 6, 10, 12, 4, 8, 18, 6
This sequence corresponds to A002326.

function a002326(n){
    var a=1;
    var m=0;
    do{
        a*=2;
        a%=2*n+1;
        m++;        
    }while(a>1);
    return m;
}
for(var n=1;n<11;n++){
    m=a002326(n);
    console.log(m);
}

Let’s visualize this with a graph. The X-axis represents the number of elements in the array (divided by 2), and the Y-axis shows the number of iterations. White dots mark the initial states (origin at the upper left):

The dots appear chaotically distributed.

Another intriguing property is that after a specific number of iterations, the array’s order may reverse. For a 12-element array, this occurs on the 6th iteration. Here’s a graph marking arrays that reverse:

The distribution remains chaotic.

Observing element movement is captivating. Consider filling the first half of the array with zeros and the second with ones (zeros as black pixels, ones as white). Each row (X) represents the array’s state, with iterations along the Y-axis.

For arrays with 142–150 elements:

For 288 elements (144*2):

For 610 elements:

Exploring Other Array Lengths

A small script (shuffle1d.html in the repository) generates these patterns. Enter the array size (half the total length) and click "Draw."

The results remain chaotic. Let’s track a single element’s trajectory in a 610-element array, marking only the first element:

The position of an element after each iteration can be calculated without shuffling the entire array:

x[n+1] = {
    x[n]/2 + y/2, if x[n] % 2 == 0;
    (x[n] - 1)/2, if x[n] % 2 == 1;
}

(where y is the array length)

For arrays of lengths 2 to 22, the upper part of the image shows the array and the element’s position per iteration, while the lower part shows all possible positions:

The trajectory is unpredictable. Stacking these trajectories produces:

(Clickable)

Black pixels indicate positions the first element never occupies during shuffling.

Let’s extend this to matrices. First, two array-based tricks:

Reverse the order of elements in the second half at each iteration:

Invert the elements in the second half:

These operations will be useful later.

Now, matrices. Fill a matrix with a diagonal pattern (first row: first element, second row: second element, etc.):

Shuffle the columns using Perfect Shuffle to track element movement. For an 80x80 matrix:

Adding another diagonal reveals a moiré pattern:

For a 242x242 matrix, iterations 27 to 35:

Now, shuffle both rows and columns. Here’s the modified function:

function shufflecr(array) {
    var half = array.length / 2;
    var temparray = [];

    for (var x = 0; x < half; x++) {
        temparray[x * 2] = [];
        temparray[x * 2 + 1] = [];
        for (var y = 0; y < half; y++) {
            temparray[x * 2][y * 2] = array[x + half][y + half];
            temparray[x * 2 + 1][y * 2] = array[x][y + half];
            temparray[x * 2][y * 2 + 1] = array[x + half][y];
            temparray[x * 2 + 1][y * 2 + 1] = array[x][y];
        }
    }
    return temparray;
}

This function shuffles both rows and columns. Visually:

Shuffling diagonals is uneventful:

Shifting the diagonals one pixel right yields:

The function works, but the result is lackluster. Let’s try a random square in an 80x80 matrix (chosen because the order reverses on the 27th iteration):

Animate the square’s movement:

Now, draw a circle:

Notice the edge effects: the matrix behaves like a torus, especially visible on the second iteration:

Trigonometric Connections

Consider the 288-element array pattern:

Write it into a matrix and shuffle:

Or, take the original pattern:

Keep pixels at coordinates x*4+i, y*4+j (for i=2, j=2):

Remove empty rows and columns:

For other i and j, this resembles a perfect unshuffle, dividing the matrix into 16 parts:

These patterns resemble those generated by trigonometric functions like z = sin(n*x/y) or z = sin(n*x*y):

For z = sin(3*x*y) (white if z >= 0, black if z < 0):

For z = sin(13*x*y):

Another pattern:

With a fill to highlight closed areas:

Reversing the order of elements produces striking patterns:

With a fill:

Inversion creates intriguing results. Start with an empty matrix, divide it into four parts, invert two, and shuffle:

The first iteration is unremarkable, but subsequent iterations reveal patterns:

To understand this, return to arrays. Instead of fixed lengths, grow the array iteratively:

Create an array of n elements.
Make a copy and invert it.
Shuffle the copy with the original to get an array of n * 2 elements.
Repeat from step 2.

Starting with [1, 0]:

[1, 0] (original)
[1, 0, 0, 1] (array and inverted copy)
[0, 1, 1, 0] (shuffled)
[0, 1, 1, 0, 1, 0, 0, 1] (array and inverted copy)
[1, 0, 0, 1, 0, 1, 1, 0]
[0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0]
And so on.

This generates the Morse-Thue sequence (A010060), a simple fractal.

Back to matrices. Instead of shuffling four parts, copy the matrix, apply transformations, and shuffle with the copies. Possible transformations include:

Action 0: Unchanged matrix
Actions 1–7: Rotations
Actions 8–15: Inversions

Visually:

Apply actions 0, 0, 7, 7 (action 7 rotates by 90°):

Next iteration:

Eighth iteration:

This produces a unique fractal.

Alternative Fractal Method

Without Perfect Shuffle, apply actions 0, 0, 7, 7. First iteration:

Second iteration:

Eighth iteration:

Two additional tricks:

Overlay a copy with transparent white pixels, rotated 90°:

Rotate the copy vertically:

By combining these actions, we can generate 16^4 = 65,536 fractals.


2.1.15.14 grey	2.1.15.14 rgb	2.1.15.14 binary

6.14.0.0 grey	6.14.0.0 rgb	6.14.0.0 binary

9.4.0.11 grey	9.4.0.11 rgb	9.4.0.11 binary

4.8.4.7 grey	4.8.4.7 rgb	4.8.4.7 binary

13.14.0.3 grey	13.14.0.3 rgb	13.14.0.3 binary

14.1.11.10 grey	14.1.11.10 rgb	14.1.11.10 binary

4.13.15.0 grey	4.13.15.0 rgb	4.13.15.0 binary

15.6.8.3 grey	15.6.8.3 rgb	15.6.8.3 binary

4.0.2.14 grey	4.0.2.14 rgb	4.0.2.14 binary

5.0.2.15 grey	5.0.2.15 rgb	5.0.2.15 binary

2.0.2.9 grey	2.0.2.9 rgb	2.0.2.9 binary

14.3.1.0 grey	14.3.1.0 rgb	14.3.1.0 binary

4.1.11.8 grey	4.1.11.8 rgb	4.1.11.8 binary

10.0.3.0 grey	10.0.3.0 rgb	10.0.3.0 binary

Explore more with script (hybrid.html in the repository).

Here are the "card tricks" we got. I think we can stop here for now (you can shuffle three-dimensional matrices, but that's another time), (if any of the grasshoppers want to look at shuffling three-dimensional matrices - script (fractal_3d_2.html in the repository)).

Melody Generation

I experimented with a melody generator using a genetic algorithm. Each melody consists of note genes. Random shuffling produces cacophony, but Perfect Shuffle creates intriguing melodies.

Bonus: A Connection Between Perfect Shuffle and Billiard Fractals

If you enjoyed this exploration of shuffling arrays, matrices, and fractals, there’s a surprising connection waiting around the corner.

The Perfect Shuffle sequence turns out to be deeply related to the patterns generated by discrete billiard reflections inside a rectangle — two systems that appear completely unrelated, yet produce identical structures when you look at them symbolically.

Curious?

Read the companion article:
Billiard Fractals — The Infinite Patterns Hidden in a Rectangle

Final Note

"No cats were shuffled in the making of this fractal. Any resemblance to space chicks is purely coincidental."

About the Author

Serhii Herasymov

Programmer & Math Enthusiast (Ukraine → Ireland)

Email: sergeygerasimofff@gmail.com

GitHub: xcontcom

X (Twitter): @xcontcom

Feedback welcome.

Evolving Cellular Automata

Serhii Herasymov — Sun, 07 Dec 2025 13:32:37 +0000

Evolving Cellular Automata

Repository

Introduction

Hello, curious minds! Let’s explore the fascinating intersection of cellular automata and genetic algorithms to uncover emergent patterns and behaviors.

This article is a translation and adaptation of my 2019 piece, originally published in Russian on Habr. Here, we’ll dive into the mechanics of cellular automata, evolve them using genetic algorithms, and showcase some intriguing results.

Cellular Automata Rules

The simplest form of cellular automata is the one-dimensional variant. (While zero-dimensional oscillators exist, we’ll set them aside for now.) In a one-dimensional cellular automaton, we start with a single array representing the initial state, where each cell holds a binary value (0 or 1). The next state of each cell depends on its current state and those of its two immediate neighbors, determined by a predefined rule.

With three cells (the cell itself and its two neighbors), there are 2^3 = 8 possible configurations:

For each configuration, we define the cell’s next state (0 or 1), forming an 8-bit rule, known as the Wolfram code. This results in 2^8 = 256 possible one-dimensional cellular automata.

Manually inspecting all 256 rules is feasible, but let’s scale up to two-dimensional cellular automata, where things get exponentially more complex.

In a two-dimensional automaton, we use a grid (matrix) instead of an array. Each cell has eight neighbors in a Moore neighborhood, including horizontal, vertical, and diagonal neighbors. (The Von Neumann neighborhood, excluding diagonals, is less commonly used.)

For consistency, we order the neighbors as follows:

With a cell and its eight neighbors, there are 2⁹ = 512 possible configurations, and the rules are encoded as a 512-bit string. This yields 2⁵¹² ≈ 1.34 × 10¹⁵⁴ possible two-dimensional automata — a number far exceeding the estimated ~10⁸⁰ atoms in the observable universe.

Exploring all these rules manually is impossible. If you checked one rule per second since the Big Bang, you’d have only covered (4.35 \times 10^{17}) by now. Enter genetic algorithms, which allow us to search for automata that meet specific criteria efficiently.

Two-Dimensional Cellular Automata

Let’s build a two-dimensional cellular automaton. We’ll start by generating a random 512-bit rule:

const rulesize = 512;
const rule = new Array(rulesize).fill(0).map(() => Math.round(Math.random()));

Next, we initialize an 89x89 grid (chosen for its odd dimensions to maintain dynamic behavior and computational efficiency):

const sizex = 89;
const sizey = 89;
const size = 2;
const a = [];

for (let x = 0; x < sizex; x++) {
    a[x] = [];
    for (let y = 0; y < sizey; y++) {
        a[x][y] = Math.round(Math.random());
        if (a[x][y]) context.fillRect(x * size, y * size, size, size);
    }
}

To compute the next state, we use periodic boundary conditions (wrapping the grid into a torus shape to avoid edge effects):

function countpoints() {
    const temp = new Array(sizex);
    for (let x = 0; x < sizex; x++) {
        temp[x] = new Array(sizey);
        const xm = (x - 1 + sizex) % sizex; // Wrap around using modulo
        const xp = (x + 1) % sizex;

        for (let y = 0; y < sizey; y++) {
            const ym = (y - 1 + sizey) % sizey;
            const yp = (y + 1) % sizey;

            // Calculate the bitmask in a single step
            let q = (
                (a[xm][ym] << 8) |
                (a[x][ym] << 7) |
                (a[xp][ym] << 6) |
                (a[xm][y] << 5) |
                (a[x][y] << 4) |
                (a[xp][y] << 3) |
                (a[xm][yp] << 2) |
                (a[x][yp] << 1) |
                a[xp][yp]
            );

            temp[x][y] = rule[q];
            if (temp[x][y]) context.fillRect(x * size, y * size, size, size);
        }
    }
    a = temp;
}

Running this automaton with a random rule often produces chaotic, white-noise-like behavior, akin to static on an untuned TV:

Let’s “tune” this TV using a genetic algorithm to find more structured patterns.

Genetic Algorithm

Think of the genetic algorithm as tuning a TV to find a clear signal. Instead of manually adjusting knobs, we define desired properties, and the algorithm searches for rules that best match them. The vast number of possible automata ensures diverse outcomes with each run.

We’ll use a population of 200 automata, each defined by a 512-bit rule, stored in a population[200][512] array. We also track fitness scores for each individual:

const PopulationSize = 200;
const rulesize = 512;
const population = [];
const fitness = [];

for (let n = 0; n < PopulationSize; n++) {
    population[n] = new Array(rulesize).fill(0).map(() => Math.round(Math.random()));
    fitness[n] = 0;
}

The genetic algorithm involves two key processes: selection and evolution.

Selection

We evaluate each automaton’s fitness based on how well it meets our criteria. The top 50% (100 individuals) with the highest fitness survive, while the rest are discarded.

Evolution

Evolution consists of crossover and mutation:

Crossover: We pair surviving individuals to create offspring. For each pair, we select a random crossover point and swap genes to produce two new descendants:

Mutation: We apply a 5% mutation rate, where each gene in an individual’s rule has a 5% chance of flipping (0 to 1 or vice versa). Higher mutation rates can introduce beneficial diversity but risk destabilizing successful traits.

Here’s the evolution function:

function evolute() {
    const sizehalf = PopulationSize / 2;
    const sizequarter = sizehalf / 2;
    const arrayt = population.map((p, i) => [p, fitness[i]]);
    arrayt.sort((a, b) => b[1] - a[1]); // Sort by fitness, descending
    arrayt.length = sizehalf; // Keep top half
    const newPopulation = [];
    const newFitness = [];

    for (let i = 0; i < sizequarter; i++) {
        const i0 = i * 4;
        const i1 = i * 4 + 1;
        const i2 = i * 4 + 2;
        const i3 = i * 4 + 3;

        const removed1 = Math.floor(Math.random() * arrayt.length);
        const parent1f = arrayt.splice(removed1, 1)[0];
        const parent1 = parent1f[0];
        const removed2 = Math.floor(Math.random() * arrayt.length);
        const parent2f = arrayt.splice(removed2, 1)[0];
        const parent2 = parent2f[0];

        const qen = Math.floor(Math.random() * rulesize);
        const child1 = parent1.slice(0, qen).concat(parent2.slice(qen));
        const child2 = parent2.slice(0, qen).concat(parent1.slice(qen));

        newPopulation[i0] = parent1;
        newPopulation[i1] = parent2;
        newPopulation[i2] = child1;
        newPopulation[i3] = child2;

        newFitness[i0] = 0;
        newFitness[i1] = 0;
        newFitness[i2] = 0;
        newFitness[i3] = 0;
    }

    const mutation = document.getElementById("mutatepercent").value * 1;
    const m = 100 / mutation;
    for (let i = 0; i < PopulationSize; i++) {
        if (Math.random() * m < 1) {
            const gen = Math.floor(Math.random() * rulesize);
            newPopulation[i][gen] = newPopulation[i][gen] ? 0 : 1;
        }
    }

    population = newPopulation;
    fitness = newFitness;
}

Note: The original article used a single crossover point, but experiments suggest that multiple random crossover points or higher mutation rates (e.g., 25%) can improve evolution speed and quality. However, we’ll stick with the 5% mutation rate and single-point crossover for consistency with the original experiments.

Experiments

To guide natural selection, we define fitness criteria for desirable automata and let the genetic algorithm optimize for them.

Experiment 1: Static Patterns

Chaotic automata flicker like static. Let’s seek stable patterns by comparing the 99th and 100th states and counting unchanged cells as fitness. To avoid trivial solutions (e.g., all 0s or all 1s), we add a constraint: the difference between the number of 0s and 1s must be less than 100.

function countfitness(array1, array2) {
    let sum = 0;
    let a0 = 0;
    let a1 = 0;
    for (let x = 0; x < sizex; x++) {
        for (let y = 0; y < sizey; y++) {
            if (array1[x][y] === array2[x][y]) sum++;
            if (array1[x][y] === 0) a0++;
            else a1++;
        }
    }
    return Math.abs(a0 - a1) < 100 ? sum : 0;
}

After 421 epochs, we found an optimal solution where all (89 \times 89 = 7921) cells remain unchanged:

Blue dots: Best individuals per epoch.
Red dots: Worst individuals meeting the second criterion.
Black dots: Average fitness (including those failing the second criterion).

Gene pool (Y: individuals, X: genes):

The resulting automaton:

Running the algorithm again yields different stable patterns:

Experiment 2: Pattern Matching

Let’s search for automata that produce specific patterns in a second-order Moore neighborhood (25 cells). We start with this pattern:

Since each cell’s state depends on only 9 neighbors, matching a 25-cell pattern is challenging. A random automaton at the 100th iteration shows no matches:

To make the search feasible, we:

Allow one mistake in pattern matching.
Check the last 50 states (iterations 51–100).
Drop the 0s vs. 1s constraint.

After 3569 epochs:

Y-axis: Arbitrary fitness scale (~30 patterns found, compared to 7921 in the static case).
X-axis: Epochs, with white lines at 500 and 1000.
Blue/red/black dots: Best, worst, and average fitness.

Solutions at 500, 1000, and 3569 epochs:

Gene pool at 3569 epochs:

Dynamic behavior:

A single cell evolving into an oscillator:

Marking genes that form the oscillator (grey lines indicate unused genes):

Adding a constraint (difference between 0s and 1s ≤ 400) yields:

First 100 iterations:

Final iteration:

Tightening the constraint to ≤ 100 after 14865 epochs:

Scaled:

Dynamic behavior (50 iterations):

50th iteration:

Experiment 3: Another Pattern

Targeting this pattern with exact matching:

After 3000 epochs:

Dynamic behavior (100 iterations):

100th iteration:

Experiment 4: Precise Pattern Matching

For this pattern, we require exact Matches:

After 4549 epochs:

The solution at 4549 epochs, shown dynamically over 100 iterations:

After 500–2000 iterations, the automaton’s state becomes nearly fully ordered. The 89x89 grid size (chosen as an odd number) prevents complete ordering, maintaining some dynamic behavior:

However, with a 90x90 grid, the same automaton achieves full ordering:

Let’s examine the intermediate solution at 2500 epochs:

This automaton transforms a chaotic initial state into an ordered final state, characterized by a pattern shift left along the X-axis plus a few oscillator cells.

The ordering speed varies by grid size:

A 16x16 automaton orders in approximately 100 iterations:

A 32x32 automaton orders in about 1000 iterations:

A 64x64 automaton orders in roughly 6000 iterations:

A 90x90 automaton orders in about 370,000 iterations:

An 11x11 (even-sized) automaton orders in approximately 178,700 iterations:

A 13x13 automaton did not order within a reasonable timeframe.

On a 256x256 grid at the 100,000th iteration, the automaton displays remarkable self-organization:

Observing this self-organization process on a large field is captivating:

Interactive visualization

Rerunning the evolution yields different solutions, such as this one:

Experiment 5: Next Pattern

Targeting this pattern, allowing one mistake in matching:

After 5788 epochs, the fitness graph (arbitrary scale):

Dynamic behavior of the automaton:

The ordered state consists of a pattern shifted upward along the Y-axis with a few oscillator cells:

Mutants in this population are particularly intriguing. Here’s one on a 256x256 grid over 200 iterations:

Interactive visualization

Experiment 6: Complex Pattern (“habr”)

Original article was written on habr. For habr site I looked for specific pattern - letters "habr":

This pattern is a bit irrelevant here, but I made some experiments with this pattern and they quite interesting.

In previous experiments, we calculated the sum of cells around which the pattern is formed (if there was one error, 1 was added to the sum; if there were no errors, 2 were added). The resulting sum was used as a fitness coefficient for the genetic algorithm.

This method does not work for a more complex pattern. An automaton in which a smaller number of cells more accurately correspond to a given criterion (the number of cells that match the pattern in the vicinity of a cell) will lose every time to an automaton in which a larger number of cells less accurately correspond to a given criterion. As in the example with squares above:

For the desired pattern, on the 100th iteration of each automaton in the population, in the environment of each cell we will count the number of cells that match the pattern. We will take only the best result for each automaton. The number of cells that match the pattern will be used as the fitness coefficient. The pattern consists of 7x17=119 cells. This number will be considered the optimal solution. 6000 evolution cycles allowed us to find an automaton that draws a pattern with 5 errors (114 cells match the pattern).

Graph:

With 25% mutations, the gene pool is diverse:

Dynamic behavior:

Best solution at 100th iteration:

Original vs. found pattern:

Here another graph, that shows how the system evolves with different percent of mutations:

Red dots is average fitness, black dots - all individuals. Even with 100% mutants system evolves. But 25% gives us the best evolution process.

Let's make another experiment.

The same pattern as before. 5% of mutations, 1-8 genes mutate (random amount between 1 and 8). 100 epochs of evolution:

The graph is a heat map. The size of a point on the graph is 5 pixels. The origin (0,0) is the upper left corner.
The Y axis (from 0 to 100) is the evolution cycles. The X axis (from 0 to 119) is the number of cells matching the pattern (for each individual in the population, we take the best result). The brightness of the point is the number of individuals with the specified (X coordinate) result.

Let's run the genetic algorithm 4 times with the same parameters (100 cycles, 5% mutations, up to 8 genes mutate). The graph shows all 5 runs:

Experiment 7: Improved Crossover

The original single-point crossover:

A more effective multi-point crossover:

Results with 5% and 25% mutations:

Using genetic algorithms to evolve cellular automata reveals a vast landscape of patterns and behaviors. From static grids to dynamic oscillators and complex shapes like “habr,” the interplay of selection, crossover, and mutation uncovers solutions that would be impossible to find manually. Experimenting with mutation rates and crossover methods highlights the robustness and flexibility of this approach, offering endless possibilities for discovery.

Try running these experiments yourself to find new patterns — each run is a new adventure!

Repository

What’s Next: Artificial Selection (Part 2)

That’s all for natural selection. For more experiments with artificial selection and second-order cellular automata:

Full article (Part 2)

Conclusion

In this project we used a genetic algorithm to explore the enormous space of two-dimensional cellular automata.

Already for first-order, two-state, two-dimensional automata there are 2^512 possible rules, and for second-order automata there are 2^1024.

This is an astronomically large space: far more automata than there are atoms in the observable universe.

The famous Conway Game of Life is just one particular point in this space, and in fact belongs to a tiny, very special subclass of rules.

When people hear “cellular automata”, they almost always think only of Conway’s automaton, without realizing how many other rules are possible and how rich their behaviour can be.

The experiments in this article illustrate that even simple genetic selection, based on informal visual criteria (“less chaotic”, “darker”, “slower growth”, “more interesting oscillators”), is enough to discover a wide variety of non-trivial automata.

In this work we used the rule of the automaton itself as the genotype and tried to evolve rules that match some chosen criterion.

In the next project we will fix the rule (for example, Conway’s Game of Life or any other chosen rule) and instead use the initial configuration of the automaton as the genotype:

https://github.com/xcontcom/initial-state-evolution

On a common toroidal grid we place two automata, A and B.

The field is wrapped around horizontally for each automaton, and vertically the upper and lower boundaries of the two automata are shared—so each automaton interacts with the other along a common border.

As a fitness function for automaton A we use the state of automaton B at iteration n, and as the fitness function for B we use the state of A at iteration n.

In other words, each automaton’s goal is to influence the configuration of its opponent (or partner).

Conway’s Game of Life is known to be Turing-complete; therefore, in principle, any algorithm or system of arbitrary complexity can be implemented within it.

By co-evolving automata that compete (or cooperate) with each other, we hope to obtain increasingly complex patterns and behaviours.

Unlike the current project, we will not prescribe a fixed target or end state.

Instead, the “goal” of the system is the ongoing competition (or collaboration) between automata.

This makes the evolutionary process more open-ended and may lead to structures that are unexpected and qualitatively different from what we would obtain by optimizing a manually chosen criterion.

Billiard Fractals: The Infinite Patterns Hidden in a Rectangle

Serhii Herasymov — Fri, 28 Nov 2025 03:05:30 +0000

Complex systems often appear chaotic or incomprehensible, yet closer examination reveals that such complexity can frequently be reduced to a simple underlying mechanism. By systematically removing layers of emergent behavior, one can uncover a fundamental rule or equation from which the entire system originates.

While the system described in this article may appear trivial at first glance, the resulting patterns exhibit quasi-fractal behavior that can be analyzed, encoded, and even predicted through symbolic methods. The work presented here was developed independently through direct observation, rather than derived from prior literature.

A useful way to motivate this exploration is by analogy with a common physical phenomenon - wave interference. Consider waves on the surface of a river: a wavefront moves toward the shore, reflects, and overlaps with itself. Do these reflections contain an underlying order? Is it possible to extract structure from the interference?

To investigate this, we simplify the system. Rather than modeling the full wave, we consider only the motion vector - essentially, a ray. We also smooth the “shoreline” and discretize the environment into a rectangular grid. From this setup emerges the core construction of this article.

":)"

The example of waves on the surface of a river serves as a real, intuitive starting point - an accessible physical system that demonstrates how simple rules, such as reflection and interference, can produce complex behavior. It illustrates the central idea: that what appears chaotic often emerges from deterministic structure.

The initial motivation was driven by the conviction that apparent disorder is not randomness, but the result of unresolved or hidden structure. Any system that seems chaotic is governed by rules - its complexity a consequence of perspective, not unpredictability.

To explore this further, attention turned to constructing the simplest possible system that could look chaotic yet remain fully deterministic.

One such system involved a sine wave originating from the corner of a rectangle and reflecting off its boundaries. The nonlinearity of the sine function causes it to intersect itself in complex and unintuitive ways. However, due to limited tools available at the time, the model was simplified even further.

Instead of a sine wave, a straight line was used. The line was made periodic (dashed), and the system was designed to be reproducible using only a pencil and a sheet of graph paper. Despite its simplicity, this construction revealed intricate and structured patterns-forming the foundation of what would later be described as the “billiard fractals.”

Visualizing the Billiard Algorithm

The following sequence illustrates the core mechanism of the discrete billiard system:

An animated version:

Output pattern:

A selection of patterns generated from rectangles with various dimensions:

JavaScript implementation of this algorithm

pattern.js - source code

Fibonacci Numbers and Pattern Refinement

The patterns generated by this system exhibit fractal structure: they are self-similar across scales, recursive in construction, and can be compressed symbolically. As the rectangle dimensions increase following the Fibonacci sequence, the patterns reveal increasingly detailed versions of the same underlying structure.

This refinement process is simple: given rectangle dimensions (m, n), new dimensions are generated by Fibonacci summation. For example, starting with 8×13:

13 + 8 → 21 → becomes 13×21
Then 21 + 13 → 34 → becomes 21×34, and so on

Each step increases resolution while preserving the underlying structure.

8×13	13×21	21×34	34×55	55×89

"233×377 preview comparison"

The article’s header image corresponds to the 233×377 pattern. Its structure can be directly compared with the earlier 13×21 case.

When constructing these patterns using Fibonacci-based dimensions, we are effectively approximating a rectangle with side lengths in the golden ratio - that is, a ratio approaching (1 : φ). With each step, the approximation improves, and the pattern gains additional structure and resolution.

Although the overall structure of the pattern remains consistent during Fibonacci-based refinement, certain symmetries within the pattern depend on the parity of the rectangle's side lengths. Specifically, when both the width and height are odd integers, the resulting pattern exhibits clear diagonal, horizontal, and vertical symmetry. This occurs because the billiard path, under these conditions, terminates in the corner diagonally opposite from its starting point. In contrast, when one or both sides are even, the path terminates elsewhere, and the resulting pattern loses this precise symmetry - although the underlying recursive structure remains unchanged.

Boundary Analysis and Recursive Symmetry

To understand why the structure persists under Fibonacci expansion, consider cutting a square from the pattern. The boundary behavior reveals that the path (i.e., the “billiard ball”) always returns to its entry point:

Moreover, the path always (except for diagonal cases) crosses an even number of cells. This ensures that the pattern remains consistent across such subdivisions.

By recursively separating square regions from the larger pattern, the symbolic seed of the system can be exposed:

Binary Representation and Symbolic Extraction

The path traced by the billiard ball through the grid can be encoded as a binary sequence. As the ball moves from cell to cell, its internal state alternates according to a fixed rule. We can label these alternating states with binary values - for example, assigning 0 to one state and 1 to the other. This produces a binary field that can be visualized directly.

For example:

The top row of the binary field can be viewed as a symbolic boundary - a compact representation of the billiard system's behavior along a single edge. By studying the structure of the full 2D pattern and recursively extracting square sections from it, we arrive at symbolic generation rules. These rules allow us to reconstruct the boundary sequences using only binary operations.

Two core recursive generators are presented below:

function invers(array) {
    var temp = [];
    for (let i = 0; i < array.length; i++)
        temp[i] = array[i] === 0 ? 1 : 0;
    return temp;
}

function revers(array, s) {
    var temp = [];
    for (let i = 0; i < s; i++)
        temp[i] = array[array.length - i - 1];
    return temp;
}

function sequence(fn, fn1) {
    if (fn1 === 3) return [1];
    fn1 = fn - fn1;
    fn = fn - fn1;
    var array = sequence(fn, fn1);
    var a0 = invers(array);
    var a1 = (fn1 % 2 === 0) ? [1] : [];
    var a2 = revers(array, Math.floor((fn - fn1) / 2));
    return a0.concat(a1, a2);
}

function sequenceFibonacci(iterations) {
    let f0a = [0];
    let f1a = [0];
    for (let i = 0; i < iterations; i++) {
        let f0 = f0a.length;
        let a2 = revers(f1a, f0);
        if (f1a.length % 2 === 0) a2[0] ^= 1;
        f0a = f1a;
        f1a = f1a.concat(a2);
    }
    let array = [];
    for (let i = 0; i < Math.floor(f1a.length / 2); i++)
        array[i] = f1a[i * 2];
    return array;
}

These constructions reproduce the symbolic edges of Fibonacci-based patterns and can be interpreted as recursive encoding schemes derived directly from the observed geometry.

"Toward Generalization"

While the above generators are constructed specifically for Fibonacci-sized rectangles, preliminary experiments suggest that similar structures may emerge for other co-prime pairs. These systems may obey different symbolic transformation rules, but exhibit comparable recursive or compressible traits. A formal generalization of these behaviors remains an open area of exploration.

One of the central challenges that motivated the progression from the original 2013 construction to the deeper analysis in 2019 was the question of irrational proportions: what happens when the rectangle’s side lengths form a truly irrational ratio, such as (1 : φ), rather than an integer-based approximation like 13 : 21?

While recursive generators such as sequence(fn, fn1) accurately reproduced the symbolic boundary sequences for Fibonacci-based rectangles, they were inherently tied to integer dimensions. The challenge was clear: how can one generate the same structures when no exact grid alignment is possible - when the trajectory no longer closes?

This question defines the next stage of the investigation. To answer it, we will analyze the boundary sequences themselves - the so-called fractal sequences - and show how they encode the entire 2D pattern. We will show that these sequences - far from being edge artifacts - contain enough information to deterministically reconstruct the entire 2D pattern. This finding enables a powerful dimensional reduction: the entire billiard system can be expressed as a 1D sequence.

Binary Billiards

We now shift from the dashed-line visualization to a binary representation. Instead of drawing the trajectory, we color the cells the ball passes through, alternating black (0) and white (1) with each step.

Given a rectangle with side lengths
$M$ and $N$ , the ball is launched from a corner and follows diagonal motion, reflecting off the walls. Each step alternates the internal binary state.

The reflection rule causes the trajectory to shift by one cell after each bounce. This alternation creates a consistent visual structure.

When
$M$ and $N$ are coprime, the trajectory visits every cell exactly once:

JavaScript implementation

binarypattern.js

"Gif"

JavaScript implementation

binarypattern_point.js

If the dimensions share a common divisor ( $g c d (M, N) > 1$ ), the trajectory terminates at a corner before filling all cells:

In this case, the system is equivalent to a billiard in a reduced rectangle with dimensions ( $MGCD\frac{M}{GCD}$ , $NGCD\frac{N}{GCD}$ ):

Boundary Behavior and Symmetry

In the coprime case, the ball crosses every row and column. Notably, each pass between the top and left walls consists of an even number of steps.

From this, we can observe a critical symmetry: the left column contains the inverted bits of the top row, excluding the initial bit.

Furthermore, every second bit ( $2_{n-1}$ ) in the top sequence is the inverted version of its neighbor ( $2_{n}$ ). Therefore, we can discard every second bit and retain full pattern information:

For example, with dimensions $M = 21, N = 13$ , the resulting sequence is: 1010010110

This sequence is unique for every coprime pair $(M, N)$ . It encodes all necessary information about the pattern.

Sequence Interpretation

The trajectory between two reflections from the upper wall is always $2 N$ cells long. Each such pass begins with a black cell (bit = 0) and ends with a white cell (bit = 1):

More formally:

A bit of 1 indicates that the ball arrived from a reflection off the right wall
A bit of 0 indicates it came from the left wall

This mapping gives the sequence its meaning. In the diagram below, the trajectory is colored black when moving right and white when moving left:

Encoding Rational Division via Billiards

A curious side effect of the billiard construction is that it naturally encodes binary division of two numbers. Specifically, by tracking the direction of the billiard ball’s movement at each wall collision, and sampling this information at exponentially increasing intervals, one can extract the binary digits of a rational fraction.

Let the billiard table have side lengths $M$ and $N$ , and let the ball bounce between corners. At each collision with the top or bottom wall:

If the ball is moving to the right, record a 0

If moving to the left, record a 1

Then, at every $2^n$ -th collision (i.e., 1st, 2nd, 4th, 8th, …), we record the state.

Example: for a table of size $M = 21, N = 13$ , we obtain:

1st  (bottom, →): 0  
2nd  (top, ←):    1  
4th  (top, →):    0  
8th  (top, →):    0  
16th (top, ←):    1  
32nd (top, ←):    1  
...

This produces the binary expansion:

0.1001111001111001111…
Which is precisely the binary representation of $1321\frac{13}{21}$ .

Reconstruction from the Sequence

The full billiard pattern can be reconstructed from this single boundary sequence. Even extrapolation beyond the grid is possible.

Let us begin by placing the bits along the top edge of a square grid of width $M$ . Bits are spaced every 2 cells - these are the points where the ball would hit the upper wall:

Then:

If the bit is 1, we extend a diagonal to the left
If the bit is 0, we extend it to the right

The first bit (bit 0) is treated specially - it begins the pattern:

The reconstruction produces the exact original pattern:

JavaScript implementation

visualizer.js

This result shows that the 1D sequence contains complete information about the original 2D billiard pattern.

But we're not done.

From the surface of the river, we reduced the system to a rectangular billiard with a dashed diagonal trajectory. Then we reduced it further - to a binary field generated by alternating internal states. Now, we push the reduction one step further: we collapse the entire 2D billiard into a one-dimensional rule. A symbolic system with no geometry left - only structure.

This is where we begin to uncover the origin of these fractals.

One-dimensional billiards

On the $X$ number axis, we take two points: $0$ and $M$ .

Moving from one point to another, we measure the distances $N$ :

We marked a point. We continue to measure the distance from this point, maintaining the direction. If we reach point $0$ or $M$ , we change the direction:

As you can see in the pictures above, the first point shows the place where the ball touches the bottom wall of the billiard table. We are not interested in this point. We will only mark the points $2 k N$ for $k = 0, 1, 2, \dots$ .

How to mark these points? Let's unfold our billiard table on the $X$ axis. Let's mark the points $0, M, 2 M, 3 M, \dots$ . Now, having reached point $M$ , we do not change the direction of movement, but continue moving to point $2 M$ .

Points that are multiples of $M$ divide our axis into segments. We will conditionally mark these segments with ones and zeros (alternating). On the segments marked with zeros, the ball (in rectangular billiards) moves from left to right. On the segments marked with ones, it moves from right to left. Or more simply: the ball moves from left to right if $Q_k=0$ , for

$2);k=0,1,2,…Q_k=\lfloor \frac{2kN}{M} \rfloor \; (\textrm{mod} \; 2); \quad k=0,1,2,…$

It is easy to see that the point at which the ball touched the upper wall of the billiard table is the remainder of dividing $2 k N$ by $M$ . In this case, we don't need to record the movement of the ball in the opposite direction. We take the integer part of dividing $2 k N$ by $M$ , if it is even - we calculate the remainder of dividing $2 k N$ by $M$ . We divide the resulting remainder by 2 (the distance between adjacent touch points is two cells). This gives us the indices of the array elements that correspond to rightward motion (zeros). All other entries - representing leftward trajectories - are filled with ones.

Sequence length = $M2\frac{M}{2}$ .

function sequence(m,n){
    var md=m/2;
    var array=[];
    for(var k=0;k<md;k++) array[k]=1;
    for(var k=0;k<md;k++) if(Math.floor(2*k*n/m)%2==0) array[((2*k*n)%m)/2]=0;
    return array;
}

Now we can build a binary sequence for billiards with any sides $M$ and $N$ (natural numbers). Some examples:
144x89 (Fibonacci numbers):
010100101101001011010110100101101001010010110101001011010100101

169x70 (Pell numbers):
010101101010010101101010010101101010100101011010101001010110101001010101010010101101010010

233x55 (odd Fibonacci numbers $F_n$ and $F_{n-3}$ ):
010010011011011001001101101100100100110110010010011011001001001101101100100110110110010010011011001001001101100100100

This is interesting

Very curious graphs are obtained if you take a billiard table with width M and construct sequences for each N from 0 to M. Then these sequences are stacked.

    var array;
    for(var y=1;y<m;y++){
        array=sequence(m,y);
        for(var x=0;x<array.length;x++){
            if(array[x]==0) context.fillRect (x, y, 1, 1);
        }
    }

Some examples.

M=610	M=611	M=612	M=613	M=614

JavaScript implementation

sequences.js

We have sequences. How else can we visualize binary sequences? With Turtle graphics.

Turtle graphics

The sequence length determines the complexity of the curve. The more irrational the ratio between M and N, the more non-periodic and fractal - like the structure becomes.

Draw a segment. Then take bits from our sequence one by one. If bit = 1 - rotate the segment relative to the previous one by $60∘60^{\circ}$ (clockwise). If bit = 0 - rotate the segment by $−60∘-60^{\circ}$ . The beginning of the next segment is the end of the previous one.

Take two large enough Fibonacci numbers: $F_{29}=514229$ and $F_{28}=317811$ . This ensures that the pattern is long enough for the fractal structure to become apparent, but still bounded enough for visualization.

We built the sequence:
0010110100101101001010010110100101101010010110100101101001010010110100101… (257114 symbols plus a zero bit).

Visualize using turtle graphics. The size of the initial segment is 1 pixels (the initial segment is in the lower right corner):

The next example is Pell numbers.

P_n = \begin{cases} 0, & n=0; \ 1, & n=1; \ 2P_{n-1} + P_{n-2}, & n>1. \end{cases}

We take $P_{16}=470832$ and $P_{15}=195025$ .

The sequence is:
00101001010110101001010101010100101011010100101010110101001010101101 (235415 symbols plus a zero bit).

The size of the initial segment is 1 pixel:

Another example is the odd Fibonacci numbers $F_n$ and $F_{n-3}$ .
Let's take $F_{28}=317811$ and $F_{25}=75025$ .
The sequence is:
00110110010010011011001001101101100100110110110010011011011001001… (158905 plus a zero bit).
Instead of the angles $60∘60^{\circ}$ and $−60∘-60^{\circ}$ , we will use the angles $90∘90^{\circ}$ and $−90∘-90^{\circ}$ .
The size of the initial segment is 1 pixels:

This curve is called "Fibonacci Word Fractal". The Hausdorff dimension for this curve is known:

$D=3log⁡Φlog⁡(1+2)=1.6379;Φ=1+52D=3{\frac {\log \Phi }{\log(1+{\sqrt {2}})}}=1.6379; \quad \Phi =\frac {1+{\sqrt {5}}}{2}$

JavaScript implementation

turtle.js

The problem unsolved yet.

Is it possible to draw a pattern for billiards, the sides of which are incommensurable (one of the sides is an irrational number)? At first glance, the task seems impossible. Trying to solve this problem, we will face a number of questions:

If the sides are incommensurable, we cannot tile the billiards with cells of the same size.
If the sides are incommensurable, the ball will infinitely reflect and will never get into the corner.
Sequences in billiards are not filled in order, but chaotically.

The first two questions obviously have no solution. But if there were a way to fill the sequence in order, then we could, moving along the sequence from left to right, restore the pattern in the way we used above. And thus see what the pattern looks like in the upper left corner of the billiard table whose sides are incommensurable.

o_O

Let's take a billiard table, the sides of which are equal to the Fibonacci numbers (this trick may not work with other numbers). Let's throw a ball into it and record the number of times the ball touches the upper wall. We'll paint the numbers white if the ball moved from right to left and black if the ball moved from left to right:

White corresponds to the number one in the sequence, black to zero. Now let's arrange the numbers in order:

We've got exactly the same sequence of ones and zeros.

21(1), 13(0), 8(1), 26(0), 5(0), 16(1), 18(0), 3(1), 24(1), 10(0), 11(1), 23(0), 2(0), 19(1), 15(0), 6(1), 27(1), 7(0), 14(1), 20(0), 1(1), 22(1), 12(0), 9(1), 25(0), 4(0), 17(1)

1(1), 2(0), 3(1), 4(0), 5(0), 6(1), 7(0), 8(1), 9(1), 10(0), 11(1), 12(0), 13(0), 14(1), 15(0), 16(1), 17(1), 18(0), 19(1), 20(0), 21(1), 22(1), 23(0), 24(1), 25(0), 26(0), 27(1)

"For other numbers"

The origin is the upper left corner. The X axis is the width of the billiard table M. The Y axis is the height of the billiard table N. The numbers for which the sequences match are marked with white dots:

Numbers for which the sequence is inverted:

JavaScript implementation:

The first line is the mouse coordinates, which are used as the width and height of the billiard table.
The second line is the first 100 bits of the sequence obtained through the remainders of the division.
The third line is the first 100 bits of the sequence obtained through the parity of the integer part.

Black color - Visualization of the first sequence using Turtle graphics.
Purple - visualization of the second sequence.

JavaScript implementation

turtle_dynamic.js

In fact, in some cases, we do not need to take the remainder from the division. For Fibonacci numbers, it is enough to check the parity of the integer part of the division of $2 k N$ by $M$ :

$2);k=0,1,2,…Q_k=\lfloor \frac{2kN}{M} \rfloor \; (\textrm{mod} \; 2); \quad k=0,1,2,…$

In the numerator we have $F_{n}$ . In the denominator - $F_{n+1}$ .

As is known:

$lim⁡n→∞FnFn+1=1Φ\lim_{n\to\infty} \frac{F_n}{F_{n+1}}= \frac{1}{\Phi}$

This gives us a bridge: rational billiards ( $F_{n}$ , $F_{n+1}$ ) converge toward an irrational limit - $Φ\Phi$ - allowing us to define an infinite symbolic sequence.

$Φ\Phi$ is the Golden Ratio. An irrational number. Now we can write our formula as:

$2);k=0,1,2,…Q_k=\lfloor \frac{2k}{\Phi} \rfloor \; (\textrm{mod} \; 2); \quad k=0,1,2,…$

We have obtained a formula with which we can fill in the sequence for a billiard table, the width of which is $Φ\Phi$ and the height is $1$ . The length of the sequence = $∞\infty$ , but we can restore part of the pattern by moving from left to right along the sequence and looking into the upper left corner of the billiard table. It remains to figure out how to calculate $Φ\Phi$

One divided by the golden ratio can be rewritten as:

$1Φ=−1+52\frac{1}{\Phi}=\frac {-1+{\sqrt {5}}}{2}$

We can get rid of the two:

$2kΦ=2k(−1+5)2=k5−k\frac{2k}{\Phi}=\frac {2k(-1+{\sqrt {5}})}{2}=k\sqrt{5}-k$

Our formula takes the form:

$2);k=0,1,2,…Q_k=\lfloor k\sqrt{5}-k \rfloor \; (\textrm{mod} \; 2); \quad k=0,1,2,…$

Now we can draw part of the billiard pattern with sides $1$ and $Φ\Phi$ :

If we do not subtract $k$ each time, then every second bit in the sequence is inverted. We get the general formula:

$2);k=0,1,2,…Q_k=\lfloor k\sqrt{x} \rfloor \; (\textrm{mod} \; 2); \quad k=0,1,2,…$

Let's build a sequence for $k2k\sqrt{2}$

var x=2;
var q=[];
for(var k=0;k<256000;k++) q[k]=Math.floor(k*Math.sqrt(x))%2;

The first few bits of the sequence (A083035):
01001101100100110011001101100100110110011001001101100100110110…

Angles are $90∘90^{\circ}$ and $−90∘-90^{\circ}$ . The size of the initial segment is 5 pixels:

"This is interesting"

From this curve, we can reconstruct the "billiard pattern" and see what is around the curve:

It would be interesting to find M and N for this pattern.

Number of segments in the repeating part of the curve = $P_n$ (Pell numbers: 0, 1, 2, 5, 12, 29, 70, 169, 408, 985, 2378, … ).

$2=lim⁡n→∞Pn−1+PnPn\sqrt{2} = \lim_{n\to\infty} \tfrac{P_{n-1}+P_n}{P_n}$

Someone may doubt that the parity of the integer part of $k2k\sqrt{2}$ gives a fractal sequence. Let's visualize part of this sequence with the visualizer described earlier:

For clarity, we colored the longest curve in the resulting pattern:

This curve has a name - "Fibonacci word fractal".

Thus, by gradually reducing billiard geometry through symbolic encodings, we arrive at a powerful realization: even irrational systems, which defy spatial tiling and corner reflection, can still produce deterministic, fractal structure - using only integer math.

Filling the Fractal: From Symbolic Structure to Visual Density

In the previous sections, we showed how symbolic sequences can generate complex boundary structures in discrete 2D space. All of these patterns - whether generated by rational billiards or floor-based symbolic systems - form enclosed regions.
Some of these regions close against the boundaries of the rectangle, while others are fully self-contained, looping within the grid. In either case, the resulting trajectories always define fully enclosed cells.

Consider the following example:

A binary sequence generated by the floor function:

$2);n=0,1,2,…Q_n=\lfloor n\sqrt{2} \rfloor \; (\textrm{mod} \; 2); \quad n=0,1,2,…$

0100110110010011001001101100

Using the visualizer in diagonal mode, we generate a familiar fractal boundary.

However, switching to horizontal-vertical visualization, we invert every even-indexed bit and plot dashed lines accordingly. For bits with value 0, the lines are offset by one unit.

Vertical lines:

Horizontal lines:

Merged:

This technique resembles Hitomezashi stitching, a traditional Japanese method of generating textile patterns. While historically used to create decorative designs, this form of boundary generation bears striking similarity to the symbolic outlines produced in our system.

However, Hitomezashi does not solve the filling problem. It provides only the skeleton - the structure of edges - but no mechanism for interior construction.

A Symbolic Filling Algorithm

To solve this, we developed a symbolic method to automatically fill the interior regions - using only the original sequence. No region detection, no geometry, no search algorithms.

The method works as follows:

Construct a cumulative array a[n], based on the bitwise value of the fractal sequence:

We define:

Q_n = \lfloor n\sqrt{x} \rfloor \bmod 2, \quad n=0,1,2,\dots \

and

a_n = \begin{cases} a_{n-1} + 1, & Q_n = 1;\ a_{n-1} - 1, & Q_n = 0. \end{cases}

In code:

var a = [0];
for (var i = 1; i < size; i++) {
    if (Math.floor(i * Math.sqrt(2)) % 2 == 1)
        a[i] = a[i - 1] + 1;
    else
        a[i] = a[i - 1] - 1;
}

Then, for each cell (x, y), compute:

q = (a[x] + a[y] + 512) % 4;
if (q === 0 || q === 1)
    context.fillRect(x, y, 1, 1);

The result is a filled pattern - not derived from pixel analysis or region marking, but emerging directly from the same symbolic system that generated the boundaries.

Full Algorithm

var a=[0];
for(var i=1;i<size;i++){
    if(Math.floor(i*Math.sqrt(2))%2==1)
        a[i]=a[i-1]+1;
    else
        a[i]=a[i-1]-1;
}
for(var x=0;x<size;x++){
    for(var y=0;y<size;y++){
        q=(a[x]+a[y]+512)%4;
        if(q==0 || q==1) context.fillRect(x, y, 1, 1);
    }
}

Hitomezashi and Symbolic Filling Algorithm

Examples

Fractal fill based on:

$2);n=0,1,2,…Q_n=\lfloor n\sqrt{2} \rfloor \; (\textrm{mod} \; 2); \quad n=0,1,2,…$

"Gif"

$2);n=0,1,2,…Q_n=\lfloor n(\sqrt{5}+1) \rfloor \; (\textrm{mod} \; 2); \quad n=0,1,2,…$

(Fibonacci-based)

"Gif"

Interactive

Visualizer - Fractal Fill from Irrationals

fractal.js

Dynamic Visualization - Using Rational Approximations

fractal_dynamic.js

Perfect Shuffle Comparison

Pattern Matching via Perfect Shuffles

We can compare fractal sequences from floor-based systems with sequences produced by a perfect shuffle function.

Here's our shuffle logic:

function shuffle(array, shiftAmount) {
    let len = array.length;
    let shuffled = new Array(len * 2);
    for (let i = 0; i < len; i++) {
        shuffled[2 * i] = array[(i + shiftAmount) % len];
        shuffled[2 * i + 1] = array[i];
    }
    return shuffled;
}

We start with a base array:

let shuffleIterations = 2*7+1;
let powerOfTwo = 2**7;
let shiftAmount = y * powerOfTwo;
let array1 = [1, 0];
for (let i = 0; i < shuffleIterations; i++) {
    array1 = shuffle(array1, shiftAmount);
}

And compare it with this floor-based sequence:

let powerOfTwo = 2**7;
let irrationalApproximation = y / powerOfTwo;
let array2 = [];
for (let i = 0; i < sizexy; i++) {
    array2[i] = Math.floor(i * irrationalApproximation) % 2;
}

Both arrays are then visualized using our bit-pattern analysis:

for (let y = 0; y < sizexy; y++) {
    // generate array1 or array2

    let digit;
    let bits = [];
    for (let i = 0; i < len; i++) {
        digit = 0;
        for (let j = 0; j < bitLength; j++) {
            digit |= array[i + j] << (bitLength - 1 - j);
        }
        if (!bits.includes(digit)) {
            bits.push(digit);
        }
    }
    for (let i = 0; i < bits.length; i++) {
        context.fillRect(bits[i] * size, y * size, size, size);
    }
}

Here are the results:

Perfect shuffle:

Floor-based sequence:

let pow=7;
let shuffleIterations = 2*pow+1;
let powerOfTwo = 2**pow;

let map1=[];
for(let y=0;y<sizexy;y++){
    map1[y]=[];
    let shiftAmount = y * powerOfTwo;
    let array = [0, 1];
    for (let i = 0; i < shuffleIterations; i++) {
        array = shuffle(array, shiftAmount);
    }
    map1[y]=array;
}
drawMap(document.getElementById('myCanvas'), map1);

map2=[];
for(let y=0;y<sizexy;y++){
    let irrational = y / powerOfTwo;
    let array2=[];
    for (let i=0;i<sizexy;i++){
        array2[i]=Math.floor(i * irrational)%2;
    }
    map2[y]=array2;
}
drawMap(document.getElementById('myCanvas1'), map2);

We generate two binary images using two completely different methods.

In the first image, each row is created by starting with [0, 1] and applying a recursive "perfect shuffle" operation. This operation doubles the array size each time and mixes a shifted copy of the previous state with itself. The shift amount depends on the current row number. We repeat this shuffle a specific number of times for each row to reach the final size.

In the second image, each row is generated by multiplying each column index by a scaled version of the row number, then taking the floor and reducing it modulo 2.

The two outputs are not just visually similar—they are bitwise identical, pixel by pixel, after horizontal flipping.

We verified this not just through floating-point floor operations, but also with a pure integer formulation using bit shifts:

array2[i] = ((i * y) >> pow) % 2;

This version avoids any irrational approximations entirely. It shows that the mapping:

$2B_{y}(i) = \left\lfloor \dfrac{i \cdot y}{2^{\text{pow}}} \right\rfloor \bmod 2$

is equivalent to a symbolic shift-and-fold recursion using perfect shuffles:

[0, 1] is shuffled 2·pow + 1 times
Each shuffle step doubles the array length
Shift is y × 2^pow each time

We also observe that the symbolic structure of the shuffle is time-reversed compared to the floor-based sequence, which explains the horizontal flipping.

This means the perfect shuffle method is not just a metaphor or numerical coincidence — it's an algebraic encoder for binary floor-based mappings.

Theorem (Symbolic Equivalence)

Let:

$\in \mathbb{N}$
$N = 2^{2p+1}$
$Shuffley[2p+1]\text{Shuffle}_{y}^{[2p+1]}$ denote the result of applying a recursive perfect shuffle with shift $\cdot 2^p$ , starting from $[0, 1]$
$2)\text{BinaryFloor}_{y}(x) = \left( \left\lfloor \dfrac{x \cdot y}{2^p} \right\rfloor \bmod 2 \right)$

Then:

$Shuffley[2p+1](x)=1−BinaryFloory(N−x−1)\boxed{\text{Shuffle}{y}^{[2p+1]}(x) = 1 - \text{BinaryFloor}{y}(N - x - 1)}$

This relation holds exactly, using pure bitwise arithmetic.

See:
shuffle_and_billiard_isomorphism.js

Nonlinear Functions and Symbolic Surface Slicing

Symbolic Rotation and the Circle Representation

The binary sequence generated by

$2\left\lfloor k \sqrt{2} \right\rfloor \bmod 2$

can be interpreted geometrically as a rotation on the unit circle. At each step $k$ , we rotate by an irrational angle (e.g., $2\sqrt{2}$ ) and assign a 0 or 1 depending on which half of the circle the point lands in.

An alternative but mathematically equivalent formulation is:

$sin⁡(πk2)>0\sin\left( \pi k \sqrt{2} \right) > 0$

The sine function is positive in one half of the circle and negative in the other - effectively producing the same binary threshold.

This circular diagram shows the orbit of

$\sqrt{2} \bmod 1$

projected onto the unit circle. The points rotate continuously and never repeat. Each crossing of the horizontal axis (i.e., sign change of $sin⁡(πk2)\sin\left( \pi k \sqrt{2} \right)$ ) causes a binary switch. The result is a Sturmian sequence - a symbolic encoding of irrational rotation - visualized as a point traveling clockwise and recording its hemisphere.

Interestingly, even a sequence like

$sin⁡(k)>0\sin\left( k \right) > 0$

generates a similarly complex structure. Though the function seems simple, stepping by exactly 1 radian per iteration is irrational with respect to the sine wave’s natural period of $\pi$ . Thus:

$2\sin(k) > 0 \quad \Longleftrightarrow \quad \left\lfloor \frac{k}{\pi} \right\rfloor \bmod 2$

Which means: irrational step size alone is enough to create symbolic fractals. It’s not the shape of the function - it’s the incommensurability that matters.

Geometric Interpretation of Linear Discretization

We begin with a geometric interpretation of the function

$2\left\lfloor k \sqrt{2} \right\rfloor \bmod 2$

This can be visualized as a straight line $\sqrt{2}$ traversing a periodic 2D space where the $y$ -axis alternates between bands labeled 0 and 1. The value of the function at each step depends on the integer part of $k2k\sqrt{2}$ , determining whether the line intersects an even or odd band.

Alternatively, we may interpret the system as a function $y = x$ with $2\sqrt{2}$ acting as a discretization step along the vertical axis. This perspective reduces the behavior of the original billiard system to a symbolic sampling of a continuous linear function.

Both

$2andsin⁡(πk2)>0\left\lfloor k \sqrt{2} \right\rfloor \bmod 2 \quad \text{and} \quad \sin\left( \pi k \sqrt{2} \right) > 0$

serve as symbolic discretization methods - one yielding binary values via integer floor division, the other through sign thresholding on a continuous sinusoid.

Transition to Quadratic Functions

Having reduced the system to a linear function and its discretization, we now explore the effects of replacing the base function $y = x$ with a nonlinear alternative.

Linear rotation yields symbolic Sturmian sequences. What happens if we increase the curvature? Quadratic growth is the simplest next step.

We now substitute the linear function with:

$y = x^2$

We construct the new sequence using:

$2);x=0,1,2,…Q_x=\lfloor x^2\sqrt{2} \rfloor \; (\textrm{mod} \; 2); \quad x=0,1,2,…$

q[x] = Math.floor(x * x * Math.sqrt(2)) % 2;

This substitution shifts us from uniform linear growth to a system governed by accelerating curvature, producing richer and more complex symbolic behavior.

Visualizing this sequence with the Symbolic Filling Algorithm yields a disordered pattern.

"Gif"

Similarly, applying Turtle graphics results in "chaotic" outputs

suggesting that the underlying symbolic dynamics differ significantly from the linear case. This observation motivated us to try a simpler visual encoding.

This prompts a reconsideration of visualization methods. Rather than relying on complex binary sequence encodings (e.g., symbolic filling or turtle graphics), we adopt a simpler approach:

Place the first k elements of the sequence q[x] in the initial row.
For each subsequent row along the y-axis, use q[x+ky].
The result is a striking array of structured, interference-like patterns.

k=35:

k=661:

This reveals that a seemingly chaotic 1D sequence can exhibit coherent spatial behavior when unfolded along two dimensions. Generalizing this, we rewrite the indexing expression as:

q[x] = Math.floor((x + k * y) ** 2 * Math.sqrt(2)) % 2;

Expanding the square:

$x + k y)^2 = x^2 + 2kxy + k^2 y^2$

This quadratic form suggests a symbolic sampling of a 3D curved surface. We abstract this into a general expression:

$\left( x^2 + bxy + c y^2 \right)^d$

Here:

$a$ controls vertical scaling (the discretization frequency),
$b$ introduces diagonal shear,
$c$ modulates stretching along the $y$ -axis,
and $d$ introduces curvature nonlinearity.

For example, setting $b = 0, c = 1, d = 1$ yields:

$\left( x^2 + y^2 \right)$

an elliptical paraboloid, a classic bowl-shaped surface.

We then visualize both:

$2\left\lfloor z \sqrt{2} \right\rfloor \bmod 2$

$sin⁡(πz2)\sin\left( \pi z \sqrt{2} \right)$

...using this surface. Despite differing in output (binary vs. continuous), both render structurally equivalent patterns: the sine version produces smooth grayscale textures, while the floor function yields crisp binary segmentation. In either case, the resulting 2D images exhibit radial, interference-like motifs - strongly reminiscent of diffraction or holographic patterns.

binary	continuous

The sine-based rendering reveals finer gradients and smoother interference zones, but the underlying symbolic structure is identical to the binary version.

Holographic Analogy

In classical holography, a point source emits nested spherical wavefronts. When these are sliced by a flat recording plane, the amplitude at each (x, y) position encodes a continuous interference pattern. The spacing between wavefronts defines the wavelength - effectively a discretization step in the z-direction.

Our system inverts this paradigm: instead of slicing nested curved shells with a flat surface, we intersect a single curved surface with stacked, evenly spaced binary planes - symbolic layers representing a plane wavefront. Each (x, y) coordinate samples which symbolic layer the surface intersects, producing a binary (or thresholded) value.

While the components differ - continuous vs symbolic, curved emitter vs curved surface - the underlying mechanism is the same: a curved geometry intersected by layered structure, producing interference-like textures from simple spatial rules.

Despite differences in physical interpretation, both systems share a core structure: curved geometry meets layered slicing - and complexity emerges.

Our picture size is 400 pixels. We used $a=1400a=\frac{1}{400}$ for previous picture to fit the picture.

Here another pattern with $a=1200a=\frac{1}{200}$ :

binary	continuous

$c$ modulates stretching along the $y$ -axis

(we will use continuous patterns from now).

$c = 2$	$c = 0.2$

$b$ introduces diagonal shear

$b = - 1.5$	$b = 1.5$

For $b = 0, c = - 1, d = 1$ our equation becomes:

$z=a(x^2-y^2)$

This is a hyperbolic paraboloid - a surface with negative Gaussian curvature:

$a=1400a=\frac{1}{400}$	$a=1200a=\frac{1}{200}$

For $b = 0, c = 1, d = 1/5$ our equation becomes:

$z=a(x^2+y^2)^{1/5}$


$a = 1$	$a = 10$

$a = 100$	$a = 1000$

The most interesting patterns are obtained if we take a $d$ that differs slightly from 1. For example, for

$d=\frac{1}{1.01}$


$a=18258a=\frac{182}{58}$	$a=173165a=\frac{173}{165}$

$a=170109a=\frac{170}{109}$	$a=11170a=\frac{111}{70}$

$a=57178a=\frac{57}{178}$	$a=186119a=\frac{186}{119}$

$d=\frac{1}{1.01}$


$a=174111a=\frac{174}{111}$	$a=50152a=\frac{50}{152}$

Interactive

Static Visualization

hologram_s.js

Dynamic Visualization

hologram_dynamic.js

Hologram Experiment

hologram_reconstruction.html

Conclusion

This study began with a simple system: a billiard ball reflecting within a rectangular grid. By examining its trajectory and translating spatial behavior into symbolic sequences, we uncovered patterns exhibiting recursive structure, boundary self-similarity, and interference-like properties.

Through progressive abstraction - from 2D reflections to 1D symbolic sequences, and eventually to nonlinear surfaces - we showed that binary patterns can emerge from discretized irrational steps alone. Formulas such as:

$2andQk=sin⁡(πkx)>0Q_k=\left\lfloor k \sqrt{x} \right\rfloor \bmod 2 \quad \text{and} \quad Q_k=\sin\left( \pi k \sqrt{x} \right) > 0$

generate deterministic sequences that yield quasi-fractal geometries when visualized spatially. Though entirely discrete, these systems exhibit behaviors commonly associated with continuous wave phenomena.

Further generalization using surface equations of the form:

$\left( x^2 + bxy + c y^2 \right)^d$

reveals that curvature, shear, and nonlinearity shape the resulting symbolic slices in ways reminiscent of optical interference. The resemblance between these binary patterns and holographic textures suggests a structural parallel, despite their different origins.

No physical claims are made. However, the symbolic systems presented here demonstrate that simple integer-based operations, applied to irrational quantities or curved domains, can generate high-order structure with properties of both visual and mathematical interest. This offers a conceptual bridge between number theory, symbolic dynamics, and emergent spatial organization.

While this study has uncovered significant insights, several open questions remain to inspire further research:

Generalization to Non-Fibonacci Ratios. The recursive generators presented here excel for Fibonacci-based rectangles, but can similar symbolic encoding schemes be developed for arbitrary co-prime pairs or other number sequences (e.g., Lucas numbers, Pell numbers)? What universal properties govern the resulting fractal structures?
Irrational Billiards and Sequence Ordering. For billiard tables with incommensurable sides (e.g., 1 : φ), we partially addressed pattern generation using floor-based sequences. However, can a robust method be devised to order these sequences consistently, enabling full reconstruction of the pattern in the absence of grid alignment?
Perfect Shuffle and Floor-Based Equivalence. The striking similarity between perfect shuffle sequences and irrational floor-based sequences suggests a potential mathematical identity. Is there a formal equivalence between these systems, or do they converge only under specific conditions? What deeper structures might unify these approaches?
Nonlinear Surface Generalizations. The transition to quadratic and higher-order surfaces produced rich, holographic-like patterns. How do variations in the exponent (d) or coefficients (a), (b), and (c) affect the fractal dimension or complexity of the resulting patterns? Can these be analytically quantified?
Applications to Physical Systems. The visual resemblance to holographic interference patterns raises the question: can these symbolic systems model real-world wave phenomena, such as optical diffraction or quantum interference, in a computationally efficient way? These questions highlight the potential for further exploration at the intersection of symbolic dynamics, fractal geometry, and computational visualization. By bridging discrete mathematics with continuous phenomena, this work offers a foundation for investigating how simple rules can encode complex, emergent behaviors - potentially extending to fields like number theory, physics, and computer graphics.

GitHub repository: https://github.com/xcontcom/billiard-fractals

// Detect dark theme var iframe = document.getElementById('tweet-1941606569273155829-837'); if (document.body.className.includes('dark-theme')) { iframe.src = "https://platform.twitter.com/embed/Tweet.html?id=1941606569273155829&theme=dark" }