Analysis of Algorithms for Moser's Problems on Sums of Consecutive Primes

TL;DR

This paper analyzes two algorithms for computing the number of representations of integers as sums of consecutive primes, providing complexity bounds and generating data to address four longstanding open questions posed by Moser.

Contribution

It offers a detailed complexity analysis of two existing algorithms and applies them to generate data that resolves four open problems in number theory.

Findings

First algorithm runs in O(x log x) time with x^{2/3} space.

Second algorithm has two versions: one with O(x log x) time and x^{3/5} space, another with O(x(log x)^2) time and O(√(x log x)) space.

Data produced addresses all four open questions posed by Moser.

Abstract

In his 1963 paper on the sum of consecutive primes, Moser posed four open questions related to $f(n)$, the number of ways an integer $n$ can be written as a sum of consecutive primes. (See also problem C2 from Richard K.~Guy's \textit{Unsolved Problems in Number Theory}.) In this paper, we present and analyze two algorithms that, when given a bound $x$, construct a histogram of values of $f(n)$ for all $n\le x$. These two algorithms were described, but not analyzed, by Jean Charles Meyrignac (2000) and Michael S. Branicky (2022). We show the first algorithm takes $O(x\log x)$ time using $x^{2/3}$ space, and the second has two versions, one of which takes $O(x\log x)$ time but only $x^{3/5}$ space, and the other which takes $O(x(\log x)^2)$ time but only $O( \sqrt{x\log x})$ space. However, Meyrinac's algorithm is easier to parallelize. We then present data generated by these algorithms that address all four open questions.

Figures

Tables

Table 1. Table 1: Frequency data up to 10 14 10^{14}

Count	${10}^2$	${10}^3$	${10}^4$	${10}^5$	${10}^6$	${10}^7$	${10}^8$	${10}^9$
0	46	520	5191	51462	518001	5205110	52209546	522955756
1	38	310	3290	34538	344100	3427038	34146573	340693986
2	14	140	1213	11236	111132	1099545	10950371	109272550
3	2	28	275	2396	22916	228659	2292360	23003362
4	0	0	29	323	3409	35009	353614	3584873
5	0	2	2	44	403	4197	42946	440748
6	0	0	0	1	37	412	4205	44623
7	0	0	0	0	2	28	356	3793
8	0	0	0	0	0	2	28	299
9	0	0	0	0	0	0	1	10

Table 2. Table 3: h ​ ( k ) h(k) is heuristically Poisson with mean λ = log ⁡ 2 \lambda=\log 2

$k$	$h\left(k\right)$	$n\mathrm{Pr}\left[X=k\right]$
0	52255406573294	50000000080183
1	33983734548972	34657359083576
2	10980796355393	12011325367217
3	2351331657326	2775205437692
4	375496312243	480906456153
5	47717060499	66667790839
6	5027735200	7701765209
7	451927961	762636691
8	35376934	66077434
9	2452073	5089043
10	151480	352746
11	8546	22228
12	430	1284
13	14	68
14	1	3

Table 3. Table 4: Smallest n n with f ​ ( n ) = k f(n)=k representations

$k$	$n$
0	1
1	2
2	5
3	41
4	1 151
5	311
6	34 421
7	218 918
8	3 634 531
9	48 205 429
10	1 798 467 197
11	12 941 709 050
12	166 400 805 323
13	6 123 584 726 269
14	84 941 668 414 584

Table 4. Table 5: The 14 representations for g = 84 941 668 414 584 g=84\,941\,668\,414\,584 .

$\ell$	$p_{\min }$	$p_{\max }$
2 117 074	21 797 833	58 785 359
361 092	231 753 581	238 710 779
288 268	291 853 531	297 473 801
199 390	424 030 259	427 989 799
112 544	753 590 641	755 886 067
73 026	1 162 407 049	1 163 930 791
68 854	1 232 927 929	1 234 369 457
296	286 965 092 209	286 965 099 727
294	288 917 235 553	288 917 243 497
206	412 338 193 609	412 338 198 731
146	581 792 247 697	581 792 251 207
86	987 693 817 667	987 693 819 859
26	3 266 987 246 389	3 266 987 247 019
2	42 470 834 207 273	42 470 834 207 311

Equations

$$2357=773+787+797=461+463+467+487.$$

$$\frac{1}{x}\sum _{n=1}^{x}f(n)=\log 2+o(1)$$

$$h(k)=\#\{n\le x:f(n)=k\}.$$

$$2+3+5+\cdots +p_M\le x<2+3+5+\cdots +p_M+p_{M+1}.$$

$$\frac{x}{m+2}<p(m,x)\le \frac{2x}{m}(1+o(1)).$$

$$p_{\ell +1}+p_{\ell +2}+\cdots +p_{\ell +m}\le x<p_{\ell +1}+p_{\ell +2}+\cdots +p_{\ell +m}+p_{\ell +m+1},$$

$$x$$

$$\sum _{i=1}^u{p}_i\sim \sum _{i=1}^{u}i\log i\sim \frac{u^2\log u}{2}$$

$$x$$

$$M(x)\sim \frac{2\sqrt{x}}{\sqrt{\log x}}.$$

$$\frac{M^2\log M}{2}\le x(1+o(1))<\frac{(M+1{)}^2\log (M+1)}{2}.$$

$$n(\mathbf{slide}(c))-n(c)=p-q=\Theta (m\log (n/m))$$

$$\Delta /(m\log (x/m)),$$

$$m_{cutoff}=\Theta (x_1/\Delta ).$$

$$n(c_b)\ge x_1$$

$$n(c_a)\ge x_1$$

$$|n(c_a)-n(c_b)|=O(m\log (x/m))$$

$$n(c_a)=n({\mathbf{slide}}^s(c_b))-p(m+1,x_1).$$

$$n(c_a)=n(c_b)+\Theta (sm\log (x_1/m))-p(m+1,x_1)$$

$$p(m+1,x_1)=n(c_b)-n(c_a)+\Theta (sm\log (x_1/m))$$

$$p(m+1,x_1)=\Theta (sm\log (x_1/m))$$

$$\frac{x_1}{\log x_1}\approx s\cdot m_{cutoff}2$$

$$s=\Theta \left(\frac{\Delta }{m_{cutoff}\log x}\right),$$

$$\frac{x}{\Delta }\left(\Delta +\sum _{m=m_{cutoff}}^M\frac{\Delta }{m\log x}\right).$$

$$\frac{\Delta }{m}+m\log x$$

$$m$$

$$\frac{x}{\Delta }\sum _{m=1}^{m_{cutoff}}\left(\frac{\Delta }{m}+m\log x\right)$$

$$\Delta \gg m_{cutoff}^2.$$

$$\frac{x}{\Delta }\sum _{m=1}^{m_{cutoff}}\left(\frac{\Delta }{m}+m\log x+m+\frac{\Delta }{m\cdot \log x}\right)$$

$$x\log m_{cutoff}+m_{cutoff}^3\cdot \log x.$$

$$O\left(x\log x+\frac{(x^{1/3+u}{)}^{2}x\log x}{x^{2/3-u}}\right)=O(x^{1+3u}\log x).$$

Full text

1 Introduction

We say a positive integer $n$ is gleeful if it can be written as the sum of consecutive primes. So, for example, $n=10$ is gleeful since $10=2+3+5$. Some integers have multiple gleeful representations, or distinct ways of expressing them as such a sum; for instance, $17$ has two representations: $2+3+5+7$ and $17$ (where we consider the empty sum of a single prime to be gleeful), and for a somewhat larger example, we have

[TABLE]

Let $f(n)$ count the number of gleeful representations for the integer $n$. Thus, $f(6)=0$, $f(10)=1$, and $f(17)=f(2357)=2$.

In his 1963 paper, Moser [7] proved that the average value of $f(n)$ is asymptotically equal to $\log 2$; that is,

[TABLE]

for large $x$. He then posed the following four questions:

Is $f(n)=1$ infinitely often? 2. 2.

Is $f(n)=k$ solvable for every integer $k\geq 0$? 3. 3.

Does the set of numbers $n$ such that $f(n)=k$ have positive density for every integer $k\geq 0$? 4. 4.

Is $\limsup_{n\rightarrow\infty}f(n)=\infty$?

These questions are then reprised in problem C2 from [5], and several Online Encyclopedia of Integer Sequences (OEIS) entries address these questions: [1]:

A054859, A054845, A067381, A054996, A054997, A054998, A054999, A055000, A055001 .

To address Moser’s four questions, in this paper, we describe and analyze variations of two known algorithms that generate data on this function $f$. As we will demonstrate, that data suggests the answers are all ”yes”. In particular, given a bound $x$, our algorithms compute $f(n)$ for all $n\leq x$ and then construct a histogram $h(k,x)$ (oftentimes shorthanded to $h(k)$ where $x$ is implicit) where

[TABLE]

A first attempt at an algorithm to generate the values $h(k)$ might be like Algorithm 1.

This algorithm takes time linear in $x$, which is optimal for this problem. But it requires space linear in $x$ to hold the array $f$ (plus the prefix sums array - see [8]). This is problematic on currently available hardware if we want to compute the histogram $h$ for, say, $x=10^{14}$. Thus, the algorithm variants we present here focus on reducing memory use while keeping the running time as close to $O(x)$ as possible.

Throughout this paper, we count arithmetic operations on integers with $O(\log x)$ bits, including branching and array indexing, as taking constant time, and we count space used in bits. In practice, all of our variables are 64-bit integers or smaller.

We now state our main results and outline the rest of this paper.

In §2, we look at an algorithm due to Jean Charles Meyrignac (2000) that appeared in a puzzle solution to Problem #46, posted on the primepuzzle.net website curated by Carlos Rivera [10]. Henceforth we call this the puzzle algorithm. We analyze a variant of the puzzle algorithm, mostly by filling in the details of Meyrinac’s description. The algorithm parallelizes nicely, has an interesting time/space tradeoff, and generates values of $f(n)$ for all $n$ from an interval.

In §3, we look at an algorithm due to Michael S. Branicky (2022), which was written in Python and appears on the OEIS [1, A054845]. Again, we analyze two variants of this priority-queue-based algorithm, depending on how primes are generated. These algorithms have the property that $f(n)$ values are generated one at a time in increasing order based on $n$. Although these algorithms present a better time/space tradeoff, they resist parallelization. Nevertheless, we were able to run one of our priority-queue variants up to $x=10^{12}$ using a single core, and up to $x=10^{13}$ using 156 cores.

In §4, we present histograms for $x$ up to $10^{14}$. We also discuss our implementation details. The curious reader may find our code and data here: https://eleanorwaiss.github.io/gleeful.

Finally, in §4.1 we give conjectured answers to Moser’s questions based on our data.

2 The Puzzle Algorithm

The puzzle algorithm, as we call it, is based on the following description, written by Jean Charles Meyrignac [10]:

1. First, pre-compute a table of primes
2. Lower bound = 0
3. Higher Bound=Lower Bound + size of sliding window
4. Compute all consecutive sums of n elements for every element of the pre-computed table and counts every occurrence.
5. Display the best occurrences
6. Change the lower bound and go to 3)

