%I #20 Jul 04 2017 03:38:23
%S 1,657,3681,10809,15777,17937,24201,28521,54657,81441,122697,154881,
%T 230481,265257,336321,346041,455337,473481,547137,613017,718857,
%U 833121,898137,1161009,1226457,1274841,1305081,1510281,1584801,1618497,1695609,1752417,1846161,1965609
%N Numbers n such that the sum of the divisors of n is of the form m^2+1.
%C The corresponding values of m are 0, 31, 73, 125, 151, 161, 187, 203, 281, ...
%C For n > 1, it appears that a(n) is of the form 9p where p is prime congruent to 1 (mod 24). The corresponding primes p are 73, 409, 1201, 1753, 1993, 2689, 3169, 6073, 9049, 13633, 17209, ... (a subsequence of A107008?).
%C For n > 1, the divisors of a(n) are of the form {1, 3, 9, p, 3p, 9p} and sigma(a(n)) = 13(1 + p) = 1 + m^2.
%C Proof:
%C We use the formula (a^2 + b^2)(c^2 + d^2) = (ac - bd)^2 + (ad + bc)^2 with:
%C a = 2, b = 3, c = (3m + 2)/13 and d = (2m - 3)/13. We obtain:
%C a^2 + b^2 = 13, c^2 + d^2 = 1 + p = (m^2 + 1)/13 => p = (m^2 - 12)/13,
%C (ac - bd)^2 = 1 and (ad + bc)^2 = m^2.
%C The first terms not of the form 9p are 2493961, 7106353, 8325721, 10708297, and 14120281. Note that solving the Diophantine equation m^2 + 1 = sigma(k)*(1+p), for various numbers k, it is possible to identify other potentially infinite subsequences similar to 9p. For example, 47^2*p, 7^4*p, 3^6*p, 3^2*11^2*p, and so on. - _Giovanni Resta_, Jul 02 2017
%H Giovanni Resta, <a href="/A289290/b289290.txt">Table of n, a(n) for n = 1..10000</a>
%e 657 is in the sequence because sigma(657) = 962 = 31^2 + 1.
%p with(numtheory):nn:=10^6:
%p for n from 1 to nn do:
%p y:=sqrt(sigma(n)-1):
%p if y=floor(y) then printf(`%d, `,n):
%p else
%p fi:
%p od:
%t Select[Range[10^6], IntegerQ[Sqrt[DivisorSigma[1, #] - 1]] &] (* _Giovanni Resta_, Jul 02 2017 *)
%o (PARI) isok(n) = issquare(sigma(n) - 1); \\ _Michel Marcus_, Jul 02 2017
%Y Cf. A000203, A002522, A107008.
%K nonn
%O 1,2
%A _Michel Lagneau_, Jul 02 2017
%E Name edited by _Robert Israel_, Jul 03 2017