# Sequential addition of points on a circle, optimizing asymptotic packing radius

Suppose I have to put $N$ points $x_1, x_2, \ldots, x_N$ on the circle $S^1$ of length 1 so as to achieve the largest minimum separation (packing radius). The optimal solution is the equally spaced arrangement $x_k = k/N \bmod 1$, which achieves $\min_{i<j}\|x_i-x_j\| = 1/N$.

However, suppose I didn't know in advance the number of points $N$ I needed to put down, and I had to put the points down one by one. I am not allowed to move points that are already placed, and I want to a process that, when compared to the optimum for fixed $N$, does not do much worse asymptotically. More precisely:

Problem: what is the sequence $x_i$, $i=1,2,\ldots$, that maximizes $\mu=\lim\inf_N \min_{1\le i<j\le N} \|x_i-x_j\|/(1/N)$?

First example: greedy process. Put the first point at 0, the second at 1/2, and each subsequent one so as to bisect the largest empty interval. This gives $\mu=1/2$: when we add the $(2^k+1)$th point, the smallest separation is $1/2^{k+1}$.

Second example: Let $x_k = k\phi \bmod 1$, where $\phi = (\sqrt{5}-1)/2\approx 0.618$ is the inverse golden ratio. This process achieves $\mu = \phi$.

Conjecture: $\mu = \phi$ is optimal.

I think I know how to prove that among linear sequences $x_k = k\alpha$, that $\alpha=\phi$ (or any quadratic irrational whose continued fraction expansion is eventually all ones) is optimal, and I can start to imagine how to treat polynomial sequences, but I would like to know whether it's true that no sequence can do better than $x_k = k \phi$.

Motivation: the motivation for this problem is from phyllotaxis (the arrangement of leaves and other plant parts along the growth a plant), where the golden angle and Fibonacci spirals are commonly observed. Many researchers from D'arcy Thompson to Adil Mughal and Denis Weaire have proposed mechanism associated with packing problems that result in the appearance of the golden angle and Fibonacci spirals. From a mathematical perspective, the optimization problem I propose here seems to be the neatest and most succinct that has this pattern as a solution.

UPDATE: Christian Remling's answer, which I have accepted, provides convincing evidence that my conjecture is in fact wrong, which is very exciting. His answer gives an alternative process, which might be able to give the real optimal $\mu$. I added an answer heuristically analyzing Christian's process and arguing that $\mu=1/\log(4)\approx 0.721$ might be the optimal asymptotic separation that can be achieved with this process.

• This reminds me of the robot and apple problem, where a robot starts its search on a numberline for an apple placed somewhere on it. (Look for "cow path" in the literature where the search space looks like a branched continuum.) The line version seeks an optimal sequence of turning points so that the ratio search time (number of unit steps taken) to distance of apple from the start is minimized over all choices of distance (min of Max over n, you can get the technical definition from the literature). Gerhard "Wandering In Search Of Optimality" Paseman, 2017.07.10. – Gerhard Paseman Jul 10 '17 at 22:38
• I am reminded of ideas of discrepancy that show up in Uniform distribution theory (see the book of Kuipers and Niederreiter). Of course there are many ways to measure how uniformly distributed a sequence is. K and N look at (I believe) the maximal difference between the proportion of the first T terms of a sequence lying in I and the length of I; and then maximize over I. It's not clear how (if at all) this relates to your problem, but I would not be surprised if they have some lemmas that are useful for you. – Anthony Quas Jul 10 '17 at 23:11
• I have been calling this gradual-insertion kind of problem "peppering." I also arrived at the N * phi (mod 1) method from sunflowers! The phi method is still better than the methods below, for sampling a continuous function over an interval or sampling from unequal partitions of an interval (e.g. for timesharing). – SteveWithamDuplicate Sep 16 at 15:53

I will show $$\mu = \tfrac{1}{\log 4}$$.

