Update on 6/19/2020. This discussion led to deeper and deeper results on the topic. The last findings are described in my new post (including my two answers), here.

I came up with the following conjecture. All non-square integers $z$ can be represented as $z=x^2 + y$ where $x$ is an integer and $y$ is a prime. The exceptions are

z = 10, 34, 58, 85, 91, 130, 214, 226, 370, 526, 706, 730, 771, 1255, 1351, 1414, 1906, 2986, 3676, 9634, 21679.

Note that this is deeper than Goldbach conjecture (all even numbers are the sum of two primes) because squares are far rarer than primes. Also, few numbers are the sum of two squares, such numbers (sums of two squares) are far more abundant than primes, but their natural density is also zero. But all numbers are the sum of four squares. Surprisingly, all integers can be represented as $z = \lfloor x^c \rfloor + \lfloor y^c \rfloor$ where $x, y$ are positive integers and $c < \log_{22} 63$ is a positive constant; but this fails at $c = \log_{22} 63$ as $z=73$ becomes an exception. See section 1 in this article for details; this is also a conjecture.

Question: Can you verify if my conjecture is true up to some very large $z$? I tested it only for $0\leq z < 750000$.

Heuristics behind this conjecture

This is by no mean a proof, but rather, I explain here why I think it could be true. Let denote as $r(z)$ the number of solutions to $x^2 +y \leq z$ where $x, y$ are integers and $y$ is prime. For a fixed large $z$, we want to count the number of integer couples $(x, w)$ below the curve $z=x^2+ w\log w$, with $x, w\geq 0$, in order to approximate $r(z)$. The choice of $w \log w$ is a direct consequence of the prime number theorem, replacing primes by their approximation, for large primes. That count $r(z)$ grows faster than $O(z)$. The derivative $dr(z)/dz$ thus grows faster than $O(1)$, and it shows how the number of solutions to $z=x^2+y$ grows on average, faster than $O(1)$ as $z$ increases.

More details about the heuristic approach

Essentially, we are trying to count the number of blue points under the red curve in the plot below (in this example, $z=100$). The equation for the curve is $w \log w = z-x^2$, and $z$ is assumed to be fixed.

enter image description here

The equation can be re-written as $w = (z-x^2)/W(z-x^2)$ where $W$ is the Lambert function, which behaves asymptotically like the $\log$ function. Thus the number of points below the red curve is asymptotically (for large values of $z$) equal to

$$r(z) \sim \int_0^\sqrt{z} \frac{z-x^2}{W(z-x^2)}dz \sim \int_0^\sqrt{z} \frac{z-x^2}{\log(z-x^2)}dz = \frac{1}{2}\int_0^z \frac{u}{\sqrt{z-u}\cdot\log u}du.$$

Let us denote as $\phi(z)$ the function defined by the rightmost integral. We have $r(z) \sim \phi(z)$. I computed the exact values of $r(z)$ and $\phi(z)$ for various small and large $z$, and clearly, $r(z) \rightarrow C \cdot \phi(z)$, but I am not sure if $C=1$. See WolframAlpha computations here.

The number of solutions to $z=x^2+y$ (with $y$ prime) is thus, on average, as $z$ gets larger and larger, asymptotically equivalent to $d\phi(z) / dz$. Below is a table featuring $r(z)$ and $\phi(z)$.

enter image description here

Good asymptotic approximations for very large $z$ are

$$\phi(z)\approx\frac{2}{3}\cdot \frac{z^{3/2}}{\log z} \mbox{ and } \frac{d\phi(z)}{dz}\approx \frac{\sqrt{z}}{\log z}.$$

The last result is compatible with the one posted in the answer by Dietrich Burde, confirming that the approach I used here is sound. Note that the same methodology could be applied to sums of squares or sums of primes or any sums of integers. It is pretty generic.

Final comment

The number of solutions to $z = x^2 + y$ (with $y$ prime, $x$ an integer) is equal to $r(z)-r(z-1)$. In all cases, $r(z)$ grows slowly (polynomial at most) and thus $r(z)-r(z-1) \sim dr(z)/dz$. We could get deeper results with second- and third-order approximations in all the asymptotic results used in this article, rather than just first-order approximations.

Below is a chart featuring the distribution for the number of solutions to $z=x^2+y$ [that is, the distribution of $r(z)-r(z-1)$] for $700000\leq z < 740000$. enter image description here

For instance, there are $441$ different $z$'s between $z = 700000$ and $z = 740000$ for which $z=x^2 + y$ has exactly $50$ solutions. Below is the same chart, but for $100000\leq z < 140000$. The two distributions are strikingly similar in shap2. enter image description here

