3
$\begingroup$

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.

$\endgroup$
4
  • 1
    $\begingroup$ 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. $\endgroup$ – Devansh Kamra Jan 23 at 17:54
  • $\begingroup$ Yes, I meant perfect squares. I also edited the info. @DevanshKamra $\endgroup$ – Bili Debili Jan 23 at 17:58
  • 1
    $\begingroup$ also 50 = 25+25 $\endgroup$ – Lozenges Jan 23 at 17:58
  • $\begingroup$ @Lozenges They have to be distinct, I will also edit that in the info $\endgroup$ – Bili Debili Jan 23 at 17:59
7
$\begingroup$

(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):

enter image description here

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)
$\endgroup$
4
  • 1
    $\begingroup$ 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. $\endgroup$ – Ross Millikan Jan 24 at 15:22
  • 1
    $\begingroup$ You're right, such a proof actually exists that I just found, proofwiki.org/wiki/…. I will add edit it in. $\endgroup$ – bonsoon Jan 24 at 21:25
  • 1
    $\begingroup$ And I see, a reference of yours here! math.stackexchange.com/questions/1018801/… $\endgroup$ – bonsoon Jan 24 at 21:32
  • 1
    $\begingroup$ I had forgotten I had done it $\endgroup$ – Ross Millikan Jan 24 at 22:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.