## Find the number of $n$ such that $n!$ is a sum of three squares

6

1

I want to check the following theorem by using Mathematica:

$\textbf{Theorem}$. $\text{The estimate}$

$\# \{n \le x:n! \text{ is a sum of three squares}\}=7x/8+O(x^{2/3})$

$\text{holds.}$

(The positive integers $n$ such that $n!$ is a sum of three squares have a density which is equal to $7/8$.)

I tried:

Solve[{n! == x^2 + y^2 + z^2}, {x, y, z, n}, Integers]


but there is an error message.

I want to do something like this (find the all perfect number from $1$ to $1000$)

Select[Range, DivisorSigma[1, #] == 2 # &]


How do I find the number of $n$ such that $n!$ is a sum of three squares?

1Make your search space bounded: Solve[{n! == x^2 + y^2 + z^2, x <= y <= z, Sequence @@ (0 < # < 70 & /@ {x, y, z})}, {x, y, z, n}, Integers] (* {{x -> 1, y -> 1, z -> 2, n -> 3}, {x -> 2, y -> 2, z -> 4, n -> 4}, {x -> 2, y -> 4, z -> 10, n -> 5}, {x -> 4, y -> 20, z -> 68, n -> 7}, {x -> 8, y -> 16, z -> 20, n -> 6}, {x -> 12, y -> 36, z -> 60, n -> 7}, {x -> 20, y -> 44, z -> 52, n -> 7}} *) - this is not necessarily fast, but it works. – kirma – 2016-04-24T11:34:42.130

1@vito This was a neat problem. Where did you see this theorem btw? – Chip Hurst – 2016-04-26T02:40:40.690

@ChipHurst "The Legacy of Alladi Ramakrishnan in the Mathematical sciences". page 243 – vito – 2016-04-26T07:34:34.870

4

If you only care about counting and enumerating, why not use Legendre's three-square theorem?

SetAttributes[SumOf3SquaresQ, Listable];

SumOf3SquaresQ[n_] := Mod[n/4^IntegerExponent[n, 4], 8] != 7

FactorialSumOf3SquaresPi[x_] := Total[Boole[SumOf3SquaresQ[Range[x]!]]]

FactorialSumOf3SquaresPi

8746


Or test your conjecture by subtracting $7x/8$ and dividing by $x^{2/3}$:

With[{x = 1000},
ListLinePlot[
(Accumulate[Boole[SumOf3SquaresQ[Range[x]!]]] - 7 Range[x]/8)/Range[x]^(2/3)
]
] Or list such numbers:

Pick[Range, SumOf3SquaresQ[Range!]]

{1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44,
45, 46, 47, 50}


## Edit

If we keep a running tally (instead of computing things separately), we can get results much quicker:

FactorialSumOf3SquaresPiFast = Compile[{{x, _Integer}},
Module[{count = 0, m, offset = 1, prod = 1},
Do[
m = n;
While[BitAnd[m, 3] == 0, m = Floor[m/4]];
If[EvenQ[m], m = Floor[m/2]; offset++];
prod = BitAnd[prod * m, 7];
If[prod != 7 - Boole[EvenQ[offset]], count++],
{n, 1, x}
];
count
],
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];


Compare timings of old vs new:

FactorialSumOf3SquaresPi // AbsoluteTiming

{1.57493, 8746}

FactorialSumOf3SquaresPiFast // AbsoluteTiming

{0.000495, 8746}


Here's timings on large $x$:

FactorialSumOf3SquaresPiFast[10^7] // AbsoluteTiming

{0.488944, 8751280}

FactorialSumOf3SquaresPiFast[10^8] // AbsoluteTiming

{4.8907, 87502143}

FactorialSumOf3SquaresPiFast[10^9] // AbsoluteTiming

{49.8406, 875001736}


It's annoying that BitShiftRight[] is still uncompilable. Otherwise, nice! BTW: I'd use Quotient[m, 2] instead of Floor[m/2]. – J. M.'s ennui – 2016-04-25T18:23:49.737

@J.M. I originally used Quotient, but this way shaved off some time actually. Yeah, too bad BitShiftRight isn't compilable. – Chip Hurst – 2016-04-25T18:36:39.090

Huh, would've thought keeping everything in integers would be faster. Oh well. – J. M.'s ennui – 2016-04-25T19:02:35.677

12

Brute-force, but compact:

DeleteCases[Table[{k, PowersRepresentations[k!, 3, 2]}, {k, 10}], {___, 0, ___}, {3}]

{{1, {}}, {2, {}}, {3, {{1, 1, 2}}}, {4, {{2, 2, 4}}}, {5, {{2, 4, 10}}},
{6, {{8, 16, 20}}}, {7, {{4, 20, 68}, {12, 36, 60}, {20, 44, 52}}},
{8, {{8, 16, 200}, {8, 80, 184}, {40, 88, 176}, {72, 120, 144}, {80, 104, 152}}},
{9, {{8, 304, 520}, {16, 200, 568}, {24, 48, 600}, {24, 240, 552}, {40, 104, 592},
{40, 272, 536}, {80, 184, 568}, {80, 344, 488}, {120, 264, 528}, {136, 272, 520},
{152, 400, 424}, {176, 248, 520}, {184, 368, 440}, {200, 328, 464},
{216, 360, 432}, {240, 312, 456}, {248, 376, 400}}}, {10, {}}}


For larger factorials, it might take quite longer; e.g. Length[PowersRepresentations[11!, 3, 2]] == 126.

I wasn't aware of PowersRepresentations, good to know! (Assuming one actually finds other applications than snarky Mma.SE answers for it...) – kirma – 2016-04-24T13:18:43.250

2It's a useful decomposition sometimes, but it does become combinatorically prohibitive very quickly. – J. M.'s ennui – 2016-04-24T13:21:52.813

@J.M. +1. but how to remove empty brackets from your result (list)? I want to find only number of $n$, such $n!$ is a sum of three squares – vito – 2016-05-29T15:58:03.327

1Just modify the filtering condition in DeleteCases[], @vito: DeleteCases[(* stuff *), {_, {{___, 0, ___} ..} | {}}] – J. M.'s ennui – 2016-05-29T18:15:36.780

@J.M. I tried: Length[DeleteCases[Table[{k, PowersRepresentations[k!,3, 2]}, {k, N}], {_, {{___, 0, ___} ..} | {}}]] it works perfectly for $N=11$, but for $N=20$ it doesn't produce a result even after a few hours of running :) – vito – 2016-05-29T19:16:58.947

@vito, yes, that's what I meant by "combinatorically prohibitive". – J. M.'s ennui – 2016-05-29T19:26:52.113

7

You can Try

Table[set = {x, y, z} /.
NSolve[{n! == x^2 + y^2 + z^2, x > 0, y > 0, z > 0},{x, y, z},Integers];
set = Union[Sort /@ set];
Join[{n}, set],
{n, 3,8}]


This will give you how you can express a factorial as a sum of three squares. Now you can use further conditions (like $x\neq y \neq z$) with Select or Cases to filter them.