%I #36 Jun 04 2016 06:59:36
%S 1,65,175,1105,5425,20737,32045,70525,103685,171275,200725,207553,
%T 352529,372775,1037765,1198925,1264957,1347905,1762645,1824877,
%U 2609425,2698189,3628975,3928475,4966975,6324785,6337175,8646625,8813225,9124385,10223341,12774139,13490945
%N Numbers that divide the average of the sum of the squares of their divisors.
%H Chai Wah Wu and Giovanni Resta, <a href="/A255244/b255244.txt">Table of n, a(n) for n = 1..328</a> (terms < 10^12, first 82 terms from Chai Wah Wu)
%e Divisors of 65 are 1, 5, 13, 65. The average of the sum of their squares is (1^2 + 5^2 + 13^2 + 65^2) / 4 = (1 + 25 + 169 + 4225) / 4 = 4420 / 4 = 1105 and 1105 / 65 = 17.
%p with(numtheory); P:=proc(q) local a,b,k,n;
%p for n from 2 to q do a:=divisors(n);
%p b:=add(a[k]^2,k=1..nops(a))/nops(a);
%p if type(b/n,integer) then lprint(n);
%p fi; od; end: P(10^6);
%t Select[Range[10^6],Mod[Mean[Divisors[#]^2],#]==0&] (* _Ivan N. Ianakiev_, Mar 03 2015 *)
%o (PARI) isok(n) = (q=sumdiv(n, d, d^2)/numdiv(n)) && (type(q)=="t_INT") && ((q % n) == 0); \\ _Michel Marcus_, Feb 20 2015
%o (Python)
%o from __future__ import division
%o from sympy import factorint
%o A255244_list = []
%o for n in range(1,10**9):
%o ....s0 = s2 = 1
%o ....for p,e in factorint(n).items():
%o ........s0 *= e+1
%o ........s2 *= (p**(2*(e+1))-1)//(p**2-1)
%o ....q, r = divmod(s2,s0)
%o ....if not (r or q % n):
%o ........A255244_list.append(n) # _Chai Wah Wu_, Mar 08 2015
%Y Cf. A000203, A255245.
%K nonn
%O 1,2
%A _Paolo P. Lava_, Feb 20 2015
%E More terms from _Michel Marcus_, Feb 20 2015
%E a(31)-a(33) corrected by _Chai Wah Wu_, Mar 08 2015