|
%I #38 Apr 05 2021 20:37:05
%S 1,7,751530,4730879,5260710,33116153,37200735,187062910,226141311,
%T 259109835,260405145,370049418,522409465,836308083,1105725765,
%U 1309440370,1343713507,1582989177,1609505430,1813768845,2590345926,3039492538,3656866255
%N Numbers n such that the sum of divisors of n^3 is a square.
%C In 1657 Fermat challenged the world to find such numbers. [Dickson, Vol. 1, p. 54]
%C If n is a term and n is not divisible by 7, then 7*n is a term. - Don Dechman (dondechman_2000(AT)yahoo.com), Mar 26 2008
%D A. H. Beiler, Recreations in the Theory of Numbers, Dover, NY, 1964, p. 9.
%D L. E. Dickson, History of the Theory of Numbers. Carnegie Institute Public. 256, Washington, DC, Vol. 1, 1919; Vol. 2, 1920; Vol. 3, 1923, see vol. 1, p. 54.
%D Julian Havil, Gamma, Exploring Euler's Constant, Princeton University Press, Princeton and Oxford, 2003, page 92.
%D I. Kaplansky, The challenges of Fermat, Wallis and Ozanam (and several related challenges): I. Fermat's first challenge, Preprint, 2002.
%H Donovan Johnson, <a href="/A008849/b008849.txt">Table of n, a(n) for n = 1..40</a> (terms < 10^11)
%H Donovan Johnson, <a href="http://donovanjohnson.com/OEIS/A008849/a008849.txt">73688 terms > 10^11</a>
%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/FermatsDivisorProblem.html">Fermat's Divisor Problem.</a>
%t max = 10^11; primes = {5, 7, 11, 13, 17, 19, 23, 31, 41, 43, 47, 83, 89, 191, 193, 239, 307, 443, 463, 499, 557, 701, 743, 1087, 1487, 2309, 3583, 4373, 5087, 5507, 5807, 44179}; subs = Select[ Times @@@ Subsets[primes, 7], # < max &] // Sort; f[e2_, e3_, p_] := If[n = 2^e2*3^e3*p; IntegerQ[ Sqrt[ DivisorSigma[1, n^3]]], Print[{2^e2, 3^e3, p}]; Sow[n]]; r = Reap[ Scan[ ((f[0, 0, #]; f[0, 1, #]; f[0, 3, #]; f[1, 0, #]; f[1, 1, #]; f[1, 3, #]; f[3, 0, #]; f[3, 1, #]; f[3, 3, #])& ), subs]][[2, 1]]; Select[r, # < max &] // Union (* _Jean-François Alcover_, Sep 07 2012, after _Donovan Johnson_ *)
%o (PARI) is(n)=issquare(sigma(n^3)) \\ _Charles R Greathouse IV_, Jun 20 2013
%o (Python)
%o from functools import reduce
%o from operator import mul
%o from sympy import factorint, integer_nthroot
%o A008849_list, n = [], 1
%o while n < 10**7:
%o fs = factorint(n)
%o if integer_nthroot(reduce(mul,((p**(3*fs[p]+1)-1)//(p-1) for p in fs),1),2)[1]:
%o A008849_list.append(n)
%o n += 1 # _Chai Wah Wu_, Apr 05 2021
%Y Cf. A046872, A008850, A048948.
%K nonn,nice
%O 1,2
%A _N. J. A. Sloane_
%E _David W. Wilson_ has supplied terms a(4) = 4730879 and beyond and verified completeness up to a(3) = 751530
%E I. Kaplansky and _Will Jagy_ have verified that there are no other terms below 3.8*10^9
%E 3656866255 added by Don Dechman, Mar 26 2008
|