I first prove the upper bound $$\mu \leq \tfrac{1}{\log 4}$$. Fix a positive integer $$r$$. For $$0 \leq k \leq r-1$$, let $$N_k$$ be $$2^{k/r} N$$ rounded to the nearest integer, where $$N$$ is a large parameter. Considering the $$N_{r-1}$$ intervals present at time $$N_{r-1}$$, let $$X_k$$ be the number of them created between times $$N_{k-1}$$ and $$N_k$$; with $$X_0$$ the number created before time $$N_0$$. Then $$X_0 \frac{\mu}{N_0} + X_1 \frac{\mu}{N_1} + \cdots + X_{r-1} \frac{\mu}{N_{r-1}} \leq 1.$$ Set $$S_k=X_0+X_1+\cdots +X_k$$; we can rewrite the above equation as $$\mu \left( S_0 \left(\frac{1}{N_0} - \frac{1}{N_1}\right) + S_1 \left(\frac{1}{N_1} - \frac{1}{N_2} \right) + \cdots + S_{r-2} \left(\frac{1}{N_{r-2}} - \frac{1}{N_{r-1}}\right) + S_{r-1} \frac{1}{N_{r-1}} \right) \leq 1.$$

Now, $$S_k$$ is the number of those intervals present at final time $$N_{r-1}$$ which were created before time $$N_k$$, so we have $$S_k \geq N_k - (N_{r-1} - N_k)=2 N_k - N_{r-1}$$. The coefficient $$\tfrac{1}{N_k} - \tfrac{1}{N_{k+1}}$$ on the left hand side is positive, so we have $$\mu \left( (2N_0-N_{r-1}) \left(\frac{1}{N_0} - \frac{1}{N_1}\right)+ (2N_1-N_{r-1}) \left(\frac{1}{N_1} - \frac{1}{N_2}\right)+\cdots+ (2N_{r-2}-N_{r-1}) \left(\frac{1}{N_{r-2}} - \frac{1}{N_{r-1}}\right)+\frac{N_{r-1}}{N_{r-1}} \right) \leq 1.$$

We gather terms with the same denominator to obtain $$\mu \left( \frac{2N_0-N_{r-1}}{N_0} + \frac{2 N_1 - 2 N_0}{N_1} + \cdots + \frac{2 N_{r-1} - 2 N_{r-2}}{N_{r-1}} \right) \leq 1$$ or $$\mu \left( 2r - \frac{N_{r-1}}{N_0} - \frac{2 N_0}{N_1} - \cdots - \frac{2 N_{r-2}}{N_{r-1}} \right) \leq 1.$$ Plugging in our optimized values: $$\mu \left( 2r - r 2^{(r-1)/r} \right) \leq 1$$ so $$\mu \leq \frac{1}{2r (1-2^{-1/r})}.$$

Finally, taking the limit as $$r \to \infty$$ gives $$\mu \leq \tfrac{1}{\log 4}$$ as promised.

Now, for the lower bound. Again, fix a positive integer $$r$$. Divide the circle into $$r$$ equal arcs; each arc will be subdivided into intervals. At any point in our process there will be some integer $$k$$ such that one arc is divided into some intervals of length $$2^{-k/r}$$ and some of length $$2^{-(k+r)/r}$$; the other $$r-1$$ arcs will contain intervals of only one length, namely $$2^{-j/r}$$ for $$k. (Obviously, I am ignoring rounding issues.) As time ticks on, we subdivide the intervals in the first arc until they are all split in half, then move to the next arc, and so on.

This achieves (again, ignoring rounding) $$\mu = 2^{-(k+r)/r} \left( \frac{2^{k/r}}{r} +\frac{2^{(k+1)/r}}{r} + \cdots \frac{2^{(k+r-1)/r}}{r} \right) = \frac{1}{2r(2^{1/r}-1)}.$$ Again, sending $$r \to \infty$$ achieves $$\tfrac{1}{\log 4}$$.

Here is a slicker, though less motivated way, to achieve the lower bound. We will number the points starting with $$x_0$$ at position $$0$$. For $$n >0$$, write $$n = 2^q+r$$ with $$0 \leq r <2^q$$ and put $$x_n$$ at $$\log_2\left(\tfrac{2n+1}{2^{q+1}}\right)$$. So the first few points are $$0$$, $$\log_2(1+1/2)$$, $$\log_2(1+1/4)$$, $$\log_2(1+3/4)$$, $$\log_2(1+1/8)$$, $$\log_2(1+3/8)$$, $$\log_2(1+5/8)$$, $$\log_2(1+7/8)$$, ... When point $$x_n$$ is inserted, the smallest interval will either be the one just to the right of $$x_n$$ or else the furthest right interval. Approximating $$\log$$ by its linearization, these have lengths roughly $$\tfrac{2^{q+1}}{(2n+1) \log 2} \tfrac{1}{2^{q+1}}=\tfrac{1}{(2n+1) \log 2}$$ and $$\tfrac{1}{2 \log 2} \tfrac{1}{2^q}=\tfrac{1}{2^{q+1} \log 2}$$ respectively.

