In [266]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [267]:
%reload_ext watermark
%watermark -v -m --iversions
pandas     0.25.3
seaborn    0.10.0
matplotlib 3.1.1
numpy      1.17.2
CPython 3.7.3
IPython 7.8.0

compiler   : GCC 8.0.1 20180414 (experimental) [trunk revision 259383
system     : Linux
release    : 4.15.0-74-generic
machine    : x86_64
processor  : x86_64
CPU cores  : 4
interpreter: 64bit

Helpers

Define functions to draw the plots

In [268]:
def show_scatter(a):
    plt.figure(figsize=(8,8))
    plt.scatter(range(len(a)), a, s=10, c='black', marker=".")
    plt.show()
In [269]:
def show_bitmatrix(a):
    plt.figure(figsize=(8,8))
    plt.imshow(a, cmap='Greys',  interpolation='nearest')
    plt.xticks([])
    plt.yticks([])
    plt.show()

Binary Representation

Binary plot on mathworld.wolfram.com

A binary plot of an integer sequence is a plot of the binary representations of successive terms where each term is represented as a sequence of bits with 1s colored black and 0s colored white. Then each representation is stacked to form a table where each entry represents a bit in the sequence. To make a binary plot of a sequence we need to convert each term of the sequence into its binary representation. Then we have to put this representation in a form which make us able to build the plot easily. The following function converts the sequence in a matrix where the element i-j represents is the bit j of the element i of the sequence:

In [270]:
def seq_to_bitmatrix(s):
    # maximum number of bits used in this sequence
    maxbit = len(bin(max(s)))
    M = np.zeros((len(s),maxbit),dtype='uint8')
    for i,e in enumerate(s):
        bits = bin(e)[2:] # bin() -> 0b100, we drop 2 chars
        for j,b in enumerate(bits):
            M[i,maxbit-len(bits)+j] = int(b) # fill the matrix
    return M   

Positive integers

0, 1, 2, 3...

In [271]:
posint_a = list(range(50))
In [272]:
show_scatter(posint_a)
show_bitmatrix(seq_to_bitmatrix(posint_a).T)

Integer squares

0, 1, 4, 9, 16, 25...

In [273]:
sqint_a = [pow(n, 2) for n in range(50)]
In [274]:
show_scatter(sqint_a)
show_bitmatrix(seq_to_bitmatrix(sqint_a).T)

Integer square root

0, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3,

In [275]:
sqrtint_a = [int(np.sqrt(n)) for n in range(500)]
In [276]:
show_scatter(sqrtint_a[:50])
show_scatter(sqrtint_a)
show_bitmatrix(seq_to_bitmatrix(sqrtint_a[:50]).T)

Powers of 2

In [277]:
_ = [int(2**n) for n in range(1000)]
In [278]:
show_scatter(_)

n*(n+1)/2 sequence

n = n*(n+1)/2

0, 1, 3, 6, 10, 15, 21...

In [279]:
s = [int(n*(n+1)/2) for n in range(100)]
In [280]:
show_scatter(s)
show_bitmatrix(seq_to_bitmatrix(s).T)

Prime numbers

A000040

The prime numbers.

In [281]:
def gen_primes():
    D = {}
    q = 2
    while True:
        if q not in D:
            yield q
            D[q * q] = [q]
        else:
            for p in D[q]:
                D.setdefault(p + q, []).append(p)
            del D[q]
        q += 1
In [282]:
f = gen_primes()
pa = [next(f) for i in range(100)]
In [283]:
show_scatter(pa)
show_bitmatrix(seq_to_bitmatrix(pa).T)

Fibonacci numbers

A000045

Generates the Fibonacci numbers, starting with 0

In [284]:
def fib(): 
    x, y = 0, 1 
    while 1: 
        yield x 
        x, y = y, x+y
In [285]:
f = fib()
fib_a = [next(f) for i in range(200)]
In [286]:
show_scatter(fib_a)
show_bitmatrix(seq_to_bitmatrix(fib_a))

Pascal's triangle

A007318

read by rows: C(n,k) = binomial(n,k) = n!/(k!*(n-k)!), 0 <= k <= n.

In [287]:
from numpy import prod;
def C(n, k): 
    return prod([(n-j)/(j+1) for j in range(k)])
In [288]:
res=1
for r in range(15):
    print(*["{}".format('{:>5}'.format(int(C(r, i)))) for i in range(res)])
    res+=1
    1
    1     1
    1     2     1
    1     3     3     1
    1     4     6     4     1
    1     5    10    10     5     1
    1     6    15    20    15     6     1
    1     7    21    35    35    21     7     1
    1     8    28    56    70    56    28     8     1
    1     9    36    84   126   126    84    36     9     1
    1    10    45   120   210   252   210   120    45    10     1
    1    11    55   165   330   461   461   329   164    54    10     0
    1    12    66   220   495   792   924   792   495   220    66    12     1
    1    13    78   286   715  1287  1716  1716  1287   715   286    78    13     1
    1    14    91   364  1001  2002  3003  3432  3003  2002  1001   364    91    14     1
In [289]:
from itertools import islice

def pascal():
    row = [1]
    while True:
        yield row
        row = [i+j for i,j in zip([0]+row, row+[0])]

for r in islice(pascal(),0,11):
    print("{:^60}".format("".join(map(lambda n: "{:>5}".format(n), r))))
                               1                            
                             1    1                         
                          1    2    1                       
                        1    3    3    1                    
                     1    4    6    4    1                  
                   1    5   10   10    5    1               
                1    6   15   20   15    6    1             
              1    7   21   35   35   21    7    1          
           1    8   28   56   70   56   28    8    1        
         1    9   36   84  126  126   84   36    9    1     
      1   10   45  120  210  252  210  120   45   10    1   
In [290]:
def C(n, k):
    if k<0 or k>n:
        return 0
    res=1
    for i in range(k):
        res=res*(n-i)//(i+1)
    return res
In [291]:
res=1
for r in range(11):
    print(*["{}".format('{:>5}'.format(int(C(r, i)))) for i in range(res)])
    res+=1
    1
    1     1
    1     2     1
    1     3     3     1
    1     4     6     4     1
    1     5    10    10     5     1
    1     6    15    20    15     6     1
    1     7    21    35    35    21     7     1
    1     8    28    56    70    56    28     8     1
    1     9    36    84   126   126    84    36     9     1
    1    10    45   120   210   252   210   120    45    10     1

Concatenate lists to look at it as number sequence

In [292]:
_ = [int(C(i, j)) for i in range(15) for j in range(15) if int(C(i, j)) != 0]
show_scatter(_)

Stern's diatomic series

A002487

Stern's diatomic series (or Stern-Brocot sequence): a(0) = 0, a(1) = 1; for n > 0: a(2n) = a(n), a(2n+1) = a(n) + a(n+1).

In [293]:
def stern_n(n): 
    return n if n<2 else stern_n(n/2) if n%2==0 else stern_n((n - 1)/2) + stern_n((n + 1)/2)
In [294]:
stern_a = [int(stern_n(i)) for i in range(5000)]
In [295]:
show_scatter(stern_a)

Fly straight, dammit series

A133058

1, 1, 4, 8, 2, 8, 4, 12, 3, 1, 12, 24, 2, 16, 8, 24, 3, 21, 7, 27, 48, 16...

a(0)=a(1)=1; for n>1, a(n) = a(n-1) + n + 1 if a(n-1) and n are coprime, otherwise a(n) = a(n-1)/gcd(a(n-1),n).

In [296]:
from math import gcd
ra = {}
for n in range(1000):
    if n<=1:
        ra[n]=1
    elif gcd(int(ra[n-1]), n)==1:
        ra[n]= int(ra[n-1]+n+1)
    else:
        ra[n] = int(ra[n-1]/gcd(int(ra[n-1]), n))
In [297]:
show_scatter(list(ra.values()))

Hofstadter Q-sequence

A005185

1, 1, 2, 3, 3, 4, 5, 5, 6, 6, 6, 8, 8, 8, 10, 9, 10, 11...

Hofstadter Q-sequence: a(1) = a(2) = 1; a(n) = a(n-a(n-1)) + a(n-a(n-2)) for n > 2.

In [298]:
hofs = {}
for n in range(10000):
    if n<=2:
        hofs[n] = 1
    else:
        hofs[n] = hofs[n-hofs[n-1]]+hofs[n-hofs[n-2]]
In [299]:
show_scatter(list(hofs.values()))

Hofstadter-Conway sequence

A004001

1, 1, 2, 2, 3, 4, 4, 4, 5, 6, 7, 7, 8, 8, 8, 8, 9, 10, 11, 12, 12, 13, 14, 14...

Hofstadter-Conway $10000 sequence: a(n) = a(a(n-1)) + a(n-a(n-1)) with a(1) = a(2) = 1.

In [300]:
hofsconw = {}
for n in range(5000):
    if n<=2:
        hofsconw[n] = 1
    else:
        hofsconw[n] = hofsconw[hofsconw[n-1]]+hofsconw[n-hofsconw[n-1]]
In [301]:
show_scatter(list(hofsconw.values()))

A chaotic cousin of the Hofstadter-Conway sequence

A055748

1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6, 7, 8, 8, 8, 8, 8, 8, 9, 10...

a(1) = 1, a(2) = 1, a(n) = a(a(n-1)) + a(n - a(n-2) - 1) for n >= 3

In [302]:
hofsconwc = {}
for n in range(1000):
    if n<=2:
        hofsconwc[n] = 1
    else:
        hofsconwc[n] = hofsconwc[hofsconwc[n-1]]+hofsconwc[n-hofsconwc[n-2]-1]
In [303]:
show_scatter(list(hofsconwc.values()))

Langton's ant

Wikipedia, Langton's ant.

A274369

Let the starting square of Langton's ant have coordinates (0, 0), with ant looking in negative x-direction. a(n) is the x-coordinate of the ant after n moves.

0, 0, 1, 1, 0, 0, -1, -1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1...

A274370

Let the starting square of Langton's ant have coordinates (0, 0), with the ant looking in negative x-direction. a(n) is the y-coordinate of the ant after n moves.

0, 1, 1, 0, 0, -1, -1, 0, 0, -1, -1, -2, -2, -1, -1, 0, 0, -1, -1, -2...

In [304]:
def ant(n):
    steps = [(1, 0), (0, 1), (-1, 0), (0, -1)]
    black = set()
    x = y = 0
    position = [(x, y)]
    direction = 2
    for _ in range(n):
        if (x, y) in black:
            black.remove((x, y))
            direction += 1
        else:
            black.add((x, y))
            direction -= 1
        (dx, dy) = steps[direction%4]
        x += dx
        y += dy
        position.append((x, y))
    return position
In [305]:
lang_x, lang_y = zip(*[ant_p for ant_p in ant(15000)])
In [306]:
show_scatter(lang_x)
In [307]:
show_scatter(lang_y)
In [308]:
plt.figure(figsize=(8,8))
plt.scatter(lang_x, lang_y, s=1, c='black', marker="p")
plt.show()
In [308]:
 

Comments

comments powered by Disqus