Pour Le Science and the anti-Sum-Product Problem

In March 2021, I got an out-of-the-blue email from OEIS editor Michel Marcus which totally delighted me. He wrote:

This afternoon I went to the library.
And I was browsing “Pour La Science” the French version of the Scientific American.
And here is what I saw.

I like the mysterious tone.

He included this photo of an excerpt of an article, which on the first line mentioned an OEIS sequence of mine (with my name slightly misspelled):

Michel Marcus’s photo of an article in Pour La Science. (Translated below)

The article, translated

Michel had sent me a snippet of the article “Addition et multiplication, des tables qui intriguent” by French Computer Science professor Jean-Paul Delahaye. I read French quite clumsily, so I sent it to my friend Alec, and he volunteered a translation.

Sloane [verb not pictured] of my sequence A337655, Peter Kagey, who must have been informed of changes to the encyclopedia, proposed this second way of accounting for addition and multiplication. This has yielded a new entry in the encyclopedia, the following A337946:

1, 3, 7, 12, 22, 30, 47, 61, 85, 113, 126, 177, 193, 246, 279, 321, 341, 428, 499, 571, 616, 686, 754, 854, 975, 1052, 1150, 1317, 1376, 1457, …

Impossible Little Tables

Let’s consider the fourth question. The response is negative this time: no, it’s not possible that the tables of addition and multiplication have simultaneously only a few distinct values. This is a result of Paul Erdős and Endre Szemerédi from 1983, which we state precisely here:

“There exist two constants \(C>0\) and \(e > 0\) such that, for the set \(E\) of real numbers of […]

Translation by Alec Jones

A Sum-Product Problem

The article—or at least the part that I could make sense of—is about the opposite of Erdos’s Sum-Product Problem. This article explores questions about the greatest number of different values in the corresponding addition and multiplication tables for different sequences of positive integers.

Addition table for the first six terms of A005282, the lexicographically earliest infinite sequence with the maximum-sum property.
Multiplication table for the first six terms of A066720, the lexicographically earliest infinite sequence with the maximum-property property.

Jean-Paul’s sequence, A337655, is the lexicographically earliest (greedy) sequence such that for all \(n\) both the addition and multiplication tables of the first \(n\) terms have \(\binom{n+1}{2}\) distinct terms (the greatest number).

Addition table for A337655.
Multiplication table for A337655.

The addition and multiplication tables share some values. For example, they both have a \(2\), a \(4\), a \(7\), a \(10\), and so on. My sequence, A337946, is similar to Jean-Paul’s but with the additional restriction that the addition and multiplication tables cannot have any values in common.

Addition table for A337946.
Multiplication table for A337946.

Some related sequences

I’ve recently added some new sequences to the OEIS related to these sequences.

One way of thinking of A066720 is that it’s the lexicographically earliest infinite sequence \(S = \(a_i\)_{i=1}^{\infty}\) such that for all \(n\) the set \(S_n = \{a_i\}_{i=1}^{n}\) consisting of the first \(n\) terms of \(S\) has the property that \(S_n \times S_n = \{xy \mid x, y \in S_n\}\) has exactly \(\binom{n+1}{2}\) elements—the most possible.

A347498: minimize the largest term

If instead of minimizing the lexicographic order, we instead minimize the size of the largest element, we get OEIS sequence A347498: for each \(n\), find the least \(m\) such that there exists a subset \(T_n \subseteq \{1, 2, \dots, m\}\) with \(n\) elements such that \(|T_n \times T_n| = \binom{n+1}{2}\). \(A347498(n)\) records the value of the \(m\)s.

For example, \(A347498(8) = 11\) as illustrated by the table below.

Multiplication table illustrating \(A347498(8) = 11\).

A347499: Examples with minimized largest element

We wish to record examples of the sequences described in the above section. There are \(A348481(n)\) subsets of \(\{1,2,…,A347499(n)\}\) with the distinct product property, but we record the lexicographically earliest one in OEIS sequence A347499.

 n | Distinct product subset of {1,2,...,A347499(n)}
 1 | {1}
 2 | {1, 2}
 3 | {1, 2, 3}
 4 | {1, 2, 3, 5}
 5 | {1, 3, 4, 5,  6}
 6 | {1, 3, 4, 5,  6,  7}
 7 | {1, 2, 5, 6,  7,  8,  9}
 8 | {1, 2, 5, 6,  7,  8,  9, 11}
 9 | {1, 2, 5, 6,  7,  8,  9, 11, 13}
10 | {1, 2, 5, 7,  8,  9, 11, 12, 13, 15}
11 | {1, 2, 5, 7,  8,  9, 11, 12, 13, 15, 17}
12 | {1, 2, 5, 7,  8,  9, 11, 12, 13, 15, 17, 19}
13 | {1, 5, 6, 7,  9, 11, 13, 14, 15, 16, 17, 19, 20}
14 | {1, 2, 5, 7, 11, 12, 13, 16, 17, 18, 19, 20, 21, 23} 

For all of the known values, the lexicographically earliest subset begins with \(1\). Is this true in general?

A347570: Sums with more terms

Perhaps instead of asking when all sums \(x + y\) are distinct for \(x \leq y\), we want to know when all sums \(x_1 + x_2 + \dots + x_n\) are distinct for \(x_1 \leq x_2\leq \dots \leq x_n\). This family of sequences are sometimes called \(B_n\) sequences and have been studied by the likes of Richard Guy.

I’ve recorded the lexicographically earliest \(B_n\) sequences in OEIS sequence A347570.

n\k | 1  2   3   4    5     6     7      8
  1 | 1, 2,  3,  4,   5,    6,    7,     8, ...
  2 | 1, 2,  4,  8,  13,   21,   31,    45, ...
  3 | 1, 2,  5, 14,  33,   72,  125,   219, ...
  4 | 1, 2,  6, 22,  56,  154,  369,   857, ...
  5 | 1, 2,  7, 32, 109,  367,  927,  2287, ...
  6 | 1, 2,  8, 44, 155,  669, 2215,  6877, ...
  7 | 1, 2,  9, 58, 257, 1154, 4182, 14181, ...
  8 | 1, 2, 10, 74, 334, 1823, 8044, 28297, ... 

I ask about this family of sequences in a recent question on Code Golf Stack Exchange.

Other ideas?

Did this post spark any related ideas? Let me know about them on Twitter @PeterKagey—I’d love to hear about them!

Zimin Words and Bifixes

One of the earliest contributions to the On-Line Encyclopedia of Integer Sequences (OEIS) was a family sequences counting the number of words that begin (or don’t begin) with a palindrome:

  • Let \(f_k(n)\) be the number of strings of length \(n\) over a \(k\)-letter alphabet that begin with a nontrivial palindrome” for various values of \(k\).
  • Let \(g_k(n)\) be the number of strings of length n over a \(k\)-letter alphabet that do not begin with a nontrivial palindrome.
  • Number of binary strings of length \(n\) that begin with an odd-length palindrome. (A254128)