Our version of this algorithm varies in several ways, but the main idea is the same. Values of $f(n)$ are computed for every integer $n$ in a segment or window, and then these are scanned to add to a running histogram. We spend some time discussing the details of how step (4) above can be done efficiently. But before we dive into the details, we need to establish some notation and prove some useful lemmas.

2.1 Notation and Lemmas

Let us introduce some notation.

Let $x$ be the limit for our computation. The interval $[1,x]$ is divided into equal-length segments (or windows) of size $\Delta$ each (the sliding window). The last segment may be shorter if $\Delta$ does not divide $x$. There are $\lceil x/\Delta\rceil$ segments total. For a particular segment, let $x_{1}$ be the left endpoint and $x_{2}$ be the right endpoint, but just outside the window. We have $\Delta=x_{2}-x_{1}$, and we assume that $\Delta$ will be a fractional power of $x$, such as $\sqrt{x}$ or $x^{2/3}$.

The basic idea is that the segment represents the values of $f(n)$ for $x_{1}\leq n<x_{2}$, which takes $\Delta$ small integers to represent (an 8-bit integer suffices in practice).

We use the term chain to indicate a string of consecutive primes, and we use $m$ to denote the number of primes, or length, of the chain. Then $n$ is the sum of those primes. If we write $c$ as a variable representing a chain, then $m(c)$ is its length and $n(c)$ is its sum. In practice, we store the start and stop prime indices, the length, and the sum in a chain variable.

