6
$\begingroup$

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[1000], DivisorSigma[1, #] == 2 # &]

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

$\endgroup$
3
  • 1
    $\begingroup$ Make 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. $\endgroup$
    – kirma
    Commented Apr 24, 2016 at 11:34
  • 1
    $\begingroup$ @vito This was a neat problem. Where did you see this theorem btw? $\endgroup$
    – Greg Hurst
    Commented Apr 26, 2016 at 2:40
  • $\begingroup$ @ChipHurst "The Legacy of Alladi Ramakrishnan in the Mathematical sciences". page 243 $\endgroup$
    – vito
    Commented Apr 26, 2016 at 7:34

3 Answers 3

4
$\begingroup$

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[10000]
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)
  ]
]

enter image description here

Or list such numbers:

Pick[Range[50], SumOf3SquaresQ[Range[50]!]]
{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[10000] // AbsoluteTiming
{1.57493, 8746}
FactorialSumOf3SquaresPiFast[10000] // 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}
$\endgroup$
3
  • $\begingroup$ It's annoying that BitShiftRight[] is still uncompilable. Otherwise, nice! BTW: I'd use Quotient[m, 2] instead of Floor[m/2]. $\endgroup$ Commented Apr 25, 2016 at 18:23
  • $\begingroup$ @J.M. I originally used Quotient, but this way shaved off some time actually. Yeah, too bad BitShiftRight isn't compilable. $\endgroup$
    – Greg Hurst
    Commented Apr 25, 2016 at 18:36
  • $\begingroup$ Huh, would've thought keeping everything in integers would be faster. Oh well. $\endgroup$ Commented Apr 25, 2016 at 19:02
12
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ I wasn't aware of PowersRepresentations, good to know! (Assuming one actually finds other applications than snarky Mma.SE answers for it...) $\endgroup$
    – kirma
    Commented Apr 24, 2016 at 13:18
  • 2
    $\begingroup$ It's a useful decomposition sometimes, but it does become combinatorically prohibitive very quickly. $\endgroup$ Commented Apr 24, 2016 at 13:21
  • $\begingroup$ @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 $\endgroup$
    – vito
    Commented May 29, 2016 at 15:58
  • 1
    $\begingroup$ Just modify the filtering condition in DeleteCases[], @vito: DeleteCases[(* stuff *), {_, {{___, 0, ___} ..} | {}}] $\endgroup$ Commented May 29, 2016 at 18:15
  • $\begingroup$ @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 :) $\endgroup$
    – vito
    Commented May 29, 2016 at 19:16
7
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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