(If I had known better, I would have published fewer sequences in favor of a table, and I would have requested contiguous blocks of A-numbers.)

I must have written some Python code to compute some small terms of this sequence, and I knew that \(g_k(n) = k^n – f_k(n)\), but I remember being at in my friend Q’s bedroom when the recursion hit me for \(f_k(n)\): $$f_k(n) = kf_k(n-1) + k^{\lceil n/2 \rceil} – f_k\big(\lceil \frac n 2 \rceil \big)$$

“Bifix-free” words

One sequence that I didn’t add to the OEIS was the “Number of binary strings of length n that begin with an even-length palindrome”—that’s because this was already in the Encyclopedia under a different name:

A094536: Number of binary words of length n that are not “bifix-free”.

0, 0, 2, 4, 10, 20, 44, 88, 182, 364, 740, 1480, 2980, 5960, …

A “bifix” is a shared prefix and suffix, so a “bifix-free” word is one such that all prefixes are different from all suffixes. More concretely, if the word is \(\alpha_1\alpha_2 \dots \alpha_n\), then \((\alpha_1, \alpha_2, \dots, \alpha_k) \neq (\alpha_{n-k+1},\alpha_{n-k+2},\dots,\alpha_n)\) for all \(k \geq 1\).

The reason why the number of binary words of length \(n\) that begin with an even length palindrome is equal to the number of binary words of length \(n\) that have a bifix is because we have a bijection between the two sets. In particular, find the shortest palindromic prefix, cut it in half, and stick the first half at the end of the word, backward. I’ve asked for a better bijection on Math Stack Exchange, so if you have any ideas, please share them with me!

In 2019–2020, Daniel Gabric, Jeffrey Shallit wrote a paper closely related to this called Borders, Palindrome Prefixes, and Square Prefixes.

Zimin words

A Zimin word can be defined recursively, but I think it’s most suggestive to see some examples:

  • \(Z_1 = A\)
  • \(Z_2 = ABA\)
  • \(Z_3 = ABACABA\)
  • \(Z_n = Z_{n-1} X Z_{n-1}\)

All Zimin words \(Z_n\) are examples of “unavoidable patterns”, because every sufficiently long string with letters in any finite alphabet contains a substring that matches the \(Z_n\) pattern.

For example the word \(0100010010111000100111000111001\) contains a substring that matches the Zimin word \(Z_3\). Namely, let \(A = 100\), \(B = 0\), and \(C = 1011\), visualized here with each \(A\) emboldened: \( 0(\mathbf{100}\,0\,\mathbf{100}\,1011\,\mathbf{100}\,0\,\mathbf{100})111000111001\).

Binary solo from Bret at 1:26.

I’ve written a Ruby script that generates a random string of length 29 and uses a regular expression to find the first instance of a substring matching the pattern \(Z_3 = ABACABA\). You can run it on TIO, the impressive (and free!) tool from Dennis Mitchell.

# Randomly generates a binary string of length 29.
random_string = 29.times.map { [0,1].sample }.join("")
p random_string
# Finds the first Zimin word ABACABA
p random_string.scan(/(.+)(.+)\1(.+)\1\2\1/)[0]
# Pattern:             A   B   A C   A B A

Why 29? Because all binary words of length 29 contain the pattern \(Z_3 = ABACABA\). However, Joshua Cooper and Danny Rorabaugh’s paper provides 48 words of length 28 that avoid that pattern (these and their reversals):



The Zimin Word \(Z_2 = ABA\) and Bifixes

The number of Zimin words of length \(n\) that match the pattern ABA is equal to the number of of words that begin with an odd-length palindrome. Analogously, the number of words with a bifix is equal to the number of words that begin with an even-length palindrome. The number of these agree when \(n\) is odd.

I’ve added OEIS sequences A342510A342512 which relate to how numbers viewed as binary strings avoid—or fail to avoid—Zimin words. I asked users to implement this on Code Golf Stack Exchange.

My Favorite Sequences: A263135

This is the fourth in my installment of My Favorite Sequences. This post discusses sequence A263135 which counts penny-to-penny connections among \(n\) pennies on the vertices of a hexagonal grid. I published this sequence in October 2015 when I was thinking about hexagonal-grid analogs to the “Not Equal” grid. The square-grid analog of this sequence is A123663.

A263135: Placing Pennies

The sequences A047932 and A263135 are about placing pennies on a hexagonal grid in such a way that maximizes the number of penny-to-penny contacts, which occurs when you place the pennies in a spiral. A047932, counts the contacts when the pennies are placed on the faces of the grid; A263135 counts the contacts with the pennies placed on the vertices.

Pattern of placing pennies in A047932.
Pattern of placing pennies in A263135.

While spiral shapes maximize the number of penny-to-penny contacts, there are sometimes non-spiral shapes that have the same number of contacts. For example, in the case of the square grid, there are \(A100092(n)\) such ways to lay down \(n\) pennies on the square grid with the maximum number of connections. Problem 108 in my Open Problems Collection asks about generalizing this OEIS sequence to other settings such as the hexagonal grid.

A047932(11) = 21
A263135(22) = 27

Comparing contacts

Notice that the “face” pennies in A047932 can have a maximum of six neighbors, while the “vertex” pennies in A263135 can have a maximum of three. In the limit, most pennies are “interior” pennies with the maximum number of contacts, so \(A047932(n) \sim 3n\) and \(A263135(n) \sim \frac32n\).

Looking at the comparative growth rates, it is natural to ask how the number of connections of \(n\) face pennies compares to the number of connections of \(2n\) vertex pennies. In October 2015 I made a conjecture on the OEIS that this difference grew like sequence A216256.

Conjecture: For \(n > 0\), \[A263135(2n) – A047932(n) = \lceil\sqrt{3n – 3/4} – 1/2\rceil = A216256(n).\]

I believe that the sequence A216256 on the right hand side appears to be the same as the sequence “n appears \(\displaystyle\left\lfloor \frac{2n+1}{3} \right\rfloor\) times,” but I’d have to crack open my Concrete Mathematics book to prove it.

This is Problem 20 in my Open Problem Collection, and I’ve placed a small, $5 bounty on solving this conjecture—so if you have an idea of how to prove this, let me know in exchange for a latte! I’ve asked about this in my Math Stack Exchange question Circle-to-circle contacts on the hexagonal grid—so feel free to answer there or let me know on Twitter, @PeterKagey.

My Favorite Sequences: “Not Equal” Grid

This is the third installment in a recurring series, My Favorite Sequences. This post discusses OEIS sequence A278299, a sequence that took over two years to compute enough terms to add to the OEIS with confidence that it was distinct.