For example, if $c:=[2,3,5,7]$, then $m(c)=4$ and $n(c)=17$. We will often wish to slide a chain, which means dropping its smallest prime and adding the next larger prime. A slide operation preserves length but increases the sum. Thus, ${\mathbf{slide}}(c)=[3,5,7,11]$. If a list of primes has been precomputed, we can implement a ${\mathbf{slide}}$ operation so that it takes constant time.

Let $M:=M(x)$ be the length of the longest chain of consecutive primes whose sum is $\leq x$. WLOG we may assume $2$ is in the chain so that

[TABLE]

For any chain of length $m$ with sum $n\leq x$, we have $m\leq M$.

Let $p(m)=p(m,x)$ be the largest prime that participates in any chain with sum $n\leq x$ of length $m$.

Lemma 2.1**.**

We have

[TABLE]

Proof.

There exists an integer $\ell>0$ such that

[TABLE]

so that $p(m,x)=p_{\ell+m}$ here.

For the lower bound, we have

[TABLE]

by Bertrand’s postulate.

For the upper bound, we note that

[TABLE]

(see Bach and Shallit [3], chapter 2). We then have

[TABLE]

by the prime number theorem. ∎

We have the following estimate for $M=M(x)$:

Lemma 2.2**.**

[TABLE]

Proof.

Consider a chain of length $M$ that consists of the first $M$ primes. From (2), we have

[TABLE]

Taking logarithms, we deduce that $\log x\sim 2\log M$. The result then follows. ∎

And we have the following bound on sums of slided chains.

Lemma 2.3**.**

If $c$ is a chain of length $m$ and sum $n$, then $n({\mathbf{slide}}(c))-n(c)$ has average value $\Theta(m\log(n/m))$.

Proof.

A single ${\mathbf{slide}}$ operation changes the sum of the chain by $p-q$, where $p$ is the newest prime added to the chain and $q$ was the smallest prime just removed from the chain by sliding. By Lemma 2.1, we know $p=\Theta(n/m)$, and so by the prime number theorem, we know consecutive primes in the chain are, on average, $\log(n/m)(1+o(1))$ apart. Thus $p\approx q+m\log(n/m)$, or

[TABLE]

on average. ∎

Corollary 2.4**.**

The number of chain sums $n$ from chains of length $m$ that fall in a segment of length $\Delta$ is, on average, proportional to

[TABLE]

where $x$ is the midpoint of the interval.

