Recamán's sequence

The Recamán sequence is a famous sequence invented by the Colombian mathematician, Bernardo Recamán Santos. It is defined by the following algorithm, starting at $a_0 = 0$:

$a_n = a_{n-1} - n$ if this number is positive and has not already appeared in the sequence; otherwise $a_n = a_{n-1} + n$.

It was covered in a Numberphile YouTube video that also featured the visualization by Edmund Harriss in which sequential numbers are linked by semicircular arcs. The following Python code reproduces the resulting image.

The test for deciding the next number in the sequence uses Python's new (since 3.8) "Walrus operator". This allows the use of a named expression: (b := a - n) in the conditional statement:

if (b := a - n) > 0 and b not in seen:
    a = b

which assigns the result of a - n to the variable name b so that its value can be reused. In Python 3.7 and earlier, the assignment has to be made on a separate line:

b = a - n
if b > 0 and b not in seen:
    a = b

I'm not wild about named expressions, and wouldn't usually use one in my own code here; I find the second approach clearer and less error prone: if you forget the parentheses, the line if b := a - n > 0 and b not in seen: is syntactically correct but assigns the boolean result of a - n > 0 to b because of the higher precedence of the > operator...

import numpy as np
import matplotlib.pyplot as plt

# Equal aspect ratio Figure with black background and no axes.
fig, ax = plt.subplots(facecolor='k')
ax.axis('equal')
ax.axis('off')

# Colour the lines sequentially from a Matplotlib colormap.
cm = plt.cm.get_cmap('Spectral')

def add_to_plot(lasta, a, n):
    """Add a semi-circular arc from a to lasta.

    Arcs alternate to be above and below the x-axis according to whether
    n is even or odd.

    """

    # Arc centre and radius.
    c, r = (a + lasta) / 2, abs(a - lasta) / 2
    x = np.linspace(-r, r, 1000)
    y = np.sqrt(r**2 - x**2) * (-1)**n
    color = cm(n/max_terms)
    ax.plot(x+c, y, c=color, lw=1)

# Keep track of which numbers have been "visited" in this set.
seen = set()
n = a = 0
lasta = None
max_terms = 60
while n < max_terms:
    lasta = a
    if (b := a - n) > 0 and b not in seen:
        a = b
    else:
        a = a + n
    seen.add(a)
    n += 1
    print(a)
    add_to_plot(lasta, a, n)

plt.show()

If you like your Recamán sequence visualized diagonally, modify the add_to_plot function to rotate by 45°:

def add_to_plot(lasta, a, n):
    c = (a + lasta) / 2
    r = abs(a - lasta) / 2
    x = np.linspace(-r, r, 1000)
    y = np.sqrt(r**2 - x**2) * (-1)**n
    color = cm(n/max_terms)
    X = x+c
    f = np.sqrt(2)
    ax.plot((X-y)*f, (X+y)*f, c=color, lw=1)

Recamán's sequence

Posted by: christian on 29 Feb 2020

Comments

Nathan 4 years, 9 months ago

christian 4 years, 9 months ago

New Comment