A111901 - OEIS

%I #19 Apr 06 2017 09:46:49

%S 1,0,1,1,2,2,3,4,5,7,8,11,13,17,20,25,30,37,44,53,63,75,89,105,123,

%T 145,169,197,229,266,307,355,408,469,538,615,703,801,912,1035,1175,

%U 1330,1504,1698,1914,2155,2423,2721,3051,3418,3824,4273,4770,5319,5925

%N Number of partitions of n into parts that are primes or squares of primes.

%H Alois P. Heinz, <a href="/A111901/b111901.txt">Table of n, a(n) for n = 0..10000</a>

%F G.f.: Product_{k>=1} 1/((1 - x^prime(k))*(1 - x^(prime(k)^2))). - _Ilya Gutkovskiy_, Dec 26 2016

%e G.f. = 1 + x^2 + x^3 + 2*x^4 + 2*x^5 + 3*x^6 + 4*x^7 + 5*x^8 + 7*x^9 + 8*x^10 + ...

%e a(10) = #{ 7+3, 5+5, 5+3+2, 2^2+2^2+2, 2^2+3+3, 2^2+2+2+2, 3+3+2+2, 2+2+2+2+2 } = 8.

%p with(numtheory):

%p a:= proc(n) option remember; `if`(n=0, 1, add(a(n-j)*add(

%p      `if`(tau(d) in [2, 3], d, 0), d=divisors(j)), j=1..n)/n)

%p     end:

%p seq(a(n), n=0..60);  # _Alois P. Heinz_, Mar 30 2017

%t a[n_] := a[n] = If[n == 0, 1, Sum[a[n - j]*DivisorSum[j, If[2 <= DivisorSigma[0, #] <= 3, #, 0]&], {j, 1, n}]/n];

%t Table[a[n], {n, 0, 60}] (* _Jean-François Alcover_, Apr 06 2017, after _Alois P. Heinz_ *)

%o (PARI) {a(n) = if(n < 0, 0, polcoeff( 1 / prod(k=1, primepi(n), (1 - x^prime(k)^2 + x*O(x^n)) * (1 - x^prime(k))), n))}; /* _Michael Somos_, Dec 26 2016 */

%Y Cf. A000607, A090677, A023893, A111902.

%K nonn

%O 0,5

%A _Reinhard Zumkeller_, Aug 20 2005