This sequence is discussed in Problem #23 of my Open Problems Collection, which asks for the smallest polyomino (by number of cells) whose cells you can color with \(n\) different colors such that any two different colors are adjacent somewhere in the polyomino. As illustrated below, when there are \(n=5\) colors (say, green, brown, blue, purple, and magenta) there is a \(13\)-cell polyomino which has a green cell adjacent to a blue cell and a purple cell adjacent to a brown cell and so on for every color combination. This is the smallest polyomino with the \(5\)-coloring property.

Five colors of blocks, where any two different colors of blocks are adjacent somewhere in the polyomino.

The Genesis: Unequal Chains

The summer after my third undergraduate year, I decided to switch my major to Math and still try to graduate on time. Due to degree requirements, I had to go back and take some lower-division classes that I was a bit over-prepared for. One of these classes—and surely my favorite—was Bill Bogley‘s linear algebra class, where I half-way paid attention and half-way mused about other things.

Bill wrote something simple on the board that sparked inspiration for me: $$a \neq b \neq c \neq a.$$ He wrote this to indicate that \(a\), \(b\), and \(c\) were all distinct, and this got me thinking: if we have to write a string of four variables in order to say that three variables are distinct, how many would we have to write down to say that four variables were distinct? It turns out that \(8\) will do the trick, with one redundancy: $$a\neq b \neq c \neq d \neq b \color{red}{\neq} c \neq a.$$ Five variables? \(11\): $$a_1 \neq a_2 \neq a_3 \neq a_4 \neq a_5 \neq a_3 \neq a_1 \neq a_4 \neq a_2 \neq a_5.$$ What about \(n\) variables?

My colleague and the then-President of the OSU Math Club, Tommy Pitts, made quick work of this problem. He pointed out that “not equal” is a symmetric, non-transitive, non-reflexive relation. This means that we can model this with a complete graph on \(n\) vertices, where each edge is a relation. Then the number of variables needed in the expression is the number of edges in the complete graph, plus the minimum number of Eulerian paths that we can split the graph into. Searching for this in the OEIS yields sequence A053439. $$A053439^*(n) = \begin{cases} \binom{n}{2} + 1 & n \text{ is odd} \\ \binom{n}{2} + \frac n 2 & n \text{ is even}\end{cases}$$

A Generalization: Unequal Chainmail

This was around 2014, at which time I was writing letters to my friend Alec Jones whenever I—rather frequently!—stumbled upon a new math problem that interested me. In the exchange of letters, he suggested a 2D version of this puzzle. Write the \(n\) variables in the square grid, and say that two variables are unequal if they’re adjacent.

While Tommy solved the 1D version of the problem quickly, the 2D version was much more stubborn! However we were able to make some progress. We found some upper bounds (e.g. the 1D solution) and some lower bounds, and we were able to prove that some small configurations were optimal. Finally, in November 2016, we had ten terms: enough to prove that this sequence was not in the OEIS. We added it as A278299.

\(a(n)\) is the tile count of the smallest polyomino with an \(n\)-coloring such that every color is adjacent to every other distinct color at least once.

OEIS sequence A278299.

(In May 2019, Alec’s student Ryan Lee found the \(11\)th term: \(A278299(11) = 34\). \(A278299(12)\) is still unknown.)

A screenshot from my game illustrating the largest known term: \(A278299(14) = 56\). Every number is connected to every other number. The red edges refer to redundant connections.

We found these terms by establishing some lower bounds (as explained below) and then implementing a Javascript game (which you can play here) with a Ruby on Rails backend to allow people to submit their hand-crafted attempts. Each solution was constructive proof of an upper bound, so when a user submitted a solution that matched the lower bound, we were able to confirm that term of the sequence.

(One heuristic for making minimal configurations is to start with the construction in OEIS sequence A260643 and add cells as necessary in an ad hoc fashion.)

Lower bounds

There are a few different ways of proving lower bounds.

  • We know that there needs to be at least \(\binom{n}{2}\) relations, one between each pair of variables. OEIS sequence A123663 gives the “number of shared edges in a spiral of n unit squares,” which can be used to compute a lower bound: $$A039823(n) = \left\lceil \frac{n^2+n+2}{4}\right\rceil$$
  • Every number needs to be in contact with at least (n-1) other numbers, and each occurrence can be in contact with at most (4) others. So each number needs to occur at least \(\lceil \frac{n-1}{4}\rceil\) times, for a total of \(n\lceil \frac{n-1}{4}\rceil\) occurrences. This bound is usually weaker than the above bound.
  • For the cases of \(n = 5\) and \(n=9\), the lower bounds were proved using ad hoc methods, by looking at how many cells would need to have a given number of neighbors.

Upper Bounds

Besides the upper bound that comes from the 1-dimensional version of the problem, that only upper bounds that I know of come from hand-crafted submissions on my Javascript game on my website.

Do you have any ideas for an explicit and efficient algorithm for constructing such solutions? If so, let me know on Twitter @PeterKagey.