Finally, among the first 750,000 $z$'s, we have:

  • $z = 78754$ is the last one to admit only one decomposition as $z = x^2+y$
  • $z = 101794$ is the last one to admit exactly two decompositions
  • $z = 339634$ is the last one to admit exactly three decompositions
  • $z = 438166$ is the last one to admit exactly four decompositions
  • $z = 383839$ is the last one to admit exactly five decompositions

The $z$'s that admit only one decomposition are listed below. I searched for this sequence to see if it had been discovered, but could not find any reference.

z = 2, 5, 8, 13, 15, 22, 24, 26, 31, 37, 40, 46, 50, 55, 61, 70, 74, 76, 82, 94, 99, 106, 115, 120, 127, 133, 136, 142, 145, 154, 159, 166, 170, 178, 184, 202, 205, 219, 221, 235, 246, 250, 253, 265, 268, 274, 295, 298, 301, 310, 316, 319, 325, 328, 334, 340, 346, 379, 391, 394, 399, 412, 424, 436, 439, 442, 445, 469, 490, 505, 511, 559, 562, 571, 574, 586, 589, 610, 616, 646, 694, 781, 793, 799, 829, 834, 835, 874, 914, 922, 946, 949, 970, 979, 991, 994, 1030, 1045, 1066, 1090, 1105, 1164, 1204, 1219, 1243, 1324, 1354, 1366, 1384, 1411, 1450, 1501, 1549, 1555, 1642, 1717, 1726, 1765, 1786, 1810, 1885, 1981, 1990, 2041, 2059, 2074, 2146, 2167, 2245, 2266, 2284, 2344, 2410, 2416, 2479, 2650, 2806, 2821, 2854, 2899, 2926, 3004, 3094, 3151, 3166, 3184, 3319, 3418, 3502, 3811, 3859, 3865, 3964, 3991, 4216, 4222, 4279, 4330, 4414, 4504, 4510, 4645, 4654, 4711, 4930, 5482, 5506, 5545, 5986, 6031, 6049, 6274, 6439, 7009, 7081, 7441, 7549, 7954, 8086, 8584, 8824, 9214, 9571, 10165, 10774, 11509, 11806, 13834, 15106, 15334, 15565, 16081, 16186, 23851, 31879, 33205, 44536, 78754
  • 3
    $\begingroup$ “Stronger” is the wrong word here, because “A is stronger than B” means “A implies B.” I see why you’d use the term, but it is not the correct word here. $\endgroup$ Jun 7, 2020 at 18:58
  • $\begingroup$ Can you suggest a better term? I agree with you. Also I inadvertently posted it before completing my post (I thought you could save as a draft). I'll be adding more stuff soon.. $\endgroup$ Jun 7, 2020 at 19:01
  • $\begingroup$ Using the word "deeper" instead. $\endgroup$ Jun 7, 2020 at 19:04

2 Answers 2


This is Hardy and Littlewood's Conjecture $H$. It says that this sequence $a(n)= 10,34,58,85,\ldots$ is finite and that the number of representations of $n$ as the sum of a prime and a square is asymptotically $$ \frac{\sqrt{n}}{\log (n)} \cdot \prod_{p > 2}\left( 1 - \frac{(n / p)}{p - 1}\right)$$

where $(n / p)$ is the Legendre symbol.

References: Mikawa, Narkiewicz, Wang

The conjecture is tested up to $10^{11}$ so far, i.e., it is known that $a(22) > 10^{11}$, if it exists.

  • $\begingroup$ Here is a link to the sequence in question: oeis.org/A020495 $\endgroup$ Jun 7, 2020 at 19:21
  • 1
    $\begingroup$ There is also book The Hardy-Littlewood Method, by R. C. Vaughan. $\endgroup$
    – Will Jagy
    Jun 7, 2020 at 20:38
  • $\begingroup$ It's interesting that it hasn't been verified further. It should be possible to test all numbers in that range in less than an hour on a normal computer (checking 240 squares will find a prime in all but one in a million cases, and that can be done very efficiently). $\endgroup$
    – gnasher729
    Jun 12, 2020 at 11:14

I found it interesting that this has been verified only up to $10^{11}$, so I wrote a small(ish) program veryifying the conjecture up to $10^{11}$ in just 400 seconds, and up to $1.6 \cdot 10^{13}$ in under one day. At that size it checks 1 billion numbers in about 5 seconds on a single three year old, not top-of-the-range computer.

