This is a front-end for the Online Encyclopedia of Integer Sequences, made by Christian Perfect. The idea is to provide OEIS entries in non-ancient HTML, and then to think about how they're presented visually. The source code is on GitHub.
%I A014597 #65 Jun 17 2025 10:27:54 %S A014597 1,3,5,11,12,27,29,71,72,213,215,603,635,1917,1183893 %N A014597 Numbers k such that k^2 is a sum of distinct factorials. %C A014597 a(16)^2 > 48! (about 1.24139*10^61), if it exists. - _Jon E. Schoenfield_, Aug 04 2006 %C A014597 A197183(a(n)) = 1. - _Reinhard Zumkeller_, Dec 04 2011 %C A014597 a(16) > 4.3*10^55 if it exists. - _Bert Dobbelaere_, Sep 16 2020 %D A014597 Posting by _Dan Hoey_ to math-fun mailing list. %H A014597 Shyam Sunder Gupta, <a href="https://doi.org/10.1007/978-981-97-2465-9_16">Fascinating Factorials</a>, Exploring the Beauty of Fascinating Numbers, Springer (2025) Ch. 16, 411-442. %H A014597 Eric Weisstein's World of Mathematics, <a href="https://mathworld.wolfram.com/Factorial.html">Factorial</a> %e A014597 1183893^2 = 1! + 2! + 3! + 7! + 8! + 9! + 10! + 11! + 12! + 13! + 14! + 15!. %e A014597 2 is not a member since 4 is not a sum of distinct factorials. %t A014597 ok[n_] := (k=1; ff={}; While[k! < n^2, AppendTo[ff, k!]; k++]; xx = Array[x, Length[ff]]; Reduce[And @@ (0 <= # <= 1 & /@ xx) && n^2 == xx.ff, xx, Integers] =!= False); ok[1] = True; Reap[Do[If[ok[n], Print[n]; Sow[n]], {n, 1, 2*10^6}]][[2, 1]] (* _Jean-François Alcover_, Jul 16 2012 *) %o A014597 (Haskell) %o A014597 import Data.List (elemIndices) %o A014597 a014597 n = a014597_list !! (n-1) %o A014597 a014597_list = tail $ elemIndices 1 $ map a197183 [0..] %o A014597 -- _Reinhard Zumkeller_, Dec 04 2011 %o A014597 (Python) %o A014597 from math import factorial, isqrt %o A014597 from itertools import chain, combinations %o A014597 from sympy.ntheory.primetest import is_square %o A014597 fac =[factorial(n) for n in range(1, 16)] # raise 16 to search higher %o A014597 def powerset(s): # skipping empty set %o A014597 return chain.from_iterable(combinations(s, r) for r in range(1, len(s)+1)) %o A014597 gen = (isqrt(sum(s)) for s in powerset(fac) if is_square(sum(s))) %o A014597 print(sorted(set(gen))) # _Michael S. Branicky_, Jan 03 2021 %Y A014597 Cf. A025494, A051761, A059589. %K A014597 nonn,more,hard,nice %O A014597 1,2 %A A014597 _Eric W. Weisstein_ %E A014597 15th term from _Jud McCranie_, who remarks that there are no others involving terms < 21!.