%I #33 Apr 28 2022 13:35:55
%S 1,2,136,156,328,444,584,600,712,732,776,876,904,1096,1164,1176,1308,
%T 1544,1864,1884,1928,2056,2172,2248,2316,2504,2601,2696,2748,2824,
%U 2892,2904,3208,3240,3249,3272,3324,3464,3592,3656,3756,4044,4056,4168,4188,4476,4552,4616
%N Numbers N such that gcd(N - d, N*d) >= d^2, where d = A000005(N) is the number of divisors of N.
%C As d^2 | N-d we have N = k*d^2 + d for some k >= 0 and d > 1. So gcd(k*d^2 + d - d, (N*d^2 + d)*d) = gcd(k*d^2, k*d^3 + d^2) = gcd(k*d^2, d^2) = d^2. So for any N such that d^2 | gcd(N - d, N*d) we have gcd(N - d, N*d) = d^2. - _David A. Corneth_, Apr 20 2022
%C Since gcd(N - d, N*d) is never larger than d^2 (if N = n*g, d = f*g with gcd(n,f) = 1, then gcd(N - d, N*d) = g*gcd(n-f,n*f*g) = g*gcd(n-f, f*f*g) <= g*g, since by assumption, no factor of f divides n), so one can also replace "=" by ">=" in the definition.
%H Michael De Vlieger, <a href="/A353012/b353012.txt">Table of n, a(n) for n = 1..10000</a>
%F For all m in A033950, the sequence contains all numbers m*p^k with k = m/d(m) - 1, and p^k == 1 (mod m), in particular 8*A007519 and 12*A068228 (k = 1, m = 8 and 12), 9*A129805^2, 18*A129805^2 and 24*A215848^2 (k = 2, m = 9, 18 and 24, A^2 = {x^2, x in A}), etc.
%e N = 1 is in the sequence because d(N) = 1, gcd(1 - 1, 1*1) = 1 = d^2.
%e N = 2 is in the sequence because d(N) = 2, gcd(2 - 2, 2*2) = 4 = d^2.
%e N = 136 = 8*17 is in the sequence because d(N) = 4*2 = 8, gcd(8*17 - 8, 8*17*8) = gcd(8*16, 8*8*17) = 8*8 = d^2. Similarly for N = 8*p with any prime p = 8*k + 1.
%e N = 156 = 2^2*3*13 is in the sequence because d(n) = 3*2*2 = 12, gcd(12*13 - 12, 12*13*12) = gcd(12*12, 12*12*13) = 12*12 = d^2. Similarly for any N = 12*p with prime p = 12*k + 1.
%e More generally, when N = m*p^k with p^k == 1 (mod m) and m = (k+1)*d(m), then d(N) = d(m)*(k+1) = m and gcd(n - d, n*d) = gcd(m*p^k - m, m*p^k*m) = m*gcd(p^k - 1, p^k*m) = m^2. This holds for m = 8 and 12 with k = 1, for m = 9, 18 and 24 with k = 2, etc: see sequence A033950 for the m-values.
%t Select[Range[4650], GCD[#1 - #2, #1 #2] == #2^2 & @@ {#, DivisorSigma[0, #]} &] (* _Michael De Vlieger_, Apr 21 2022 *)
%o (PARI) select( {is(n, d=numdiv(n))=gcd(n-d,d^2)==d^2}, [1..10^4])
%Y Cf. A000005 (number of divisors), A352483 (numerator of (n-d)/(n*d)), A352482 (denominator), A049820 (n - d), A146566 (n*d is divisible by n-d), A033950 (refactorable or tau numbers: d(n) | n, supersequence of this).
%K nonn
%O 1,2
%A _M. F. Hasler_, Apr 15 2022