OFFSET
1,2
COMMENTS
The corresponding values of m are 0, 31, 73, 125, 151, 161, 187, 203, 281, ...
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?).
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.
Proof:
We use the formula (a^2 + b^2)(c^2 + d^2) = (ac - bd)^2 + (ad + bc)^2 with:
a = 2, b = 3, c = (3m + 2)/13 and d = (2m - 3)/13. We obtain:
a^2 + b^2 = 13, c^2 + d^2 = 1 + p = (m^2 + 1)/13 => p = (m^2 - 12)/13,
(ac - bd)^2 = 1 and (ad + bc)^2 = m^2.
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
LINKS
Giovanni Resta, Table of n, a(n) for n = 1..10000
EXAMPLE
657 is in the sequence because sigma(657) = 962 = 31^2 + 1.
MAPLE
with(numtheory):nn:=10^6:
for n from 1 to nn do:
y:=sqrt(sigma(n)-1):
if y=floor(y) then printf(`%d, `, n):
else
fi:
od:
MATHEMATICA
Select[Range[10^6], IntegerQ[Sqrt[DivisorSigma[1, #] - 1]] &] (* Giovanni Resta, Jul 02 2017 *)
PROG
(PARI) isok(n) = issquare(sigma(n) - 1); \\ Michel Marcus, Jul 02 2017
CROSSREFS
KEYWORD
nonn
AUTHOR
Michel Lagneau, Jul 02 2017
EXTENSIONS
Name edited by Robert Israel, Jul 03 2017
STATUS
approved