2.2 Algorithm Description

We are now ready to describe our version of the puzzle algorithm. The input is a bound $x$, and the output is an array representing the histogram $h$, where $h(k)=\#\{n\leq x:f(n)=k\}$.

In step (4) of Meyrignac’s description, the basic idea is to work through the chains from maximum length $M$ down to length 1, generating all relevant chains of the current length $m$ before moving on. If the length $m$ is ”large”, we build the first chain of length $m$ from the previous first chain of length $m+1$ by dropping the largest prime and then sliding as needed. (If we were to build the chain of length $m$ from scratch, this takes $\sum m\approx M^{2}$ time per segment, which is way too slow.)

The original chain of length $M$ is constructed by finding primes until the sum is maximal but less than $x_{2}$.

When $m$ gets small enough, the first chain of length $m$ has no overlap with the previous one of length $m+1$, and may even be far enough away that making the first $m$-length chain from scratch is more efficient. This is how we define ”small” values of $m$. We define ${m_{\mathit{cutoff}}}$ as the largest ”small” value of $m$.

To set a more precise value for ${m_{\mathit{cutoff}}}$, we set the work threshold to the the time needed to process the segment for that value of $m$ after that first chain is found. By Corollary 2.4, this is $O(\Delta/({m_{\mathit{cutoff}}}\log x))$ operations.

A key in the analysis is estimating this cutoff point, which we do now.

Lemma 2.5**.**

On average, we have

[TABLE]

Proof.

Let $c_{a}$ be the first chain of length $m$ we use (after), and $c_{b}$ be the first chain of length $m+1$ (before). Both chains have sums near $x_{1}$, as in

[TABLE]

are both true. By Lemma 2.3, we can deduce that

[TABLE]

on average.

By Lemma 2.1, the largest prime in $c_{b}$ is $p(m+1,x_{1})=\Theta(x_{1}/m)$. Let $s$ denote the number of ${\mathbf{slide}}$ operations to obtain $c_{a}$ from $c_{b}$ after $p(m+1,x_{1})$ is removed from $c_{b}$. Then we have

[TABLE]

Note that removing $p(m+1,x_{1})$ before or after sliding gives the same sum.

Sliding a chain of length $m$ with sum near $x$ changes its sum by $\Theta(m\log(x/m))$ by Lemma 2.3. This gives us

[TABLE]

or

[TABLE]

on average. Using Equation (3) then gives

[TABLE]

or, with Lemma 2.1,

[TABLE]

on average for when $m={m_{\mathit{cutoff}}}$. But, we have

[TABLE]

as the number of slide operations matches the work for the segment for that value of $m$ by Corollary 2.4. A bit of simplification completes the proof. ∎

With an estimate for ${m_{\mathit{cutoff}}}$ in hand, we are now ready to describe the Puzzle Algorithm is detail (Algorithm 2).

2.3 Time and space analysis

In our analysis, we assume that $\Delta$ is a fractional power of $x$ and, in fact, at least $M(x)$. We assume that $x_{2}\approx x$. The number of loop iterations when $x_{2}$ is substantially smaller than $x$ is limited, and we can assume, therefore, that this case contributes negligibly.

The space used by our algorithm will be at least $\Delta$ bits. If possible, we want to keep space usage as low as possible while keeping our running time as close to linear in $x$ as we can.

Since the average value of $f(n)$ is constant, most values of $f(n)$ can be represented in a constant number of bits, with larger values stored in a secondary data structure to keep the overall space use linear. In practice, we simply used an array of 8-bit integers. We have not yet found an integer $n$ with $f(n)>14$.

Throughout, we mention the use of average-case analysis. For example, consecutive primes near $x$ differ by $\log x$ on average, by the prime number theorem. For a particular step on a particular segment, this will vary, but over the course of the entire algorithm, significant variations from the average must cancel out, as the algorithm examines all primes and all chains up to $x$.

2.3.1 Preprocessing

It only takes $O(\Delta)$ time to do steps 1 to 2. In particular, finding the primes up to $O(\Delta)$ can be done in linear time using any number of prime sieves, including the sieve of Eratosthenes with a wheel or the Atkin-Bernstein sieve [2]. In practice, we use an incremental sieve [12] to find more primes as needed, since we don’t know the implied constant in $O(\Delta)$.

2.3.2 Main Loop

We move on to the the main loop, Line 2, which repeats $O(x/\Delta)$ times.

Inside the main loop, all the steps outside of processing large and small values of $m$ can be easily seen to take at most $O(\Delta)$ time, especially when we recall that $M\ll\Delta$.

2.3.3 Large $m$

Processing chains with large $m$ begins on Line 2. Everything, aside from the loops starting on Lines 2 and 2, takes constant time. Summing over all large $m$ gives us $O(M)=O(\Delta)$.