Multiplying these by $$n+1$$ (the number of intervals) gives $$\tfrac{n+1}{(2n+1) \log 2}$$ and $$\tfrac{n+1}{2^{q+1} \log 2}$$. The former approaches $$\tfrac{1}{\log 4}$$ while the latter oscillates between $$\tfrac{1}{\log 2}$$ and $$\tfrac{1}{\log 4}$$.

• Thank you for this excellent analysis! I do not doubt that the issues of rounding and the issue of the inequality (*) only applying asymptotically (i.e. in terms of liminf, and not in terms of an inequality for large enough N) are easy to deal with. – Yoav Kallus Aug 9 '17 at 6:24
• A slightly simpler version: Start n at 1. Point n is at log_2(2 n -1), and the smallest interval is log_2(2 n ) - log_2(2 n -1). @David's method has better numerical precision. Two nice features of this method: 1) The ratio between the largest and smallest interval is always less than 2. 2) It resembles how a Huffman code is built (always combining the smallest two), but in reverse. – SteveWithamDuplicate Sep 16 at 15:32

This is really a comment, but it's getting a bit long for the comment box.

I want to point out that in addition to what you call the "greedy process," there's another obvious attempt, which could be called the "stingy process."

Let a target $\mu$ be given. There's basically one rule, namely don't squander unused space unless absolutely forced to. More formally, if I have already chosen points $0\le x_1, \ldots , x_N\le a$, then these define $N$ intervals $I_j$. I must now place the next point such that the minimum distance is $\ge\mu/(N+1)$ (I will keep the whole sequence $\ge\mu$, not just the $\liminf$). If I have an interval of length at least twice this distance, then I put my point inside such an interval (let's say inside the smallest one that works and I will also make one distance equal to $\mu/(N+1)$, though that might not be optimal). If not, then I grudgingly add an interval of that length to the right end of the current configuration.

This gives me intervals that get subdivided, and every once in a while a new interval gets added at the current right endpoint. The whole procedure will be a success if the limit of these right endpoints is $\le 1$.

I fooled around some with this; let's take $\mu=1$. Then the first few points are $$0, 1/2, 1/2+1/3, 1/4, 1/2+1/3+1/5, 1/2+1/6, \ldots$$ It seems that $\lim a_N =1/2+1/3+1/5+1/7+1/11+1/14+\ldots \simeq 1.4$. If this is correct (but I'm not making any strong claims it is), then it means that we can get $\mu\simeq 1/1.4\simeq 0.7$, which is better than the golden mean.

In any event, it should be easy (for someone more computer savvy) to do the numerics and get a more reliable estimate.

• Remarkable. I (rather inefficiently, to start) coded up your stingy process in C++ (github.com/ykallus/stingy-point-sequence). The process terminates after only 26 points for $\mu=0.72$. For $\mu=0.715$, the process was able to place $> 40000$ points before I killed it. So, clearly, my original conjecture was wrong. – Yoav Kallus Jul 11 '17 at 15:35
• Here is a plot of the first 10000 points for $\mu = 0.71$: i.imgur.com/Wm1jx9h.png – Yoav Kallus Jul 11 '17 at 15:45
• Though the stingy process from a cold start terminates for $\mu=0.72$, using various schemes for warm starting, I'm able to get the process to continue seemingly indefinitely for $\mu=0.72$ and even slightly above. – Yoav Kallus Jul 11 '17 at 20:13
• It is not clear why the last interval (between say $1/2+1/3+1/5+1/7+1/11+1/14$ and $1/\mu$) satisfy the minimum distance condition. – Alexey Ustinov Jul 12 '17 at 4:22
• @AlexeyUstinov: Well, I was just throwing out an idea in this post, not trying to get all small details addressed. This, however, is clearly not a problem: we can increase the interval by $\epsilon$, and then the interval between $a_n$ and $1+\epsilon$ will always be long enough for large $n$ as long as $a_n\le 1$. – Christian Remling Jul 12 '17 at 16:34

