# In how many ways can a number be decomposed in a sum of distinct perfect squares

Is there a formula for how many ways can a number be decomposed in sum of distinct perfect squares.

Example: $$50 = 49 + 1 = 1 + 4 + 9 + 36 = 9 + 16 + 25$$

The number 50 can be decomposed in 3 different ways.

• If you are talking about distinct squares, then not every random number can be decomposed into squares. For example, $6$ cannot be decomposed into perfect squares. – Devansh Kamra Jan 23 at 17:54
• Yes, I meant perfect squares. I also edited the info. @DevanshKamra – Bili Debili Jan 23 at 17:58
• also 50 = 25+25 – Lozenges Jan 23 at 17:58
• @Lozenges They have to be distinct, I will also edit that in the info – Bili Debili Jan 23 at 17:59

(Edit: Added a plot of $$a_n$$ and its asymptotics, corrected numerics. Also, a short answer: There's not really a formula, but people can work out its growth.)

Denote $$a_n =$$ number of ways to write $$n$$ as a sum of distinct squares.

Here we demonstrate, though not quite a formula (as you might expect), how we can compute $$a_n$$ by using the generating function for $$a_n$$. Indeed, we have $$A =\sum_{n=0}^{\infty} a_n q^n = \prod_{k=1}^\infty(1+q^{k^2})=(1+q)(1+q^4)(1+q^9)(1+q^{16})\cdots \qquad(\dagger)$$ (Why does this work? Pretend you are expanding this product and keep track of the powers of $$q$$ that you collected.)

You can play around to see what these values are (choose a large enough $$N$$) with SageMath (https://sagecell.sagemath.org/):

var('q')
N = 10
A = 1
for k in range(1,N):
A = A*(1+q^(k^2))
A.expand()

This computes the first $$N=10$$ product in $$(\dagger)$$, and gives

...+ 4*q^65 + q^64 + q^63 + 3*q^62 + 2*q^61 + q^59 + q^58 + q^57 + q^56 + q^55 + 2*q^54 + 2*q^53 + q^52 + q^51 + 3*q^50 + 2*q^49 + 2*q^46 + 2*q^45 + q^42 + 2*q^41 + q^40 + q^39 + q^38 + q^37 + q^36 + q^35 + q^34 + 2*q^30 + 2*q^29 + 2*q^26 + 2*q^25 + q^21 + q^20 + q^17 + q^16 + q^14 + q^13 + q^10 + q^9 + q^5 + q^4 + q + 1

(truncated)

Note the coefficient in front of $$q^{50}$$ is $$3$$, so there are three such decompositions for 50. As another example, the coefficient of $$q^{65}$$ is 4, so we have four such decompositions for 65. The term $$q^{60}$$ is missing, this means there is no such decomposition for 60.

In any case, you can compute $$a_n$$ by using its generating function $$A$$ and look at the coefficients of $$q^n$$ in $$A$$. (And Use a suitably large $$N>\sqrt n$$, computationally, so you don't miss any decompositions. See update edit below for a SageMath implementation of this.)

Some observation is that it seems quite "sparse", these decompositions are rare as the power of $$q$$ in the factor $$(1+q^{k^2})$$ grows quite quickly. However this does not imply every integer $$n$$ has small $$a_n$$. Indeed, as $$n$$ increases so will $$a_n$$. For instance $$a_{1000} = 1269,a_{10000} = 3296089777$$, and $$a_{100000} = 2759256389896728737285379$$. Here is a plot of $$a_n$$ (in $$\color{purple} {\rm purple}$$) for $$n\le 1000$$ using SageMath (code to plot this at the end):

Now, with these several terms of $$a_n$$, we can search on OEIS which gives this reference here: http://oeis.org/A033461 . There doesn't seem to have a nice closed formula for $$a_n$$, but people have estimated its asymptotic growth (Vaclav Kotesovec, Dec 09 2016):

$$a_n \sim \exp(3 * 2^{-5/3} \pi^{1/3} ((\sqrt{2}-1)\zeta(3/2))^{2/3} n^{1/3}) \frac{(\sqrt{2}-1)\zeta(3/2))^{1/3} }{ 2^{4/3} \sqrt{3} \pi^{1/3} n^{5/6}},$$ where $$\zeta(3/2)\approx 2.612375$$ is the Zeta function at $$3/2$$ . This asymptotic is ploted in $$\color{yellow} {\rm yellow}$$ above. Such an asymptotic estimate is often found by doing some (complex) analysis on the generating function $$A(q)$$ mentioned above (a typical method is called saddle-point method. A great reference is Analytic Combinatorics by Philippe Flajolet and Robert Sedgewick). I do not know the exact detail of this particular asymptotic, however.

An simple question to ask is: Is $$a_n\neq 0$$ for all sufficiently large $$n$$? Or are there "gaps", that $$a_n = 0$$ infinitely often? One would guess it is the former by looking at the plot above, and indeed it is true:

Theorem. The number $$128$$ is the largest integer that cannot be written as a sum of distinct squares. And a complete list of positive integers $$n$$ that cannot be written as a sum of distinct squares are: $$2, 3, 6, 7, 8, 11, 12, 15, 18, 19, 22, 23, 24, 27, 28, 31, 32, 33,\\ 43, 44, 47, 48, 60, 67, 72, 76, 92, 96, 108, 112, 128.$$

Proofs and references of this can be found here: https://proofwiki.org/wiki/Numbers_not_Sum_of_Distinct_Squares#:~:text=Then%20k%20can%20be%20expressed,distinct%20squares%20by%20Induction%20Hypothesis and here: Numbers which are not the sum of distinct squares .

Postscript. Even though we didn't get a "formula" as you might traditionally expect, this generating function $$A$$ encodes all the information of these distinct square decompositions, to which we can extract the appropriate coefficient. In fact, we do have (formally) $$a_n = \frac{A^{(n)}(0)}{n!},$$ where $$A^{(n)}(0)$$ is the $$n$$-th derivative of $$A(q)$$ evaluated at $$q=0$$. In this vein we have to ask, "what is an answer" and "what is formula"...

Update. One can easily fashion this SageMath code to compute $$a_n$$ (though not the most efficient):

var('q')
n = 1000 # choose your value n to output a_n = number of decomp of n into distinct squares.
N = 1+ceil(sqrt(n))
A = 1
for k in range(1,N):
A = A*(1+q^(k^2))
A.coefficients(sparse=False)[n]

---
output:
1269

Here $$n=1000$$, and we get $$a_{1000} = 1269$$. Also for reference $$a_{10000} = 3296089777$$ and $$a_{128} = 0$$. For a more efficient algorithm you can search for subset sum algorithms (say https://stackoverflow.com/questions/18305843/find-all-subsets-that-sum-to-a-particular-value).

Added. SageMath to plot the first few $$a_n$$ values and its asymptotic estimate:

var('q')
n = 1000 # choose your value n to output a_n = number of decomp of n into distinct squares.
N = 1+ceil(sqrt(n))
A = 1
for k in range(1,N):
A = A*(1+q^(k^2))
w = A.coefficients(sparse=False)
z = 2.6123753486854883433485675679240716305708006524000634075733282488149277676882728609962438681263119523829763587721497556981576329
h = list_plot([w[i] for i in range(1,n)],plotjoined=True, color='purple')+list_plot([exp(3 * 2^(-5/3) * pi^(1/3) * ((sqrt(2)-1)*z)^(2/3) * x^(1/3)) * ((sqrt(2)-1)*z)^(1/3) / (2^(4/3) * sqrt(3) * pi^(1/3) * x^(5/6)) for x in range(1,n)],plotjoined=True, color='yellow')
show(h)
• How far above $128$ did you search looking for numbers that cannot be represented? If the range is enough it is often possible to extend it into a proof, but I don't know how much is needed. If you had checked to $10,000$, say, you might be able to show that any higher number can be brought into the range $129-10,000$ using no squares less than $101^2$. You could then add the squares that bring you to that range to the solution for whichever number you arrive at. – Ross Millikan Jan 24 at 15:22
• You're right, such a proof actually exists that I just found, proofwiki.org/wiki/…. I will add edit it in. – bonsoon Jan 24 at 21:25
• And I see, a reference of yours here! math.stackexchange.com/questions/1018801/… – bonsoon Jan 24 at 21:32
• I had forgotten I had done it – Ross Millikan Jan 24 at 22:25