By Corollary 2.4, the loop at Line 2 takes at most $O(\Delta/(m\log(x/m)))$ time. We can simplify this to $O(\Delta/(m\log x))$ because $1\leq m\leq M=\Theta(\sqrt{x/\log x}$ by Lemma 2.2, so that $\log x=\Theta(\log(x/m))$.

The value of $s$ in Line 2 is chosen specifically to make sure that the cost of the loop at Line 2 is not more than the cost of the one at Line 2.

For large $m$, then, we have a cost proportional to

[TABLE]

By Lemma 2.2, this simplifies to $O(x)$.

2.3.4 Small $m$

Processing chains with small $m$ begins on Line 2. For small $m$, the largest prime of any chain of length $m$ will range from $p(m,x_{1})$ to $p(m,x_{2})$. This explains the need for Line 2. By Lemma 2.1 this is from roughly $x_{1}/m$ to $x_{2}/m$, or an interval of size $O(\Delta/m)$. The smaller primes of any chain either fall in that interval, or just below it by roughly $m\log x$ by the prime number theorem.

Now, the Atkin-Bernstein sieve will find primes of size $y$ on an interval of size $y^{1/3}$ or larger in time linear in the interval length, using Galway’s improvements [4]. Here, $y=p(m,x_{2})=O(x/m)$, so we need

[TABLE]

Using Lemma 2.5 and a bit of work, this gives us an upper bound of $x^{2/5}$ for ${m_{\mathit{cutoff}}}$, as $m\leq{m_{\mathit{cutoff}}}$ for ”small” $m$. As we will see later, ${m_{\mathit{cutoff}}}\approx x^{1/3}$, and so the Atkin-Bernstein sieve with Galway’s improvements will work fine here. (Note that in practice we used something akin to the pseudosquares prime sieve [11], but with Miller-Rabin prime tests using particular bases to achieve deterministic primality testing.)

Summing this over all segments and over all $m\leq{m_{\mathit{cutoff}}}$ gives a time bound of

[TABLE]

To keep this close to linear, we want

[TABLE]

Line 2 can be done in $O(m)$ time now that the necessary primes are available.

Line 2 takes at most $O(\Delta/(m\log x))$ time by Corollary 2.4, since the primes are available so that ${\mathbf{slide}}$ operations take constant time.

Theorem 2.6**.**

On input $x>0$, the Puzzle Algorithm computes a histogram $h$ as described above in $O(x\log x)$ arithmetic operations and $O(x^{2/3})$ space.

Proof.

From our discussion above, everything outside processing small values of $m$ takes $O(x)$ time.

The equation to describe the asymptotic running time, up to a constant multiplicative factor, for small values of $m$, is

[TABLE]

This argues for making ${m_{\mathit{cutoff}}}$ as small as possible, but we have Lemma 2.5. Multiplying the numerator and denominator of the second term by ${m_{\mathit{cutoff}}}$ and writing $\Delta\cdot{m_{\mathit{cutoff}}}=\Theta(x)$ gives

[TABLE]

We obtain that ${m_{\mathit{cutoff}}}=x^{1/3}$ and $\Delta=x^{2/3}$ and the proof is complete.

∎

On current hardware, $x^{2/3}$ space is problematic if we want to compute the histogram $h$ up to, say, $x=10^{15}$, as this means using $10^{10}$ space, which is too much. What happens if we allow our runtime to creep up - does this allow for less space use? The answer is yes, as the runtime is inversely proportional to $\Delta$. If we set $\Delta=x^{2/3-u}$, then by Lemma 2.5 we have ${m_{\mathit{cutoff}}}=O(x/\Delta)=x^{1/3+u}$, then the runtime becomes

[TABLE]

So for example, setting $\Delta=\sqrt{x}$, or $u=1/6$, gives a running time of $O(x^{3/2}\log x)$.

3 The Priority Queue Algorithm

The main idea behind the priority queue algorithm is to generate all the chains with sums between $1$ and $x$, our upper bound, and insert them into a priority queue (a min-heap) using their sums as their priority. Thus, the front of the queue will have the chains of smallest sums. Then we can find all chains that sum to 1 from the front of the queue, counting them as we remove them from the priority queue, and incrementing the histogram as appropriate. Rinse and repeat with $2,3,4,\ldots,x$. This works, but there are too many chains to store all at once – $O(x)$ according to Moser. However, observe that a particular value of $n$ has at most one representation of any given length. This means we need only store one chain of each length in the priority queue. By Lemma 2.2 we know that this is $O(\sqrt{x/\log x})$, much better than $O(x)$. When a chain is removed from the priority queue, we simply perform a ${\mathbf{slide}}$ and then insert it back into the queue.

This is precisely the idea in Michael S. Branicky’s algorithm, described as Python code from [1, A054845] in 2022:

(Python) # alternate version for going to large n import heapq from sympy import prime from itertools import islice def agen(): # generator of terms p = v = prime(1); h = [(p, 1, 1)]; nextcount = 2; oldv = ways = 0 while True: (v, s, l) = heapq.heappop(h) if v == oldv: ways += 1 else: yield ways for n in range(oldv+1, v): yield 0 ways = 1 if v >= p: p += prime(nextcount) heapq.heappush(h, (p, 1, nextcount)) nextcount += 1 oldv = v v -= prime(s); s += 1; l += 1; v += prime(l) heapq.heappush(h, (v, s, l)) print(list(islice(agen(), 102))) # Michael S. Branicky, Feb 17 2022

The reader need not examine this Python program too closely, as we will restate the algorithm, using the notation in this paper, elaborating on details, and analyzing the running time as we go.

3.1 Priorty Queue Algorithm

Let the interval $[x_{1},x_{2})$ be the input to the algorithm. For ease of notation, let $\Delta:=x_{2}-x_{1}$. As in the last section, the output is a histogram $h$ where $h[k]=\#\{n:x_{1}\leq n<x_{2}\mbox{ and }f(n)=k\}$.

It is fine to keep in mind $[1,x+1)$ as the input for now, but we are setting the stage for parallelization later.

3.1.1 Setup

A chain object is stored as $3$ integers: the smallest prime of the chain (pmin), the largest prime of the chain (pmax), and its sum (sum). We store chain objects in an array, and the array position indicates its length, so the length of the chain is stored implicitly, and array position zero is left empty.

To create this array of chain objects, we start by computing the longest chain $c$ that includes $2$, with maximal sum $<x_{2}$. Let $M$ be the length of $c$. Allocate the array ${\texttt{happy}}[\,]$ of length $M+1$, and set ${\texttt{happy}}[M]:=c$.

We then fill the array from position $M-1$ down to $1$. For long chains of length $m$, we copy the previous chain of length $m+1$, drop the largest prime, and then ${\mathbf{slide}}$ it as necessary until its sum is minimal, but $\geq x_{1}$. For short chains of length $m$, we estimate the largest prime of the chain to be $x_{1}/m$ (see Lemma 2.1), find consecutive primes below that to get $m$ of them in total, and then ${\mathbf{slide}}$ as necessary to get a sum that is minimal but $\geq x_{1}$.

We prove below that if we split long from short chains at $x_{2}^{1/3}$, the number of ${\mathbf{slide}}$ operations it takes to generate all the chains is $O(x_{2}^{2/3})$. This also tells us that we want $x_{2}-x_{1}\gg x_{2}^{2/3}$ for the setup time to take negligible time compared to our main loop. Unlike the puzzle algorithm, we do not need any space for the interval.

Lemma 3.1**.**

A list of chains, one of each length from $1$ to $M(x_{2})$, where the chain of length $m$ has minimal sum $\geq x_{1}$, can be constructed in a total number of ${\mathbf{slide}}$ operations that does not exceed $O(\sqrt{x_{2}/\log x_{2}}+x_{1}^{2/3})$.

Proof.

Let $m_{c}:=x_{1}^{1/3}$ be the cutoff length between long and short chains as described above.

The first chain of length $M(x_{2})$ can be made with $M(x_{2})=O(\sqrt{x_{2}/\log x_{2}})$ ${\mathbf{slide}}$ operations or the equivalent by Lemma 2.2.

Each subsequent long chain of length $m$ is created from the previous chain of length $m+1$ as described above, and this takes

[TABLE]

${\mathbf{slide}}$ operations, where $\pi(x)$ is the prime counting function. This telescopes when summed, giving

[TABLE]

by Lemma 2.1 and the prime number theorem.

A short chain of length $m$ needs $m$ primes near $x_{1}/m$, the equivalent of $O(m)$ ${\mathbf{slide}}$ operations. Summing over all short lengths $m$ gives $O(m_{c}^{2})$ total.

We then optimize the cost $O(M+x_{1}/m_{c}+m_{c}^{2})$ giving $m_{c}=x_{1}^{1/3}$ and our proof is complete. ∎

We conclude the setup by inserting the integers $1$ to $M$ into the priority queue, where the priority of the integer $i$ is ${\texttt{happy}}[i].{\texttt{sum}}$.

Later, we discuss two possible implementations of the priority queue: as a binary heap, or using a hash table. In either case, the initial queue can be constructed in $O(M)$ time: We can use a build operation for the heap, or $O(M)$ hash table insertions of chain objects, depending on our choice for data structure.

3.1.2 Main Loop

We describe the main loop of our algorithm below (Algorithm 3).

The inner While-loop, which starts at Line 3, iterates once per chain, which is once per gleeful representation, which is linear, or $O(\Delta)$ on average, over the course of the algorithm, due to Moser’s theorem. We have one delete-Min or dequeue operation, at most one insertion, and one ${\mathbf{slide}}$.

If we use a hash table to store the priority queue (where the key is the chain sum), insertions and deletions are constant time on average. If we use a standard binary heap, these operations take $O(\log M)$ time each. As $M<\sqrt{x_{2}}$, this is $O(\log x_{2})$, for a total of $O(\Delta\log x_{2})$ time for the interval.

Next, we discuss the cost of ${\mathbf{slide}}$ operations. If a list of primes is available, we can store the index of the prime in our chain objects, and then finding the next prime (for both pmax and pmin) takes constant time. However, we need primes as large as $x_{2}$, and finding all primes up to $x_{2}$ takes too much time, and storing them takes too much space. But we can pick a bound $B$, and compute and store all primes up to $B$. For primes needed beyond $B$, we can generate them either using an incremental prime sieve (see [12]) or using prime tests, with some support to make this as efficient as possible (see [11]).

Let us now derive optimal values for $B$, based on whether we use an incremental sieve or prime tests.

We begin with an incremental sieve. Galway [4] showed how to modify the Atkin-Bernstein sieve [2] to find primes in an interval $[a,b]$ using at most $O(b-a+b^{1/3})$ time and $O(b^{1/3})$ space. This sieve can, theoretically, be constructed to generate the primes as needed [12]. We divide the chains stored in our priority queue into those that use primes up to $B$, and those that use this incremental sieve to generate needed primes. Thus, each chain object that uses the sieve, would have attached to it its own copy of an Atkin-Bernstein-Galway incremental sieve to generate primes as needed for sliding. We now prove that the optimal value for $B$ is $x_{2}^{3/5}$ in this case.

Theorem 3.2**.**

After setup, the priority queue algorithm takes $O(\Delta\log x_{2}+x_{2}^{3/5})$ time to find and count all representations of integers in the interval $[x_{1},x_{2})$ where $\Delta=x_{2}-x_{1}$. The space used is $O(x_{2}^{3/5})$.

In practice, we want $\Delta\gg x_{1}^{2/3}$ so that the setup cost is negligible. Thus, for the interval $(1,x]$ we get a running time of $O(x\log x)$ using $O(x^{3/5})$ space. We can parallelize this keeping an optimal amount of work using up to $O(x^{1/3})$ processors, where each processor gets its own interval of length $x^{2/3}$.

Proof.

Set $B:=x_{2}^{3/5}$. The algorithm starts by finding all primes up to $O(B)$. Any linear-time prime sieve is fine for this. The primes can be stored in $O(B)$ bits.

By Lemma 3.1 the time for initializing the priority queue is negligible.

For ${\mathbf{slide}}$ operations on chains whose primes fall below $B$, ${\mathbf{slide}}$ operations take constant time. No extra space beyond the $O(B)$ to store the prime list is required. Let $m$ be the length of any such chain. Then $p(m,x_{2})\leq B$. By Lemma 2.1, we have $p(m,x_{2})=\Theta(x_{2}/m)$, giving $m\geq\Theta(x_{2}/B)=\Theta(x_{2}^{2/5})$.

Chains of length $m$ with $m\leq\Theta(x_{2}^{2/5})$, then, must have an incremental sieve attached. Each one uses space $O(p(m,x_{2})^{1/3})$. Again using Lemma 2.1, the total space used is proportional to

[TABLE]

For these shorter chains, ${\mathbf{slide}}$ operations take, on average, $O(\log x_{2})$ time each, as the Atkin-Bernstein-Galway sieve takes linear time to find primes up to a given bound, and primes are, on average, $O(\log x_{2})$ apart. ∎

If an overwhelming majority of chains are long, we might pull the average cost per ${\mathbf{slide}}$ below $\log x_{2}$, but this is not the case. The number of short chains (for the input $[1,x)$) is

[TABLE]

Given the running time, either data structure, binary heap or hash table, is reasonable for this approach.

Next, we consider generating primes as needed using a prime test. In particular, given a starting point $n$, we find the smallest prime larger than $n$ as follows:

• •

If $n<B$ we use the precomputed list of primes to find the next prime in constant time. This implies storing prime indices instead of prime values. We distinquish indices from prime values by storing indices as negative numbers (so the prime $p_{3}=5$ is stored as $-3$).

• •

For larger $n$ we follow the ideas in [11]:

Sieve an interval of size $2\log n$ using primes up to $\log n$. This takes $O(\log n)$ time leaving $O(\log n/\log\log n)$ possible primes in the interval $(n,n+2\log n]$. 2. 2.

Scan this interval for a possible prime, and then perform a base-$2$ strong pseudoprime test. 3. 3.

If the test is passed, the number is likely prime, so do a prime test. In theory, we use a pseudosquares prime test which takes $O((\log n)^{2})$ time in theory. In practice, we could use a Baille-PSW prime test, which works on 64-bit integers, or a Miller-Rabin style test on a fixed, short list of appropriate bases.

This is repeated as needed until a prime is found.

The average time to find a prime is dominated by the cost of the final prime test, at $O((\log n)^{2})$ time. In practice, we used a fixed wheel of size 30, performed some fixed trial divisions, and then used the Miller-Rabin test as described above, taking something closer to $O((\log n)^{2}/\log\log n)$ heuristic time.

We also need a previous prime function, but the cost is the same. We leave the details to the reader.

Theorem 3.3**.**

After setup, the priority queue algorithm takes $O(\Delta(\log x_{2})^{2})$ time to find and count all representations of integers in the interval $[x_{1},x_{2})$ where $\Delta=x_{2}-x_{1}$. The space used is $O(\sqrt{x_{2}\log x_{2}})$.

Proof.

After setup, the number of ${\mathbf{slide}}$ operations is $O(\Delta)$. From the discussion above, using the pseudosquares prime test we can perform each ${\mathbf{slide}}$ in $O((\log x_{2})^{2})$ time. The only significant use of space is to store the $O(M)$ chain objects, each of which uses $O(\log x_{2})$ bits. This completes the proof. ∎

4 Computational Results

Several different practical implementations of the algorithms presented in §2 and §3 were created and run over several upper bounds, and we present both their results and timing information below. Values of the histogram function $h(k,x)$, where $x$ ranges over powers of ten111The authors are aware that Tables 1 and 3 uses $x=10^{14}+145300$ in place of $x=10^{14}$. This will be corrected in a future revision of this paper., are given in Table 1. We also give the average density of $f$ on these intervals, in Figure 2.

The puzzle algorithm from §2 was implemented using Miller-Rabin deterministic primality testing, and ran in parallel on 560 AMD EPYC9734 cores at 2.20GHz. The interval $[10^{13},10^{14}]$ took 3 weeks, and the program needed 4 cumulative days to compute the interval $[1,10^{13}]$. We estimate that it would take on the order of one year to run this implementation on the interval $[10^{14},10^{15}]$, where we can heuristically expect to find the first integer with 15 gleeful representations.

The priority queue algorithm from §3 was implemented using the C++ standard library unordered_map as a hash table, and generated primes using an incremental sieve based on the sieve of Eratosthenes. This program ran on a single thread of an Intel Xeon E5-2650 v4 at 2.20GHz, running with an upper bound of $x=10^{12}$ and taking just over 26 days to terminate. Lower-order results include taking 3 days for $x=10^{11}$ and 6 hours cumulative for all smaller powers of ten.

We also created a version of the priority queue algorithm using a custom min-heap implementation with primality testing, storing chain objects indexed by length in an array. The authors discussed splitting the min-heap into two separate structures, one to store short chains and the other to store long chains. This was of interest, since one could cache the short-chain array for a possible run-time improvement. We can estimate the total work done in a heap by integrating the number of gleeful representations of a specific length over the lengths included in the heap, so the short heap would require

[TABLE]

amount of work with cutoff $k$. In Moser’s original paper, he used $k=\sqrt{x/\log x}$ as the bound on the maximum possible length of a gleeful number with sum $\leq x$ to get $x\log 2$ total representations. For the purposes of our work, if $k=(x/\log x)^{c}$ for some $c\in[0,1/2]$, then (4) simplifies to $-x\log(1-c+\frac{c\log\log x}{\log x})$. We want to have reasonable hope of a cache hit for our short chains, where the proportion of work done in cache is $-\log(1-c+\frac{c\log\log x}{\log x})/\log 2$. For illustration, if we wanted to get a cache hit on half of every update cycle when running with upper bound $x=10^{15}$, we would need to hold $x^{0.326}$ (where $c=(1-2^{-0.5})/(1-\frac{\log\log 10^{15}}{\log 10^{15}})$) chains, about 78000, in cache, well beyond the scope of our machines. For this reason, the separate heap structures were never implemented.

4.1 Moser’s Questions

Recall the questions Moser posed in his 1963 paper:

Since the average value of $f(n)$ is $\log 2$ it follows that $f(n)=0$ infinitely often. The following problems, among others, suggest themselves:

Is $f(n)=1$ infinitely often? 2. 2.

Is $f(n)=k$ solvable for every $k$? 3. 3.

Do the numbers for which $f(n)=k$ have a density for every $k$? 4. 4.

Is $\limsup_{n\to\infty}f(n)=\infty$?

In the spirit of Cramér’s model on the distribution of primes, the frequencies in Table 1 give rise to a conjecture on a probabilistic model for Moser’s function $f$. Namely, although not a true random probability distribution, $f(n)$ is roughly Poisson with mean $\lambda=\log 2$. Compare the histogram data for $x=10^{14}$ (restated for convienence) with the expected counts of a Poisson variable (Table 3). In this spirit, we can heuristically answer in the affirmative to Moser’s four questions.

In particular interest to Moser’s second and fourth questions, OEIS A054859 lists the smallest integer having exactly $k$ representations. Prior to our work, this sequence was known to $k=13$; we present the first known integer with 14 representations as a sum of consecutive primes. Table 4 gives the complete updated OEIS table, and Table 5 gives the 14 representations for our submitted integer. It is particularly interesting to note that 11 of the 14 records in Table 4 are themselves prime.

With regard to Moser’s first and third questions, we can conjecture that the density of numbers $n$ for which $f(n)=k$ follows the Poisson distribution’s density: asymptotically $\frac{(\log 2)^{k}}{2\cdot k!}$. Indeed, we expect half of the integers will have no representation as a sum of consecutive primes, $\log 2/2$ (roughly 1/3) of the integers will have exactly one representation, and so forth.

5 Acknowledgements

This work was supported in part by a grant from the Holcomb Awards Committee. We thank Frank Levinson for his generous support of Butler’s research computing infrastructure.

Thanks to Ankur Gupta and Jonathan Webster for several helpful discussions.

The second author was an undergraduate student when this work was done.