In this second answer, I want to discuss an upper bound on $\mu$ (not optimal, and I don't think this argument could give an optimal bound, even after fine tuning).

Suppose we have placed $N$ points such that (as required) the corresponding $N$ intervals $I_j$ have lengths $|I_j|\ge \mu/N$. Let's now position the next $N$ points. More specifically, let's denote by $k_j\ge 0$ the number of new points that I put into $I_j$. For this to be possible, I obviously need $|I_j|\ge (k_j+1)\mu/(2N)$.

Let $1\le S\le N$ be the number $j$'s with $k_j\ge 1$ (that is, the number of intervals that I actually use). Then $$(N-S)\frac{\mu}{N} + \sum (k_j+1)\frac{\mu}{2N} \le 1 ,$$ and since $\sum k_j=N$, this gives that $$\mu\le \frac{2}{3-S/N} .\quad\quad\quad\quad (1)$$ Of course, this is not useful if $S/N\simeq 1$.

So we now need to look at this case $S\ge (1-a)N$ separately. I want to bound from below the lengths of these $S$ intervals. The worst case scenario arises when we first position the $aN$ remaining points, and then put $S$ points into that many new intervals. Then the $j$th such interval must have length $\ge 2\mu/(N+aN+j)$, for a total length of at least $$2\mu \sum_{j=1}^{(1-a)N} \frac{1}{(1+a)N+j}\simeq 2\mu \int_0^{(1-a)N} \frac{dx}{(1+a)N+x} =2\mu\log \frac{2}{1+a} .$$ Thus $\mu\le \frac{1}{a+2\log 2/(1+a)}$ (the extra $a$ comes from the length $\ge a\mu$ of the remaining $Na$ intervals), and when combined with (1), this gives the unconditional bound $$\mu\le\max \left\{ \frac{2}{2+a} , \frac{1}{a+2\log 2/(1+a)} \right\} ,$$ valid for all $0<a\le 1$.

The optimal $a$ seems to be around $a=0.3$, and then this gives $\mu\le 0.87$.

Thinking more about Christian's stingy process, I have a new conjecture for the optimal $\mu$. I motivate the conjecture by modeling the evolution of the distribution of empty interval lengths in the continuum limit. Let $t=N$ be the number of steps the process has taken, and let $t p(x,t)dx$ be the number of intervals with length between $xt$ and $(x+dx)t$. There are two normalization conditions on $p(x,t)$:

1. There are a total of $t$ intervals: $\int_0^{\infty} p(x)dx = 1$.
2. The interval lengths add up to 1: $\int_0^{\infty} xp(x)dx = 1$.

There are three processes that change $p(x,t)$ over time:

1. $x=\text{length}\times t$ increases.
2. $p = \text{count}/t$ decreases.
3. One interval of length around $2\mu/t$ (i.e. $x=2\mu$), gets removed, and replaced by two intervals of length around $\mu/t$ each.

As a PDE, this is:

$$\frac{\partial p}{\partial t} = -\frac{2p}{t} - \frac{x}{t}\frac{\partial p}{\partial x} + \frac{2}{t}\delta(x-\mu) -\frac{1}{t}\delta(x-2\mu)$$

If the stingy process continues forever, we expect it to reach a stationary distribution $p(x)$, which would be governed by the ODE:

$$\frac{dp}{dx} = \frac{1}{x}\left[2\delta(x-\mu)-\delta(x-2\mu)-2p\right]$$

The solution to the ODE away from the source terms is $p(x)=c/x^2$, so we get $p(x) = c_0 /x^2$ for $x<\mu$, $p(x) = c_1 /x^2$ for $\mu<x<2\mu$, and $p(x) = c_2 /x^2$ for $x>\mu$. For normalization not to blow up, $c_0=0$. From matching the discontinuities to the source terms, $c_1 = 2\mu$ and $c_2=0$. Finally, from the second normalization condition we get $\mu=1/\log(4) \approx 0.72135$.

This matches pretty well with the limit where I was able to get the stingy process to not terminate in my numeric experiments. So, the new conjecture is that the optimal $\mu$ is $1/\log(4)$.