Introduction
I was reading an article about how to graph the "prime parallelograms" in J using verbs, which is based on the second part of video by Numberphile.The idea is to take the nth prime and subtract the number that is obtained by inverting the order of its representation in base 2. For example, the 16th prime is 53 = 1101012, when we reverse it we get 1010112 = 43, so the value is 53 - 43 = 10. If we graph this function, parallelograms appear.
I wondered what would happen if we graph all the numbers instead of only the primes. Also, I wanted to understand where these parallelograms come from.
Additional modules
To graph them in Racket we will need to use some additional modules, so we start with#lang racket
(require math/number-theory)
(require plot)
#;(plot-new-window? #t)
I also like to change how the numbers appear on the axes. It is a technical part, and it almost doesn’t change the result. The code is at the bottom.
Rebuilding the graph
First, we define the function that reverses the order of the bits in the binary representation of the number.(define (binary-reverse x)
(let loop ([x x] [r 0])
(if (zero? x)
r
(loop (quotient x 2) (+ (* 2 r) (remainder x 2))))))
Let's see some examples
(for ([i (in-range 16)])
(define r (binary-reverse i))
(displayln (list i r (number->string i 2) (number->string r 2))))
> (0 0 0 0)
> (1 1 1 1)
> (2 1 10 1)
> (3 3 11 11)
> (4 1 100 1)
> (5 5 101 101)
> (6 3 110 11)
> (7 7 111 111)
> (8 1 1000 1)
> (9 9 1001 1001)
> (10 5 1010 101)
> (11 13 1011 1101)
> (12 3 1100 11)
> (13 11 1101 1011)
> (14 7 1110 111)
> (15 15 1111 1111)
Note: It is possible to define this function using strings, but using strings is always much slower than operating directly with integers.
To simplify the notation, let's call inv this function that reverses the order of the bits, p the function that calculates the nth prime and let's call f the function f(n) = n - inv(n).
Now we can make a list with the points (n, f(p(n))) = (n, p(n) - inv(p(n))) as vectors.
(define prim/original (for/list ([i (in-naturals 1)]
[p (in-list (next-primes 1 10000))])
(vector i (- p (binary-reverse p)))))
Note: We could write our own functions to test primality, but the build-in next-primes and prime? are better.
And we draw it with plot
(plot (points prim/original))
We fix the color and opacity to make the drawing more nice
(plot #:title "Prime Paralelograms (original)"
(points prim/original
#:sym 'dot
#:color 'blue #:alpha .5)
#:x-label "n" #:y-label "f(p(n))"
#:x-min 0)
To save the image we use plot-file
(plot-file #:title "Prime Paralelograms (original)"
(points prim/original
#:sym 'dot
#:color 'blue #:alpha .5)
"prime-paralelograms(original).png"
#:x-label "n" #:y-label "f(p(n))"
#:x-min 0)
In order to compare it with the next graph, we draw it again. But with the labels of the values of the x-axis aligned to the left (so that the edge of the graph is the edge of the image).
(parameterize ([plot-x-tick-label-anchor 'top-right])
(plot #:title "Prime Paralelograms (original)"
(points prim/original
#:sym 'dot
#:color 'blue #:alpha .5)
#:x-label "n" #:y-label "f(p(n))"
#:x-min 0))
Now we change what we graph. We change the x part of the points, instead of using in which position of the list of prime it is, we use the value of the prime, that is, instead of drawing (n, f(p(n)) we draw (p(n), f(p(n)) = (p, f(p)).
(define prim/new (for/list ([i (in-naturals 1)]
[p (in-list (next-primes 1 10000))])
(vector p (- p (binary-reverse p)))))
(parameterize ([plot-x-tick-label-anchor 'top-right])
(plot #:title "Prime Paralelograms"
(points prim/new
#:sym 'dot
#:color 'blue #:alpha .5)
#:x-label "p(n)" #:y-label "f(p(n))"
#:x-min 0))
We get a very similar drawing if we ignore the scale of the x-axis. The parallelograms are cut at slightly different positions, but they are very similar.
What happens is that if we draw the nth prime we get almost a line.
(define prim/scale (for/list ([i (in-naturals 1)]
[p (in-list (next-primes 1 10000))])
(vector i p)))
(plot #:title "Primes"
(points prim/scale
#:sym 'dot)
#:x-label "n" #:y-label "p(n) = nth prime"
#:x-min 0)
The inverse is the function that counts the number of primes that in general is called pi(n) and its slope is slowly changing, more or less like 1/log(n), but log(n) is a function that grows very slowly. So the change between the graphs is almost a linear transformation. Therefore, when we use it to change the x-axis in the previous drawings we see that the shape of the parallelograms changes very little. (We should be able to notice more distortion in the first parallelograms.)
Comparing with all numbers
Now let's see how the graph looks when we use all the numbers instead of just the primes. To be able to compare the graph better, we fix the same range for the y-axis in all the graphs. Also, I like to choose a chart size so that the line y=x has approximately 45°. (In Excel you make a graph and then with the mouse you fix all the details of the axes and scales. Here you have to put all that fixes in the program so that the scales are exactly as you want.)(define all (for/list ([i (in-range 1 131072)])
(vector i (- i (binary-reverse i)))))
(plot #:title "All Paralelograms"
(points all
#:sym 'dot
#:color 'black #:alpha .1)
#:x-label "n" #:y-label "f(n)"
#:y-min -65536 #:y-max 131072
#:x-min 0
#:width 400 #:height 600)
Let’s repeat one of the previous graph, but on these new scales. We will create the list of primes in a slightly different way.
(define prim (for/list ([i (in-range 1 131072)]
#:when (prime? i))
(vector i (- i (binary-reverse i)))))
(plot #:title "Prime Paralelograms"
(points prim
#:sym 'dot
#:color 'blue #:alpha .5)
#:x-label "p" #:y-label "f(p)"
#:y-min -65536 #:y-max 131072
#:x-min 0
#:width 400 #:height 600)
And we superimpose the last two graphs
(plot #:title "All vs Prime Paralelograms"
(list (points all
#:sym 'dot
#:color 'black #:alpha .1)
(points prim
#:sym 'dot
#:color 'blue #:alpha .5))
#:x-label "n" #:y-label "f(n)"
#:y-min -65536 #:y-max 131072
#:x-min 0
#:width 400 #:height 600)
What we see is that the prime parallelograms occupy the bottom half of the parallelograms that are formed when we graph all the numbers.
Even and odd
Almost all primes are odd. As they are odd, when you invert the order of the bits in binary you get a number with the same number of figures in binary. So the result has about the same size. You can bound the result and formalize this idea, but it is nicer to compare in a graphic what happens when we replace prime numbers with odd numbers.(define odd (for/list ([i (in-range 1 131072 2)])
(vector i (- i (binary-reverse i)))))
(plot #:title "All vs Odd Paralelograms"
(list (points all
#:sym 'dot
#:color 'black #:alpha .1)
(points odd
#:sym 'dot
#:color 'red #:alpha .1))
#:x-label "n" #:y-label "f(n)"
#:y-min -65536 #:y-max 131072
#:x-min 0
#:width 400 #:height 600)
We see that the odd number occupy essentially the same parallelograms, although there are fewer gaps.
On the other hand, the even numbers end in 0 in binary, so by reversing the order of their bits in binary they lose at least one figure, so the result is smaller. When we draw them we get this graph.
(define even (for/list ([i (in-range 2 131072 2)])
(vector i (- i (binary-reverse i)))))
(plot #:title "All vs Even Paralelograms"
(list (points all
#:sym 'dot
#:color 'black #:alpha .1)
(points even
#:sym 'dot
#:color 'green #:alpha .1))
#:x-label "n" #:y-label "f(n)"
#:y-min -65536 #:y-max 131072
#:x-min 0
#:width 400 #:height 600)
We see that the even numbers essentially occupy the upper halves, that is, they do not overlap with what is painted by the primes. So the holes in the prime parallelograms are due to non-prime odd numbers.
Some ideas to try
- What happens in base 3? Do we get the obvious generalization? In base 4, do we get flags?
- Redraw to see only the first parallelograms. The change of scale is not so linear at the beginning, so when drawing the parallelograms of the original article (n, f(p(n)) they should be more distorted.
- Redraw everything with the x-axis in logarithmic scale, so that all parallelograms have the same width and you can see the first ones easier. I don't know what scale to use for the y-axis. Maybe you have to draw each parallelogram separately?
- The prime numbers are not random, but they look quite random. It might be interesting to generate a list of fake-primes, filtering the odd numbers with a distribution similar to that of the primes and see how the graph looks.
Changing the format of the labels
I like to change how the numbers appear on the axes, I don't like that 100000 appears as 105. This part is technical and almost does not change the result, so I will not explain the details. (Seems like a good idea for a PR as an additional option for linear-ticks.)(require plot/utils)
(define ((linear-ticks-format/no-sc) min max pre-ticks)
(define digits (digits-for-range min max))
(map (lambda (pt)
(real->plot-label (pre-tick-value pt) digits #f))
pre-ticks))
(define (linear-ticks/no-sc) (ticks (linear-ticks-layout) (linear-ticks-format/no-sc)))
(plot-x-ticks (linear-ticks/no-sc)) ; almost global change
(plot-y-ticks (linear-ticks/no-sc)) ; almost global change
Full code
#lang racket(require math/number-theory)
(require plot)
#;(plot-new-window? #t)
(require plot/utils)
(define ((linear-ticks-format/no-sc) min max pre-ticks)
(define digits (digits-for-range min max))
(map (lambda (pt)
(real->plot-label (pre-tick-value pt) digits #f))
pre-ticks))
(define (linear-ticks/no-sc) (ticks (linear-ticks-layout) (linear-ticks-format/no-sc)))
(plot-x-ticks (linear-ticks/no-sc)) ; almost global change
(plot-y-ticks (linear-ticks/no-sc)) ; almost global change
(define (binary-reverse x)
(let loop ([x x] [r 0])
(if (zero? x)
r
(loop (quotient x 2) (+ (* 2 r) (remainder x 2))))))
(for ([i (in-range 16)])
(define r (binary-reverse i))
(displayln (list i r (number->string i 2) (number->string r 2))))
#;(for ([i (in-range 128)]
#:when (prime? i))
(define r (binary-reverse i))
(displayln (list i r (number->string i 2) (number->string r 2))))
;ORIGINAL
(define prim/original (for/list ([i (in-naturals 1)]
[p (in-list (next-primes 1 10000))])
(vector i (- p (binary-reverse p)))))
(plot (points prim/original))
(plot #:title "Prime Paralelograms (original)"
(points prim/original
#:sym 'dot
#:color 'blue #:alpha .5)
#:x-label "n" #:y-label "f(p(n))"
#:x-min 0)
(plot-file #:title "Prime Paralelograms (original)"
(points prim/original
#:sym 'dot
#:color 'blue #:alpha .5)
"prime-paralelograms(original).png"
#:x-label "n" #:y-label "f(p(n))"
#:x-min 0)
(parameterize ([plot-x-tick-label-anchor 'top-right])
(plot #:title "Prime Paralelograms (original)"
(points prim/original
#:sym 'dot
#:color 'blue #:alpha .5)
#:x-label "n" #:y-label "f(p(n))"
#:x-min 0))
;NEW
(define prim/new (for/list ([i (in-naturals 1)]
[p (in-list (next-primes 1 10000))])
(vector p (- p (binary-reverse p)))))
(parameterize ([plot-x-tick-label-anchor 'top-right])
(plot #:title "Prime Paralelograms"
(points prim/new
#:sym 'dot
#:color 'blue #:alpha .5)
#:x-label "p(n)" #:y-label "f(p(n))"
#:x-min 0))
;SCALE
(define prim/scale (for/list ([i (in-naturals 1)]
[p (in-list (next-primes 1 10000))])
(vector i p)))
(plot #:title "Primes"
(points prim/scale
#:sym 'dot)
#:x-label "n" #:y-label "p(n) = nth prime"
#:x-min 0)
;All
(define all (for/list ([i (in-range 1 131072)])
(vector i (- i (binary-reverse i)))))
(plot #:title "All Paralelograms"
(points all
#:sym 'dot
#:color 'black #:alpha .1)
#:x-label "n" #:y-label "f(n)"
#:y-min -65536 #:y-max 131072
#:x-min 0
#:width 400 #:height 600)
(define prim (for/list ([i (in-range 1 131072)]
#:when (prime? i))
(vector i (- i (binary-reverse i)))))
(plot #:title "Prime Paralelograms"
(points prim
#:sym 'dot
#:color 'blue #:alpha .5)
#:x-label "p" #:y-label "f(p)"
#:y-min -65536 #:y-max 131072
#:x-min 0
#:width 400 #:height 600)
(plot #:title "All vs Prime Paralelograms"
(list (points all
#:sym 'dot
#:color 'black #:alpha .1)
(points prim
#:sym 'dot
#:color 'blue #:alpha .5))
#:x-label "n" #:y-label "f(n)"
#:y-min -65536 #:y-max 131072
#:x-min 0
#:width 400 #:height 600)
; ODD and EVEN
(define odd (for/list ([i (in-range 1 131072 2)])
(vector i (- i (binary-reverse i)))))
(plot #:title "All vs Odd Paralelograms"
(list (points all
#:sym 'dot
#:color 'black #:alpha .1)
(points odd
#:sym 'dot
#:color 'red #:alpha .1))
#:x-label "n" #:y-label "f(n)"
#:y-min -65536 #:y-max 131072
#:x-min 0
#:width 400 #:height 600)
(define even (for/list ([i (in-range 2 131072 2)])
(vector i (- i (binary-reverse i)))))
(plot #:title "All vs Even Paralelograms"
(list (points all
#:sym 'dot
#:color 'black #:alpha .1)
(points even
#:sym 'dot
#:color 'green #:alpha .1))
#:x-label "n" #:y-label "f(n)"
#:y-min -65536 #:y-max 131072
#:x-min 0
#:width 400 #:height 600)
