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$
  • 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 Apr 24 '16 at 11:34
  • 1
    $\begingroup$ @vito This was a neat problem. Where did you see this theorem btw? $\endgroup$ – Chip Hurst Apr 26 '16 at 2:40
  • $\begingroup$ @ChipHurst "The Legacy of Alladi Ramakrishnan in the Mathematical sciences". page 243 $\endgroup$ – vito Apr 26 '16 at 7:34
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$
  • $\begingroup$ It's annoying that BitShiftRight[] is still uncompilable. Otherwise, nice! BTW: I'd use Quotient[m, 2] instead of Floor[m/2]. $\endgroup$ – J. M. will be back soon Apr 25 '16 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$ – Chip Hurst Apr 25 '16 at 18:36
  • $\begingroup$ Huh, would've thought keeping everything in integers would be faster. Oh well. $\endgroup$ – J. M. will be back soon Apr 25 '16 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$
  • $\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 Apr 24 '16 at 13:18
  • 2
    $\begingroup$ It's a useful decomposition sometimes, but it does become combinatorically prohibitive very quickly. $\endgroup$ – J. M. will be back soon Apr 24 '16 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 May 29 '16 at 15:58
  • 1
    $\begingroup$ Just modify the filtering condition in DeleteCases[], @vito: DeleteCases[(* stuff *), {_, {{___, 0, ___} ..} | {}}] $\endgroup$ – J. M. will be back soon May 29 '16 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 May 29 '16 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, privacy policy and cookie policy

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