The algorithm that I used:

Create a bitmap representing odd primes. Bit #i in that bitmap is set if 2i+1 is not a prime, and cleared if 2i+1 is a prime. This is stored as an array of 64 bit integers, so 64 potential primes can be processed in constant time. Increase the size of the bitmap whenever more primes are needed, since that number is hard to product.

Choose an integer w, and then examine whether the integers in [kw, (k+1)w) contain any counter examples to the conjecture, for k = 0, 1, 2 etc. Since all primes other than 2 are odd, we examine even integers and odd integers separately. (Even integers cannot be the sum of an even square and odd prime obviously, and odd integers cannot be the sum of an odd square and an odd prime).

To do this, we create a bitmap representing either the even or the odd integers in the range, and set a bit in the bitmap for each integer that might be a counterexample, and clear the bit if that integer is not a counterexample. Bit #i either represents kw + 2i, or kw + 2i + 1. We let S = largest integer such that $S^2 < (k+1)w$.

Initially all bits are set. Then we let s be the largest even/odd integer <= S, and clear the bit for $s^2$ (squares are excluded) and $s^2+2$ (sum of square and prime), plus we examine s-2, s-4 etc as well as long as their square is in the range.

Then we start with s = largest odd / even integer <= S. We can then use the bitmap of primes to very efficiently remove all sums $s^2 + p$ from the bitmap of potential counterexamples, for 64 integers at a time. Then we proceed with s-2, s-4, s-6 etc. until there are no counterexamples left or if the next s would be negative. And that's basically it; if you choose w as a multiple of 128 (to keep bit operations simple) and not too large (should fit into the first level cache of the computer, I picked something around 200,000), then this will run at a speed of a billion integers tested in a few seconds.

If you want to go further, there are two optimisations that I didn't bother to implement. First optimisation: For the first say 100 squares, most 64 bit words will contain one or more bits representing counter examples, but at some point these words get rare. At that point it will be faster to keep track of which 64 bit words represent any counter examples, and only inspect the primes for those words. This is especially beneficial when there is only a dozen or so counter examples left. My estimate is that this should make the algorithm 3 times faster.

The second optimisation is much more involved. As the algorithm proceeds, you can see that it consistently uses all the primes, and each one exactly once (for example if you examine integers around $10^{12}$, even squares are about 4 million apart, so the same primes are not reused). This means that these primes will not be cached but must be read from RAM each time, which is the worst thing we can do.

To avoid this, we need to change the order in which we examine numbers so that the same primes will be used repeatedly. We take w much smaller. Then we examine integers in a range ($s^2 + kw$, $s^2 + (k+1)w$ for k = 0, 1, 2, etc. and $s^2$ being consecutive squares. We skip ranges that would be covered by a larger s. The primes that are examined will be about kw to (k+1)w, then (4s + kw) to (4s + (k+1)w etc. When we examine the range ($(s+2)^2 + kw$ to $(s+2)^2 + (k+1)w$) we examine primes in the same range. If we keep w small enough that all these primes fit into some cache, the time to read them will be much faster. This could be an improvement from 3 to 10 times.

This also makes it worthwhile using multiple cores on a processor, making another substantial improvement possible with a much more expensive computer.

Now there is a HUGE problem with all of this: I haven't actually verified anything. I wrote a program that printed the exceptions that were suggested, and then printed whenever it didn't find more exceptions. However, to make sure that it actually did anything meaningful, someone would have to verify the source code very, very carefully to make sure that the output of the program actually verifies anything. And with some paranoia, you'd have to verify either the compiler or the compiled code.

For this particular problem, it is impossible to output anything that demonstrates the verification of the conjecture is correct. (Of course I could print how each number up to 16 trillions is the sum of a square and a prime, but that's practically impossible to verify).

PS. It seems that 78526384 might be the largest integer which is not a sure of a cube and a prime. Finding the largest integer which is not the sum of a fourth power and a prime turns out to be very difficult: First, the numbers involved are obviously a lot larger, because fourth powers are much more rare than squares or cubes. But there is another problem: It is always the case that a^4 modulo 10 = 1 or 6, except when a = 10 modulo 0 or 5. So if n = 1 modulo 10 or 6 modulo 10, then of the 5 even or odd fourth powers that we could subtract from n, four lead to a result ending in 0 or 5. Which makes it very much rarer that the difference is a prime.

  • $\begingroup$ Can you share the source code? I want to modify reuse the optimized code for a different computation but somewhat similar that I am doing $\endgroup$ Mar 14, 2023 at 2:02

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .