|
%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
|
