A252424 Numbers k such that sum of odd divisors of k equals sum of squares of primes dividing k.

18, 36, 72, 144, 234, 288, 468, 576, 936, 1152, 1872, 2304, 3744, 4608, 7488, 9216, 14976, 18432, 29952, 36864, 59904, 73728, 119808, 147456, 239616, 294912, 479232, 589824, 958464, 1179648, 1916928, 2359296, 3833856, 4718592, 7667712, 9437184, 15335424, 18874368

Offset: 1

Author

Michel Lagneau, Dec 17 2014

Keywords

• nonn

Comments

Numbers k such that A000593(k) = A005063(k).

a(n) == 0 (mod 18), and the numbers 18*2^m, m = 0,1,... are in the sequence because the odd divisors are {1, 3, 9}, the prime factors are {2, 3} => 2^2 + 3^2 = 1 + 3 + 9 = 13.

The numbers of the form 18*13*2^m are in the sequence because the odd divisors are {1, 3, 9, 13, 39, 117}, the prime factors are {2, 3, 13} => 2^2 + 3^2 + 13^2 = 1 + 3 + 9 + 13 + 39 + 117 = 182.

Examples

			18 is in the sequence because the prime factors of 18 are {2, 3}, the odd divisors of 18 are {1, 3, 9} => 2^2 + 3^2 = 1 + 3 + 9 = 13.
Or 18 => A000593(18) = A005063(18) = 13.
		

Links

• Robert G. Wilson v, Table of n, a(n) for n = 1..56

Crossrefs

Cf. A000593, A005063.

Programs

Maple

with(numtheory):nn:=10^5:
for n from 2 to nn do:
   x:=factorset(n):n0:=nops(x):
   s0:=sum('x[i]^2','i'=1..n0):
   y:=divisors(n):n1:=nops(y):
   s :=0 :
        for j from 1 to n1 do :
       if irem (y[j],2)=1 then s:=s+y[j]:
      else
      fi:
    od:
     if s=s0
    then
   printf(`%d, `,n):
   else
   fi:
od:

Mathematica

a252424[n_Integer] := Module[{f, g},
  f[x_] := Plus @@ Select[Divisors[x], OddQ[#] &];
  g[x_] := Plus @@ (First@Transpose@FactorInteger[x]^2);
Rest@Select[Range[n], f[#] == g[#] &]]; a252424[10^6] (* Michael De Vlieger, Dec 17 2014 *)
Select[Range[19*10^6],Total[Select[Divisors[#],OddQ]]==Total[ FactorInteger[ #][[All,1]]^2]&] (* Harvey P. Dale, May 11 2020 *)
f[p_, e_] := If[p == 2, 1, (p^(e + 1) - 1)/(p - 1)]; q[n_] := Times @@ f @@@ (fct = FactorInteger[n]) == Total[fct[[;; , 1]]^2]; Select[Range[2, 10^6], q] (* Amiram Eldar, Jul 09 2022 *)

PARI

isok(n) = my(f = factor(n)); sum(i=1, #f~, f[i,1]^2) == sumdiv(n, d, d*(d%2)); \\ Michel Marcus, Dec 17 2014