The lower and upper bounds show that this is asymptotically bounded between \(n^2/4\) and \(n^2/2\). It’s possible that this doesn’t have a limit at all, but it would be interesting to bound the liminf and limsup further. My intuition is that \(n^2/4\) is the right answer, can you prove or disprove this?


  • We could play this game on the triangular grid, or in the 3-dimensional cubic grid. Do you have ideas of other graphs that you could do this on?
  • This game came from Tommy’s analysis of looking at “not equal to” as a symmetric, non-reflexive, non-transitive relation. Can you do a similar analysis on other kinds of relations?
  • Is there a good way of defining what it means for two solutions to be the same? For a given number of variables, how many essentially different solutions exist? (Related: Open problem #108.)
  • What if we think of left-right connections as being different from up-down connections, and want both? Or what if we want each variable \(x\) to be neighbors with another \(x\)?

If you have ideas about these questions or any questions of your own, please share them with me by leaving a comment or letting me know on Twitter, @PeterKagey!

My Favorite Sequences: A289523

This the second post in a recurring series, My Favorite Sequences. If you like this sort of thing, check out the Integer Sequence Review from The Aperiodical!

A289523: Packing Circles of Increasing Area

A plot of the circles in A289523.

In July 2017, I added a mathematically-silly-but-visually-fun sequence, A289523. The sequence works like this: place a circle of area \(\pi\) centered at \((1,1)\), then place a circle of area \(2\pi\) centered at \((2,a(2))\) where \(a(2) = 4\), the least positive integer such that the circle does not overlap with the first circle. Next, place a circle of area \(3\pi\) centered at \((3,a(3))\) where \(a(3) = 7\) is the least positive integer such that the circle does not overlap with the first two circle. Continue this pattern ad infinitum, creating the earliest infinite sequence of positive integers such that no two circles overlap with any others, and a circle centered at \((k, a(k))\) has area \(k\pi\).

I haven’t done much mathematical analysis on this problem, but it would be interesting to see if it’s possible to compute (or put some bounds on) the packing density of the convex hull of the circles. Also, a glance at a plot of the points suggests that the sequence is bounded above by a linear function—is this the case?

The sequence begins 1, 4, 7, 1, 11, 16, 5, 21, 27, 34, 10, 1, 41, 17, 49, 25, 57, 6, 33, 66, 43, 14, ....

Finding an upper bound

The scatter plot of A289523 suggests that the centers of the circles have a linear upper bound. This is to be expected! The areas of the circles increase linearly, and the packing density is (presumably) nonzero.

What is the slope of the upper bound? And what is the packing density of these circles in the limit?

A scatter plot of A289523.

Related Construction

At the end of March, I posted a related puzzle, “Placing Circles Along a Square Spiral”, on Code Golf Stack Exchange. For the post, I made a few animated GIFs that explain the construction and tweeted about them.

Impressively, Code Golf Stack Exchange users tsh, Arnauld, and A username each wrote (deliberately terse) Javascript code that computes the placement of these circles.

In fact, they compute something strictly harder! In the challenge, after laying down all of these circles (in blue), the challenge instructed them to go back to the start and greedily fill the gaps with (red) circles of increasing area. Next, they laid down a third (yellow) generation in the same fashion, and fourth (cyan) generation, and so on.

An animation illustrating successive generations of circles on a square spiral.

Related questions

  • What is the packing density of the first (blue) generation?
  • What is the packing density of the \(k\)-th generation?
  • How many “steps” away from the origin is the smallest circle in the \(k\)-th generation?
  • Do an infinite number of blue circles touch? Do an infinite number of any circles touch? Which ones?
  • How far can a circle be from its neighbors? Which circles are maximally far from their neighbors?
  • How does this work if the path the circles follow is not the spiral? Can different paths have significantly different packing densities?

If you have thoughts or ideas about any of this—or if you just want to make animated GIFs together—leave a comment or let me know on Twitter, @PeterKagey!

My Favorite Sequences: A261865

This is the first installment in a new series, “My Favorite Sequences”. In this series, I will write about sequences from the On-Line Encyclopedia of Integer Sequences that I’ve authored or spent a lot of time thinking about.

I’ve been contributing to the On-Line Encyclopedia of Integer Sequences since I was an undergraduate. In December 2013, I submitted sequence A233421 based on problem A2 from the 2013 Putnam Exam—which is itself based on “Ron Graham’s Sequence” (A006255)—a surprising bijection from the natural numbers to the non-primes. As of today, I’ve authored over 475 sequences based on puzzles that I’ve heard about and problems that I’ve dreamed up.

A261865: Multiples of square roots

(This problem is closely related to Problem 13 in my Open Problems Collection.)

In September 2015, I submitted sequence A261865:

\(A261865(n)\) is the least integer \(k\) such that some multiple of \(\sqrt k\) falls in the interval \((n, n+1)\).

An illustration of the first dozen terms of A261865

For example, \(A261865(3) = 3\) because there is no multiple of \(\sqrt 1\) in \((3,4)\) (since \(3 \sqrt{1} \leq 3\) and \(4 \sqrt{1} \geq 4\)); there is no multiple of \(\sqrt{2}\) in \((3,4)\) (since \(2 \sqrt{2} \leq 3\) and \(3 \sqrt 2 \geq 4\)); but there is a multiple of \(\sqrt 3\) in \((3,4)\), namely \(2\sqrt 3\).

As indicated in the picture, the sequence begins $$\color{blue}{ 2,2,3,2,2},\color{red}{3},\color{blue}{2,2,2},\color{red}{3},\color{blue}{2,2},\color{red}{3},\color{blue}{2,2,2},\color{red}{3},\color{blue}{2,2},\color{red}{3},\color{blue}{2,2},\color{magenta}{7},\dots.$$

A scatterplot of \(A261865(n)\). Notice the records at \(A261865(184)=38\) and \(A261865(8091)=43\).

A conjecture about density

As the example illustrates, \(1\) does not appear in the sequence. And almost by definition, asymptotically \(1/\sqrt 2\) of the values are \(2\)s.

Let’s denote the asymptotic density of terms that are equal to \(n\) by \(d_n\). It’s easy to check that \(d_1 = 0\), (because multiples of \(\sqrt 1\) are never between any integers) and \(d_2 = 1/\sqrt 2\), because multiples of \(\sqrt 2\) are always inserted. I conjecture in Problem 13 of my Open Problem Collection that $$a_n = \begin{cases}\displaystyle\frac{1}{\sqrt n}\left(1 – \sum_{i=1}^{n-1} a_i\right) & n \text{ is squarefree}\\[5mm] 0 & \text{otherwise}\end{cases}$$

If this conjecture is true, then the following table gives approximate densities.

\(1\)\(d_1 = 0\%\)
\(2\)\(d_2 = 70.7\%\)
\(3\)\(d_3 = 16.9\%\)
\(4\)\(d_4 = 0\%\)
\(5\)\(d_5 = 5.54\%\)
\(6\)\(d_6 = 2.79\%\)
\(7\)\(d_7 = 1.53\%\)
\(10\)\(d_{10} = 0.797\%\)
\(11\)\(d_{11} = 0.519\% \)
\(399\)\(d_{399} = 3.53 \times 10^{-11} \%\)

This was computed with the Mathematica code:

d[i_] := (d[i] = If[
  N[(1 - Sum[d[j], {j, 2, i - 1}])/Sqrt[i], 50], 

Finding Large Values

I’m interested in values of \(n\) such that \(A261865(n)\) is large, and I reckon that there are clever ways to construct these, perhaps by looking at some Diophantine approximations of \(\sqrt{2}, \sqrt{3}, \sqrt{5}, \sqrt{6}, \dots\). In February, I posted a challenge on Code Golf Stack Exchange to have folks compete in writing programs that can quickly find large values of \(A261865(n)\).

Impressively, Noodle9’s C++ program won the challenge. In under a minute, this program found that the input \(n=1001313673399\) makes \(A261865\) particularly large: \(A261865(1001313673399) = 399\). Within the time limit, no other programs could find a value of \(n\) that makes \(A261865(n)\) larger.

\(n\)Order of magnitude\(A261865(n)\)Time
1 \(1 \times 10^{0}\)2(0s)
3 \(3 \times 10^{0}\)3 (0s)
23 \(2.3 \times 10^{1}\)7 (0s)
30 \(3.0 \times 10^{1}\)15 (0s)
184 \(1.84 \times 10^{2}\)38 (0s)
8091 \(8.091 \times 10^{3}\)43 (0s)
16060 \(1.606 \times 10^{4}\)46 (0s)
16907 \(1.691 \times 10^{4}\)58 (0s)
20993 \(2.099 \times 10^{4}\)61 (0s)
26286 \(2.629 \times 10^{4}\)97 (0s)
130375 \(1.304 \times 10^{5}\)118 (0s)
169819 \(1.698 \times 10^{5}\)127 (0s)
2135662 \(2.136 \times 10^{6}\)130 (0s)
2345213 \(2.345 \times 10^{6}\)187 (0s)
46272966 \(4.627 \times 10^{7}\)193 (1s)
222125822 \(2.221 \times 10^{8}\)210 (5.2s)
237941698 \(2.379 \times 10^{8}\)217 (5.7s)
257240414 \(2.572 \times 10^{8}\)227 (6.2s)
1205703469 \(1.206 \times 10^{9}\)267 (31s)
1558293414 \(1.558 \times 10^{9}\)299 (41.8s)
4641799364 \(4.642 \times 10^{9}\)303 (2.1m)
6600656102 \(6.601 \times 10^{9}\)323 (3m)
11145613453 \(1.115 \times 10^{10}\)335 (5.2m)
20641456345 \(2.064 \times 10^{10}\)354 (9.8m)
47964301877 \(4.796 \times 10^{10}\)358 (22.9m)
105991039757 \(1.06 \times 10^{11}\)385 (52m)
119034690206 \(1.19 \times 10^{11}\)397 (59.1m)
734197670865 \(7.342 \times 10^{11}\)455 (6.4h)
931392113477 \(9.314 \times 10^{11}\)501 (8.4h)
1560674332481 \(1.561 \times 10^{12}\)505 (14.2h)
A table of record values as computed by Code Golf Stack Exchange user Neil. The first 16 values agree with Jon E. Schoenfield’s computations that were added to the OEIS in September 2015

Related Ideas

Sequence \(A327953(n)\) counts the number of positive integers \(k\) such that there is some integer \(\alpha^{(n)}_k > 2\) where \(\alpha^{(n)}_k\sqrt{k} \in (n, n+1)\). It appears to grow roughly linearly like \(A327953(n) \sim 1.3n\), but I don’t know how to prove this.

  • Take any function \(f\colon\mathbb N \rightarrow \mathbb R\) that is positive, has positive first derivative, and has negative second derivative. Then, what is the least \(k\) such that some multiple of \(f(k)\) is in \((n,n+1)\)?
  • For example, what is the least integer \(k \geq 3\) such that there is a multiple of \(\ln(k)\) in \((n, n+1)\)?
  • What is the least \(k \in \mathbb N\) such that there exists \(m \in \mathbb N\) with \(k2^{1/m} \in (n,n+1)\)?
  • What is the least \(m \in \mathbb N\) such that there exists \(k \in \mathbb N\) with \(k2^{1/m} \in (n,n+1)\)?
  • A343205 is the auxiliary sequence that gives the value \(m\) such that \(m\sqrt{A261865(n)} \in (n, n+1)\). Does this sequence have an infinite limit inferior?
Scatterplot of A343205, generated in Mathematica. If the main conjecture is true, then this is not bounded below by \(\alpha n\) for any positive value of \(\alpha\).

If you can answer any of these questions, or if you spend time thinking about this, please let me know on Twitter, @PeterKagey!

Richard Guy’s Partition Sequence

Neil Sloane is the founder of the On-Line Encyclopedia of Integer Sequences (OEIS). Every year or so, he gives a talk at Rutgers in which he discusses some of his favorite recent sequences. In 2017, he spent some time talking about a 1971 letter that he got from Richard Guy, and some questions that went along with it. In response to the talk, I investigated the letter and was able to sort out Richard’s 45-year-old idea, and correct and compute some more terms of his sequence.

Richard Guy and his sequences

Richard Guy was a remarkable mathematician who lived to the remarkable age of 103 years, 5 months, and 9 days! His life was filled with friendships and collaborations with many of the giants of recreational math: folks like John Conway, Paul Erdős, Martin Gardner, Donald Knuth, and Neil Sloane. But what I love most about Richard is how much joy and wonder he found in math. (Well, that and his life-long infatuation with his wife Louise.)

Richard guy mountaineering at age 98 with a photo of his late wife, Louise.

[I’m] an amateur [mathematician], I mean I’m not a professional mathematician. I’m an amateur in the more genuine sense of the word in that I love mathematics and I would like everybody in the world to like mathematics.

Richard Guy in Fascinating Mathematical People: Interviews and Memoirs

Richard’s letter to Neil

In January 2017, Neil Sloane gave a talk at Doron Zeilberger’s Experimental Mathematics Seminar, and about six minutes in, Neil discusses a letter that Richard sent to him at Cornell—which was the forwarded to Bell Labs—in June 1971.

Richard Guy’s 1971 letter to Neil Sloane.

When I was working on the book, the 1973 Handbook of Integer Sequences, I would get letters from Richard Guy from all over the world. As he traveled around, he would collect sequences and send them to me.

Neil Sloane, Rutgers Experimental Mathematics Seminar

At 11:30, Neil discusses “sequence I” from Richard’s letter, which he added to the OEIS as sequence A279197:

Number of self-conjugate inseparable solutions of \(X + Y = 2Z\) (integer, disjoint triples from \(\{1,2,3,\dots,3n\}\)).

Neil mentioned in the seminar that he didn’t really know exactly what the definition meant. With some sleuthing and programming, I was able to make sense of the definition, write a Haskell program, correct the 7th term, and extend the sequence by a bit. The solutions for \(A279197(1)\) through \(A279197(10)\) are listed in a file I uploaded to the OEIS, and Fausto A. C. Cariboni was able to extend the sequence even further, submitting terms \(A279197(11)\)–\(A279197(17)\).

How the sequence works.

The idea here is to partition \(\{1,2,3,\dots,3n\}\) into length-3 arithmetic progressions, \(\bigl\{\{X_i,Z_i,Y_i\}\bigr\}_{i=1}^{n}\). And in particular, we want them to be inseparable and self-conjugate.

An inseparable partition is one whose “smallest” subsets are not a solution for a smaller case. For example, if \(n=3\), then the partition \[\bigl\{ \{1,3,5\}, \{2,4,6\}, \{7,8,9\} \bigr\}\] is separable, because if the subset \(\bigl\{ \{1,3,5\}, \{2,4,6\} \bigr\}\) is a solution to the \(n=2\) case.

A self-conjugate partition is one in which swapping each \(i\) with each \(3n+1-i\) gets back to what we started with. For example, \(\bigl\{\{1,3,5\}, \{2,4,6\}\bigr\}\) is self-congugate, because if we replace the \(1\) with a \(6\) and the \(2\) with a \(5\), and the \(i\) with a \(7-i\), then we get the same set: \(\bigl\{\{6,4,2\}, \{5,3,1\} \bigr\}\)

(1,3,5),  (2,7,12), (4,9,14),  (6,8,10),  (11,13,15)
(1,3,5),  (2,8,14), (4,7,10),  (6,9,12),  (11,13,15)
(1,5,9),  (2,3,4),  (6,8,10),  (7,11,15), (12,13,14)
(1,5,9),  (2,4,6),  (3,8,13),  (7,11,15), (10,12,14)
(1,6,11), (2,3,4),  (5,10,15), (7,8,9),   (12,13,14)
(1,6,11), (2,7,12), (3,8,13),  (4,9,14),  (5,10,15)
(1,7,13), (2,4,6),  (3,9,15),  (5,8,11),  (10,12,14)
(1,7,13), (2,8,14), (3,9,15),  (4,5,6),   (10,11,12)
(1,8,15), (2,3,4),  (5,6,7),   (9,10,11), (12,13,14)
(1,8,15), (2,4,6),  (3,5,7),   (9,11,13), (10,12,14)
(1,8,15), (2,4,6),  (3,7,11),  (5,9,13),  (10,12,14)
Each line shows one of the \(A279197(5) = 11\) inseparable, self-conjugate partitions of \(\{1,2,\dots,15\}\).

Generalizing Richard Guy’s idea

Of course, it’s natural to wonder about the separable solutions, or what happens if the self-conjugate restriction is dropped. In exploring these cases, I found four cases already in the OEIS, and I computed five more: A282615A282619.

Self-conjugateA282615A279197 A282616
Non-self-conjugateA282618A282617 A282619
EitherA279199A202705 A104429
Table of sequences counting the ways of partitioning a set into length-3 arithmetic progressions, subject to various restrictions.

Generalizing further

There are lots of other generalizations that might be interesting to explore. Here’s a quick list:

  • Look at partitions of \(\{1,2,\dots,kn\}\) into \(n\) parts, all of which are an arithmetic sequence of length \(k\).
  • Count partitions of \(\{1,2,\dots,n\}\) into any number of parts of (un)equal size in a way that is (non-)self-conjugate and/or (in)separable.
  • Consider at partitions of \(\{1,2,\dots,3n\}\) into \(n\) parts, all of which are an arithmetic sequence of length \(3\), and whose diagram is “non-crossing”, that is, none of the line segments overlap anywhere. (See the 6th and 11th cases in the example for \(A279197(6) = 11\).)

If explore any generalizations of this problem on your own, if you’d like to explore together, or if you have an anecdotes about Richard Guy that you’d like to share, let me know on Twitter!

A π-estimating Twitter bot: Part III

In the final part of this three-part series, I’ll give technical step-by-step instructions for how to wire up our Twitter bot, @BotfonsNeedles, to Docker and deploy it on the free tier of AWS Lambda, so that it can run until the end of time. I’ll also include some tips that I wish I knew when I got started.

If you’d like to make a Twitter bot, but find this guide intimidating, you can fork the repository and follow the README on the GitHub page for my other bot, @oeisTriangles. Or better yet, I would love to set up a call and walk you through it step-by-step! Just let me know on Twitter.

The plan

And in this final part, we will

  • build and run a Docker image so that you can run it on a local server,
  • push the Docker image up to Amazon Elastic Container Registry (ECR),
  • hook up the Docker image to an AWS Lambda function, and
  • configure the function to run on a timer in perpetuity.
My Twitter bot, @RobotWalks, deployed using the techniques in this article.
My Twitter bot, @xorTriangles, which posts every 12 hours via AWS Lambda.

Turn it into a Docker image

Next, we’ll package this up in a Docker image so that AWS Lambda has everything that it needs to run this function. Begin by downloading Docker if you don’t already have it installed.

Next, we’re going to add a new file called Dockerfile from an AWS base image for Python, which will look like this.

FROM amazon/aws-lambda-python:3.8

RUN pip install tweepy
RUN pip install Pillow -t .

COPY random_needle.py ./
COPY needle_drawer.py ./
COPY secrets.py ./
COPY twitter_accessor.py ./
COPY tweet_builder.py ./
COPY app.py ./

CMD ["app.handler"]
  • The FROM line says that we’re going to use an Amazon Linux box that has been pre-configured to have Python 3.8.
  • The RUN lines help us to install the Python libraries that we need.
  • The COPY lines say to move the corresponding files from the local directory to the current directory (./) of the Linux box.
  • The CMD line says that when you talk to the server, it should respond with the handler function from the app.py file.

Building a Docker image

Now, we’re going to build the Docker image. Make sure you’re in the proper directory and name the bot botfons-needles (or something else you’d like) by running the following command in the directory containing your Dockerfile:

docker build -t botfons-needles .

This command will probably take a while to run. It’s downloading everything to make a little Linux box that can run Python 3.8, and doing some additional tasks as specified by your Dockerfile. Once this process is done, set up a local server (on port 9000) for the bot where you can test it out by running

docker run -p 9000:8080 botfons-needles

In order to test your code, run the following cURL command in a new terminal:

curl -XPOST "http://localhost:9000/2015-03-31/functions/function/invocations" -d '{}'

If everything works, you’re ready to move on to the next step. More likely, something is a little off, and you’ll want to stop the server, and rebuild the image. To do this, find the name of the local server with

docker container ls

which will return a CONTAINER ID such as bb81431991sb. You can use this ID to stop the container, remove the container, and remove the image.

$ docker stop bb81431991sb
$ docker rm bb81431991sb
$ docker image rm botfons-needles

Then make your changes, and start again from the docker build step.

My first Twitter bot, @oeisTriangles. Read about this in my blog post, “Parity Bitmaps from the OEIS“.

Push to Amazon ECR

In this step, we’ll push our Docker image up to Amazon. So go to Amazon ECR, log in or create an account, navigate to the ECR console, and select “Create repository” in the upper right-hand corner.

This will bring up a place for you to create a repository name.

Now once there’s a repository available, we’re ready to push our local copy up. Start by downloading the AWS Command Line Interface and logging in to AWS. Notice that there are *two* references to the server location (us-east-1) and one reference to your account number (123456789).

aws ecr get-login-password --region us-east-1 \
| docker login --username AWS \
  --password-stdin 123456789.dkr.ecr.us-east-1.amazonaws.com

Now all you need to do is tag your docker image and push it up. Get your Docker image ID with docker image ls, and (supposing it’s 1111111111), tag it with the following command (again, making sure to specify the right server location and account number):

docker tag 1111111111 123456789.dkr.ecr.us-east-1.amazonaws.com/botfons-needles

Now you’re ready to push! Simply change 123456789 to your account number and us-east-1 to your server location in the following command and run it:

docker push 123456789.dkr.ecr.us-east-1.amazonaws.com/botfons-needles

Hook up to AWS Lambda

Now you’re ready to wire this up to a AWS Lambda function! Start by going to the AWS Lambda console and click “Create function”

To create a function from the AWS Lambda console, click “Create function” in the upper right-hand corner.

This will take you to a page where you’ll want to select the third option “Container image” at the top, give your function a name (e.g. myTwitterBot) and select the Container image URI by clicking “Browse images” and selecting the Docker image you pushed up.

Search for the image you just pushed up, choose the latest tag, and “Select image”.

Then the dialog will go away, and you can click “Create function”, after which your function will start to build—although it may take a while!

Next, you’ll want to test your function to make sure that it’s able to post to Twitter properly!

With the default RAM and time limit, it’s likely to time out. If the only thing that you’re using AWS for is posting this Twitter bot, then it doesn’t hurt to go to the “Configuration” tab and increase the memory and timeout under “General configuration”. (I usually increase Memory to 1024 MB and Timeout to 15 seconds, which has always been more than enough for me.)

Increase memory to 1024 MB and timeout to 15 seconds.

Run it on a timer

If things are running smoothly, then all that’s left to do is to set up a trigger. Do this by selecting “Triggers” in the “Configuration” tab, clicking “Add Trigger”, selecting “EventBridge (CloudWatch Events)”, and making a new rule with schedule expression rate(12 hours).

That’s it! AWS Lambda should trigger your bot on the interval you specified!

There’s only one step left: send me a message or tag me on Twitter @PeterKagey so that I can follow your new bot!

A π-estimating Twitter bot: Part II

This is the second part of a three part series about making the Twitter bot @BotfonsNeedles. Click here for Part I.

In this part, I’ll explain how to use the Twitter API to

  • post the images to Twitter via the Python library Tweepy, and
  • keep track of all of the Tweets to get an increasingly accurate estimate of 𝜋.

In the next part, I’ll explain how to

  • Package all of this code up into a Docker container
  • Push the Docker image to Amazon Web Services (AWS)
  • Set up a function on AWS Lambda to run the code on a timer

Get access to the Twitter API

When I made my first Twitter bot, I followed the article “How to Make a Twitter Bot With Python“.

In order to have your Python code post to your Twitter feed, you’ll need to register for a Twitter developer account, which you can do by going to https://developer.twitter.com/ and clicking apply. You’ll need to link the account to a phone number and fill out a few minutes of forms. For all four of my bots, (@oeisTriangles, @xorTriangles, @RobotWalks, and this bot) I’ve been approved right away.

Keep in mind that you can only use your phone number on two Twitter accounts—so you’ll have to use a Google Voice number or something else if you want to make more than two bots.

Once you’re approved, go to the Developer Portal, click on the Projects & Apps Overview, and click on the blue “+ New Project” button. You will be given a series of short questions, but what you submit isn’t too important.

Create a new project by clicking “+ New Project” via the Projects & Apps overview

Getting the API Keys

Once you’ve filled out the form, you should be sent to a page with an API Key and API Secret Key. This is essentially the password to your account, so don’t share these values.

A screenshot showing a (fake) API Key and API Secret Key.

We’re going to take these values and place them in a new file called secrets.py, which will look like this:

API_KEY        = "3x4MP1e4P1kEy"
API_SECRET_KEY = "5Ecr3TK3Y3x4MP1e4P1kEytH150nEi510nG"

Getting the Access Token

Once we close the API Key dialog, we’ll need to update our app to allow us to both read and write. We can do this by clicking on the gear to access our projects “App settings”.

Click the gear to access the settings.

Once you’re in, you’ll want to edit the App permissions to “Read and Write”.

Click “Edit” to update from “Read Only” to “Read and Write”.

Then go to the “Keys and Tokens” page (which you can do by clicking the key icon from the app settings page), and generate an Access Token and Secret.

Click “Generate” to make an access token and secret.

When you click “Generate” you should get an Access Token and a Access Token Secret, which you need to add to your secrets.py file.

Thus your secrets.py file should contain four lines:

API_KEY             = "3x4MP1e4P1kEy"
API_SECRET_KEY      = "5Ecr3TK3Y3x4MP1e4P1kEytH150nEi510nG"
ACCESS_TOKEN        = "202104251234567890-exTrAacC3551Nf0"


Next, we’ll hook this up to the Twitter API via tweepy, which I’ll install in the terminal using pip:

$ pip3 install tweepy

And make a file called twitter_accessor.py that looks exactly like this:

from secrets import *
import tweepy

class TwitterAccessor:
  def __init__(self):
    auth = tweepy.OAuthHandler(API_KEY, API_SECRET_KEY)
    auth.set_access_token(ACCESS_TOKEN, ACCESS_TOKEN_SECRET)
    self.api = tweepy.API(auth)

  def read_last_tweet(self):
    timeline = self.api.user_timeline(count=1, exclude_replies=True, tweet_mode='extended')
    return timeline[0].full_text

Next, we’ll check that everything is working by making a file called hello_twitter.py:

from twitter_accessor import TwitterAccessor

new_status = "Hello Twitter!"
print("Posted status: '" + new_status + "'")

Run it via the command line:

$ python3 hello_twitter.py

If something looks broken, try to fix it. (If it’s broken because of something I’ve said, let me know.)

Reading and writing Tweets

Now you can delete your hello_twitter.py file, because we’re about to do this for real! In part 3, we’re going to wire this up to AWS Lambda, which has certain preferences for how we structure things. With this in mind, I’d recommend following my naming conventions, unless you have a reason not to.

Each Tweet should have copy that looks like this:

This trial dropped 100 needles, 59 of which crossed a line. This estimates π ≈ 2*(100/59) ≈ 3.38, an error of 7.90%.

In total, 374 of 600 needles have crossed a line.
This estimates π ≈ 3.20, an error of 2.13%.

BotfonsNeedles should parse the “374 of 600”, throw 100 more needles, and report on the updated estimate of \(\pi\).

An implementation

I’ve made a file called tweet_builder.py, with five functions:

  • pi_digits_difference takes an estimate of \(\pi\) and outputs an appropriate length string. For example, if the estimate is \(3.14192919\), then it will output "3.14192", which are all of the correct digits, plus the first two that are wrong. If the estimate is \(3.20523\), then it will output “3.20".
  • error_estimate takes an estimate of \(\pi\) and computes the right number of digits to show in its percent error. For example, if the estimate is \(3.20523\) (which is \(2.0256396\%\) too big) then it will output "2.02%".
  • get_running_estimate uses the API in TwitterAccessor to look up the last tweet—then throws some needles, and outputs both the total number of needles tossed and the total number of needles that cross a line.
  • tweet_copy takes the information from get_running_estimate, formats it with pi_digits_distance and error_estimate and writes the text for the tweet.
  • post_tweet uses the API in TwitterAccessor to send the tweet to Twitter, with an image to match.

Most of these implementations are just details which can be found on Github, but I want to highlight post_tweet, the function that is likely to be the most relevant to you.

def post_tweet(self):
  file_name = self.drawer.draw_image()
  copy = self.tweet_copy()
  self.accessor.api.update_with_media(filename=file_name, status=copy)
  return copy

What’s next

In Part III, we’ll get this running in a Docker container and have it run on AWS Lambda.

If you want to get a head start, make a file called app.py with a function called handler, which AWS Lambda will call. This function should return a string, which will get logged.

from tweet_builder import TweetBuilder

def handler(event, context):
  return TweetBuilder().post_tweet()

As usual, if you have any questions or ideas, there’s nothing I love more than collaborating. If you want help getting your bot off the ground, ask me about it on Twitter, @PeterKagey!

A π-estimating Twitter bot: Part I

This is the first part of a three part series about making the Twitter bot @BotfonsNeedles. In this part, I will write a Python 3 program that

In the second part, I’ll explain how to use the Twitter API to

  • post the images to Twitter via the Python library Tweepy, and
  • keep track of all of the Tweets to get an increasingly accurate estimate of \(\pi\).

In the third part, I’ll explain how to

  • Package all of this code up into a Docker container
  • Push the Docker image to Amazon Web Services (AWS)
  • Set up a function on AWS Lambda to run the code on a timer
An example of dropping 100 needles in Buffon’s needle problem. Needles that cross a vertical line are colored red.

Buffon’s needle problem

Buffon’s needle problem is a surprising way of computing \(\pi\). It says that if you throw \(n\) needles of length \(\ell\) randomly onto a floor that has parallel lines that are a distance of \(\ell\) apart, then the expected number of needles that cross a line is \(\frac{2n}\pi\). Therefore one way to approximate (\pi) is to divide \(2n\) by the number of needles that cross a line.

I had my computer simulate 400 needle tosses, and 258 of them crossed a line. Thus this experiment approximates \(\pi \approx 2\!\left(\frac{400}{258}\right) \approx 3.101\), about a 1.3% error from the true value.

63/100 needles cross a vertical line; approximates \(\pi \approx 200/63 \approx 3.174\).
61/100 needles cross a vertical line; approximates \(\pi \approx 200/61 \approx 3.279\).
68/100 needles cross a vertical line; approximates \(\pi \approx 200/68 \approx 2.941\).
66/100 needles cross a vertical line; approximates \(\pi \approx 200/66 \approx 3.030\).

Modeling in Python

Our goal is to write a Python program that can simulate tossing needles on the floor both numerically (e.g. “258 of 400 needles crossed a line”) and graphically (i.e. creates the PNG images like in the above example).

The RandomNeedle class.

We’ll start by defining a RandomNeedle class which takes

  • a canvas_width, \(w\);
  • a canvas_height, \(h\);
  • and a line_spacing, \(\ell\).

It then initializes by choosing a random angle (\theta \in [0,\pi]) and random placement for the center of the needle in \[(x,y) \in \left[\frac{\ell}{2}, w -\,\frac{\ell}{2}\right] \times \left[\frac{\ell}{2}, h -\,\frac{\ell}{2}\right]\] in order to avoid issues with boundary conditions.

Next, it uses the angle and some plane geometry to compute the endpoints of the needle: \[\begin{bmatrix}x\\y\end{bmatrix} \pm \frac{\ell}{2}\begin{bmatrix}\cos(\theta)\\ \sin(\theta)\end{bmatrix}.\]

The class’s first method is crosses_line, which checks to see that the \(x\)-values at either end of the needle are in different “sections”. Since we know that the parallel lines occur at all multiples of \(\ell\), we can just check that \[\left\lfloor\frac{x_\text{start}}{\ell}\right\rfloor \neq \left\lfloor\frac{x_\text{end}}{\ell}\right\rfloor.\]

The class’s second method is draw which takes a drawing_context via Pillow and simply draws a line.

import math
import random

class RandomNeedle:
  def __init__(self, canvas_width, canvas_height, line_spacing):
    theta = random.random()*math.pi
    half_needle = line_spacing//2
    self.x = random.randint(half_needle, canvas_width-half_needle)
    self.y = random.randint(half_needle, canvas_height-half_needle)
    self.del_x = half_needle * math.cos(theta)
    self.del_y = half_needle * math.sin(theta)
    self.spacing = line_spacing

  def crosses_line(self):
    initial_sector = (self.x - self.del_x)//self.spacing
    terminal_sector = (self.x + self.del_x)//self.spacing
    return abs(initial_sector - terminal_sector) == 1

  def draw(self, drawing_context):
    color = "red" if self.crosses_line() else "grey"
    initial_point  = (self.x-self.del_x, self.y-self.del_y)
    terminal_point = (self.x+self.del_x, self.y+self.del_y)
    drawing_context.line([initial_point, terminal_point], color, 10)

By generating \(100\,000\) instances of the RandomNeedle class, and keeping a running estimation of (\pi) based on what percentage of the needles cross the line, you get a plot like the following:

This estimates \(\pi\approx 2\left(\frac{10000}{63681}\right) \approx 3.1407\) an error of 0.03%.

The NeedleDrawer class

The NeedleDrawer class is all about running these simulations and drawing pictures of them. In order to draw the images, we use the Python library Pillow which I installed by running

pip3 install Pillow

When an instance of the NeedleDrawer class is initialized, makes a “floor” and “tosses” 100 needles (by creating 100 instances of the RandomNeedle class).

The main function in this class is draw_image, which makes a \(4096 \times 2048\) pixel canvas, draws the vertical lines, then draws each of the RandomNeedle instances.

(It saves the files to the /tmp directory in root because that’s the only place we can write file to our Docker instance on AWS Lambda, which will be a step in part 2 of this series.)

from PIL import Image, ImageDraw
from random_needle import RandomNeedle

class NeedleDrawer:
  def __init__(self):
    self.width   = 4096
    self.height  = 2048
    self.spacing = 256
    self.random_needles = self.toss_needles(100)

  def draw_vertical_lines(self):
    for x in range(self.spacing, self.width, self.spacing):
      self.drawing_context.line([(x,0),(x,self.height)],width=10, fill="black")

  def toss_needles(self, count):
    return [RandomNeedle(self.width, self.height, self.spacing) for _ in range(count)]
  def draw_needles(self):
    for needle in self.random_needles:

  def count_needles(self):
    cross_count = sum(1 for n in self.random_needles if n.crosses_line())
    return (cross_count, len(self.random_needles))

  def draw_image(self):
    img = Image.new("RGB", (self.width, self.height), (255,255,255))
    self.drawing_context = ImageDraw.Draw(img)
    del self.drawing_context
    return self.count_needles()

Next Steps

In the next part of this series, we’re going to add a new class that uses the Twitter API to post needle-drop experiments to Twitter. In the final part of the series, we’ll wire this up to AWS Lambda to post to Twitter on a timer.