login
Search: author:ilya gutkovskiy
     Sort: relevance | references | number | modified | created      Format: long | short | data
Powerful numbers (A001694) that are not prime powers (A000961).
+10
74
36, 72, 100, 108, 144, 196, 200, 216, 225, 288, 324, 392, 400, 432, 441, 484, 500, 576, 648, 675, 676, 784, 800, 864, 900, 968, 972, 1000, 1089, 1125, 1152, 1156, 1225, 1296, 1323, 1352, 1372, 1444, 1521, 1568, 1600, 1728, 1764, 1800, 1936, 1944, 2000, 2025, 2116, 2304, 2312, 2500, 2592, 2601, 2700, 2704, 2744
OFFSET
1,1
COMMENTS
If a prime p divides a(n) then p^2 must also divide a(n) and number of distinct primes dividing a(n) > 1.
Intersection of A001694 and A024619.
LINKS
Michael De Vlieger, Table of n, a(n) for n = 1..10000 (first 5997 terms from Robert Israel)
Eric Weisstein's World of Mathematics, Prime Power.
Eric Weisstein's World of Mathematics, Powerful Number.
FORMULA
Sum_{n>=1} 1/a(n) = zeta(2)*zeta(3)/zeta(6) - Sum_{p prime} 1/(p*(p-1)) - 1 = A082695 - A136141 - 1 = 0.17043976777096407719... - Amiram Eldar, Feb 12 2021
EXAMPLE
-------------------------------
| n | a(n) | prime |
| | | factorization |
|------------------------------
| 1 | 36 | {{2, 2}, {3, 2}} |
| 2 | 72 | {{2, 3}, {3, 2}} |
| 3 | 100 | {{2, 2}, {5, 2}} |
| 4 | 108 | {{2, 2}, {3, 3}} |
| 5 | 144 | {{2, 4}, {3, 2}} |
| 6 | 196 | {{2, 2}, {7, 2}} |
| 7 | 200 | {{2, 3}, {5, 2}} |
| 8 | 216 | {{2, 3}, {3, 3}} |
| 9 | 225 | {{3, 2}, {5, 2}} |
-------------------------------
a(n) = p_1^e_1*p_2^e_2*... : {{p_1, e_1}, {p_2, e_2}, ...}.
MAPLE
N:= 10000:
S:= {1}: P:= {1}:
p:= 1:
do
p:= nextprime(p);
if p^2 > N then break fi;
S:= map(s -> (s, seq(s*p^k, k = 2 .. floor(log[p](N/s)))), S);
P:= P union {seq(p^k, k=2..floor(log[p](N)))}:
od:
sort(convert(S minus P, list)); # Robert Israel, May 14 2017
MATHEMATICA
Select[Range@2750, Min@FactorInteger[#][[All, 2]] > 1 && ! PrimePowerQ[#] &]
(* Second program *)
nn = 2^25; Select[Rest@ Union@ Flatten@ Table[a^2*b^3, {b, nn^(1/3)}, {a, Sqrt[nn/b^3]}], ! PrimePowerQ[#] &] (* Michael De Vlieger, Jun 22 2022 *)
PROG
(Python)
from sympy import primefactors, factorint
print([n for n in range(4, 2745) if len(primefactors(n)) > 1 and min(list(factorint(n).values())) > 1]) # Karl-Heinz Hofmann, Feb 07 2023
(Python)
from math import isqrt
from sympy import integer_nthroot, primepi, mobius
def A286708(n):
def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def f(x):
c, l = n+x, 0
j = isqrt(x)
while j>1:
k2 = integer_nthroot(x//j**2, 3)[0]+1
w = squarefreepi(k2-1)
c -= j*(w-l)
l, j = w, isqrt(x//k2**3)
c -= squarefreepi(integer_nthroot(x, 3)[0])-l
return c+1+sum(primepi(integer_nthroot(x, k)[0]) for k in range(2, x.bit_length()))
return bisection(f, n, n) # Chai Wah Wu, Sep 10 2024
KEYWORD
nonn
AUTHOR
_Ilya Gutkovskiy_, May 13 2017
STATUS
approved
Irregular triangle T(n,k) read by rows: first row is 0, n-th row (n > 1) lists indices of distinct primes dividing n.
+10
52
0, 1, 2, 1, 3, 1, 2, 4, 1, 2, 1, 3, 5, 1, 2, 6, 1, 4, 2, 3, 1, 7, 1, 2, 8, 1, 3, 2, 4, 1, 5, 9, 1, 2, 3, 1, 6, 2, 1, 4, 10, 1, 2, 3, 11, 1, 2, 5, 1, 7, 3, 4, 1, 2, 12, 1, 8, 2, 6, 1, 3, 13, 1, 2, 4, 14, 1, 5, 2, 3, 1, 9, 15, 1, 2, 4, 1, 3, 2, 7, 1, 6, 16, 1, 2, 3, 5, 1, 4, 2, 8, 1, 10, 17, 1, 2, 3, 18, 1, 11
OFFSET
1,3
FORMULA
T(n,k) = A000720(A027748(n,k)).
T(n,1) = A055396(n).
T(n,A001221(n)) = A061395(n).
EXAMPLE
The irregular triangle begins:
1: {0}
2: {1}
3: {2}
4: {1}
5: {3}
6: {1, 2}
7: {4}
8: {1}
9: {2}
10: {1, 3}
11: {5}
12: {1, 2}
MATHEMATICA
Flatten[Table[PrimePi[FactorInteger[n][[All, 1]]], {n, 1, 62}]]
CROSSREFS
Cf. A000040, A000720, A001221 (row lengths), A027748, A055396, A061395, A066328 (row sums), A112798, A156061 (row products), A302170.
KEYWORD
nonn,tabf
AUTHOR
_Ilya Gutkovskiy_, May 05 2018
STATUS
approved
Expansion of 1/(1 - Sum_{k>=1} q(k)*x^k), where q(k) = number of partitions of k into distinct parts (A000009).
+10
46
1, 1, 2, 5, 11, 25, 57, 129, 292, 662, 1500, 3398, 7699, 17443, 39519, 89536, 202855, 459593, 1041267, 2359122, 5344889, 12109524, 27435660, 62158961, 140828999, 319065932, 722884274, 1637785870, 3710611298, 8406859805, 19046805534, 43152950024, 97768473163
OFFSET
0,3
COMMENTS
Invert transform of A000009.
From Gus Wiseman, Jul 31 2022: (Start)
Also the number of ways to choose a multiset partition into distinct constant multisets of a multiset of length n that covers an initial interval of positive integers. This interpretation involves only multisets, not sequences. For example, the a(1) = 1 through a(4) = 11 multiset partitions are:
{{1}} {{1,1}} {{1,1,1}} {{1,1,1,1}}
{{1},{2}} {{1},{1,1}} {{1},{1,1,1}}
{{1},{2,2}} {{1,1},{2,2}}
{{2},{1,1}} {{1},{2,2,2}}
{{1},{2},{3}} {{2},{1,1,1}}
{{1},{2},{1,1}}
{{1},{2},{2,2}}
{{1},{2},{3,3}}
{{1},{3},{2,2}}
{{2},{3},{1,1}}
{{1},{2},{3},{4}}
The non-strict version is A055887.
The strongly normal non-strict version is A063834.
The strongly normal version is A270995.
(End)
FORMULA
G.f.: 1/(1 - Sum_{k>=1} A000009(k)*x^k).
G.f.: 1/(2 - Product_{k>=1} (1 + x^k)).
G.f.: 1/(2 - Product_{k>=1} 1/(1 - x^(2*k-1))).
G.f.: 1/(2 - exp(Sum_{k>=1} (-1)^(k+1)*x^k/(k*(1 - x^k)))).
a(n) ~ c / r^n, where r = 0.441378990861652015438479635503868737167721352874... is the root of the equation QPochhammer[-1, r] = 4 and c = 0.4208931614610039677452560636348863586180784719323982664940444607322... - Vaclav Kotesovec, May 23 2018
EXAMPLE
From Gus Wiseman, Jul 31 2022: (Start)
a(n) is the number of ways to choose a strict integer partition of each part of an integer composition of n. The a(1) = 1 through a(4) = 11 choices are:
((1)) ((2)) ((3)) ((4))
((1)(1)) ((21)) ((31))
((1)(2)) ((1)(3))
((2)(1)) ((2)(2))
((1)(1)(1)) ((3)(1))
((1)(21))
((21)(1))
((1)(1)(2))
((1)(2)(1))
((2)(1)(1))
((1)(1)(1)(1))
(End)
MAPLE
b:= proc(n) option remember; `if`(n=0, 1, add(b(n-j)*add(
`if`(d::odd, d, 0), d=numtheory[divisors](j)), j=1..n)/n)
end:
a:= proc(n) option remember; `if`(n=0, 1,
add(b(j)*a(n-j), j=1..n))
end:
seq(a(n), n=0..40); # Alois P. Heinz, May 22 2018
MATHEMATICA
nmax = 32; CoefficientList[Series[1/(1 - Sum[PartitionsQ[k] x^k, {k, 1, nmax}]), {x, 0, nmax}], x]
nmax = 32; CoefficientList[Series[1/(2 - Product[1 + x^k, {k, 1, nmax}]), {x, 0, nmax}], x]
nmax = 32; CoefficientList[Series[1/(2 - 1/QPochhammer[x, x^2]), {x, 0, nmax}], x]
a[0] = 1; a[n_] := a[n] = Sum[PartitionsQ[k] a[n - k], {k, 1, n}]; Table[a[n], {n, 0, 32}]
CROSSREFS
Row sums of A308680.
The unordered version is A089259, non-strict A001970 (row-sums of A061260).
For partitions instead of compositions we have A270995, non-strict A063834.
A000041 counts integer partitions, strict A000009.
A072233 counts partitions by sum and length.
Cf. A279784.
KEYWORD
nonn
AUTHOR
_Ilya Gutkovskiy_, May 22 2018
STATUS
approved
Number of distinct prime divisors of n that are < sqrt(n).
+10
38
0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 2, 0, 1, 1, 1, 0, 2, 0, 1, 1, 1, 0, 2, 0, 1, 1, 1, 0, 3, 0, 1, 1, 1, 1, 2, 0, 1, 1, 2, 0, 2, 0, 1, 2, 1, 0, 2, 0, 2, 1, 1, 0, 2, 1, 2, 1, 1, 0, 3, 0, 1, 2, 1, 1, 2, 0, 1, 1, 3, 0, 2, 0, 1, 2, 1, 1, 2, 0, 2, 1, 1, 0, 3, 1, 1, 1, 1, 0, 3
OFFSET
1,12
COMMENTS
a(n) = 0 if and only if n = p^k where p is prime and k is 0, 1, or 2. - Charles R Greathouse IV, Apr 07 2020
LINKS
FORMULA
G.f.: Sum_{k>=1} x^(prime(k)*(prime(k) + 1)) / (1 - x^prime(k)).
a(k*n) >= a(n) for k > 0. a(n^e) = A001221(n) for e > 2. - Richard Peterson, Dec 19 2024
MAPLE
N:= 100: # for a(1)..a(N)
V:= Vector(N):
p:= 1:
do
p:= nextprime(p);
if p^2 >= N then break fi;
L:= [seq(p*k, k=p+1..N/p)];
V[L]:= V[L]+~1
od:
convert(V, list); # Robert Israel, Apr 07 2020
MATHEMATICA
Table[DivisorSum[n, 1 &, # < Sqrt[n] && PrimeQ[#] &], {n, 1, 90}]
nmax = 90; CoefficientList[Series[Sum[x^(Prime[k] (Prime[k] + 1))/(1 - x^Prime[k]), {k, 1, nmax}], {x, 0, nmax}], x] // Rest
PROG
(PARI) a(n)=my(f=factor(n)[, 1]); sum(i=1, #f, f[i]^2<n) \\ Charles R Greathouse IV, Apr 07 2020
CROSSREFS
KEYWORD
nonn
AUTHOR
_Ilya Gutkovskiy_, Apr 05 2020
STATUS
approved
Expansion of Product_{k>=1} (1 + 2^(k-1)*x^k).
+10
37
1, 1, 2, 6, 12, 32, 72, 176, 384, 960, 2112, 4992, 11264, 26112, 58368, 136192, 301056, 688128, 1548288, 3489792, 7766016, 17596416, 38993920, 87293952, 194248704, 432537600, 957349888, 2132803584, 4699717632, 10406068224, 23001563136, 50683969536, 111434268672, 245819768832
OFFSET
0,3
COMMENTS
Number of compositions of partitions of n into distinct parts. a(3) = 6: 3, 21, 12, 111, 2|1, 11|1. - Alois P. Heinz, Sep 16 2019
Also the number of ways to split a composition of n into contiguous subsequences with strictly decreasing sums. - Gus Wiseman, Jul 13 2020
This sequence is obtained from the generalized Euler transform in A266964 by taking f(n) = -1, g(n) = (-1) * 2^(n-1). - Seiichi Manyama, Aug 22 2020
FORMULA
G.f.: Product_{k>=1} (1 + A011782(k)*x^k).
a(n) ~ 2^n * exp(2*sqrt(-polylog(2, -1/2)*n)) * (-polylog(2, -1/2))^(1/4) / (sqrt(6*Pi) * n^(3/4)). - Vaclav Kotesovec, Sep 19 2019
EXAMPLE
From Gus Wiseman, Jul 13 2020: (Start)
The a(0) = 1 through a(4) = 12 splittings:
() (1) (2) (3) (4)
(1,1) (1,2) (1,3)
(2,1) (2,2)
(1,1,1) (3,1)
(2),(1) (1,1,2)
(1,1),(1) (1,2,1)
(2,1,1)
(3),(1)
(1,1,1,1)
(1,2),(1)
(2,1),(1)
(1,1,1),(1)
(End)
MATHEMATICA
nmax = 33; CoefficientList[Series[Product[(1 + 2^(k - 1) x^k), {k, 1, nmax}], {x, 0, nmax}], x]
PROG
(PARI) N=40; x='x+O('x^N); Vec(prod(k=1, N, 1+2^(k-1)*x^k)) \\ Seiichi Manyama, Aug 22 2020
CROSSREFS
The non-strict version is A075900.
Starting with a reversed partition gives A323583.
Starting with a partition gives A336134.
Partitions of partitions are A001970.
Splittings with equal sums are A074854.
Splittings of compositions are A133494.
Splittings with distinct sums are A336127.
KEYWORD
nonn
AUTHOR
_Ilya Gutkovskiy_, May 22 2018
STATUS
approved
Number of odd divisors of n that are < sqrt(n).
+10
36
0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 3, 1, 1, 2, 1, 2, 2, 1, 1, 2, 2, 1, 2, 1, 1, 3, 1, 1, 2, 1, 2, 2, 1, 1, 2, 2, 2, 2, 1, 1, 3, 1, 1, 3, 1, 2, 2, 1, 1, 2, 3, 1, 2, 1, 1, 3, 1, 2, 2, 1, 2, 2, 1, 1, 3, 2, 1, 2, 1, 1, 4, 2, 1, 2, 1, 2, 2, 1, 2, 3, 2, 1, 2, 1, 1, 4
OFFSET
1,12
COMMENTS
If we define a divisor d|n to be strictly inferior if d < n/d, then strictly inferior divisors are counted by A056924 and listed by A341674. This sequence counts strictly inferior odd divisors. - Gus Wiseman, Feb 26 2021
LINKS
FORMULA
G.f.: Sum_{k>=1} x^(2*k*(2*k - 1)) / (1 - x^(2*k - 1)).
EXAMPLE
The strictly inferior odd divisors of 945 are 1, 3, 5, 7, 9, 15, 21, 27, so a(945) = 8. - Gus Wiseman, Feb 27 2021
MATHEMATICA
Table[DivisorSum[n, 1 &, # < Sqrt[n] && OddQ[#] &], {n, 1, 90}]
nmax = 90; CoefficientList[Series[Sum[x^(2 k (2 k - 1))/(1 - x^(2 k - 1)), {k, 1, nmax}], {x, 0, nmax}], x] // Rest
PROG
(PARI) A333805(n) = sumdiv(n, d, (d%2)&&((d*d)<n)); \\ Antti Karttunen, Nov 02 2022
CROSSREFS
Dominated by A001227 (number of odd divisors).
Strictly inferior divisors (not just odd) are counted by A056924.
The non-strict version is A069288.
These divisors add up to A070039.
The case of prime divisors is A333806.
The strictly superior version is A341594.
The case of squarefree divisors is A341596.
The superior version is A341675.
The case of prime-power divisors is A341677.
A006530 selects the greatest prime factor.
A020639 selects the smallest prime factor.
- Odd -
A000009 counts partitions into odd parts, ranked by A066208.
A026424 lists numbers with odd Omega.
A027193 counts odd-length partitions.
A067659 counts strict partitions of odd length, ranked by A030059.
- Inferior divisors -
A033676 selects the greatest inferior divisor.
A033677 selects the smallest superior divisor.
A038548 counts superior (or inferior) divisors.
A060775 selects the greatest strictly inferior divisor.
A341674 lists strictly inferior divisors.
KEYWORD
nonn
AUTHOR
_Ilya Gutkovskiy_, Apr 05 2020
EXTENSIONS
Data section extended up to a(105) by Antti Karttunen, Nov 02 2022
STATUS
approved
Square array A(n,k), n>=0, k>=0, read by antidiagonals, where column k is the expansion of Product_{j>=1} (1 + x^j)^k.
+10
35
1, 1, 0, 1, 1, 0, 1, 2, 1, 0, 1, 3, 3, 2, 0, 1, 4, 6, 6, 2, 0, 1, 5, 10, 13, 9, 3, 0, 1, 6, 15, 24, 24, 14, 4, 0, 1, 7, 21, 40, 51, 42, 22, 5, 0, 1, 8, 28, 62, 95, 100, 73, 32, 6, 0, 1, 9, 36, 91, 162, 206, 190, 120, 46, 8, 0, 1, 10, 45, 128, 259, 384, 425, 344, 192, 66, 10, 0
OFFSET
0,8
COMMENTS
A(n,k) is the number of partitions of n into distinct parts (or odd parts) with k types of each part.
LINKS
N. J. A. Sloane, Transforms
FORMULA
G.f. of column k: Product_{j>=1} (1 + x^j)^k.
A(n,k) = Sum_{i=0..k} binomial(k,i) * A308680(n,k-i). - Alois P. Heinz, Aug 29 2019
EXAMPLE
A(3,2) = 6 because we have [3], [3'], [2, 1], [2', 1], [2, 1'] and [2', 1'] (partitions of 3 into distinct parts with 2 types of each part).
Also A(3,2) = 6 because we have [3], [3'], [1, 1, 1], [1, 1, 1'], [1, 1', 1'] and [1', 1', 1'] (partitions of 3 into odd parts with 2 types of each part).
Square array begins:
1, 1, 1, 1, 1, 1, ...
0, 1, 2, 3, 4, 5, ...
0, 1, 3, 6, 10, 15, ...
0, 2, 6, 13, 24, 40, ...
0, 2, 9, 24, 51, 95, ...
0, 3, 14, 42, 100, 206, ...
MAPLE
b:= proc(n, i, k) option remember; `if`(n=0, 1, `if`(i<1, 0, add(
(t-> b(t, min(t, i-1), k)*binomial(k, j))(n-i*j), j=0..n/i)))
end:
A:= (n, k)-> b(n$2, k):
seq(seq(A(n, d-n), n=0..d), d=0..12); # Alois P. Heinz, Aug 29 2019
MATHEMATICA
Table[Function[k, SeriesCoefficient[Product[(1 + x^i)^k , {i, Infinity}], {x, 0, n}]][j - n], {j, 0, 11}, {n, 0, j}] // Flatten
CROSSREFS
Columns k=0-32 give: A000007, A000009, A022567-A022596.
Rows n=0-2 give: A000012, A001477, A000217.
Main diagonal gives A270913.
Antidiagonal sums give A299106.
KEYWORD
nonn,tabl
AUTHOR
_Ilya Gutkovskiy_, May 07 2017
STATUS
approved
Expansion of Sum_{k>=2} x^(k*(k+3)/2) / Product_{j=1..k} (1 - x^j).
+10
34
0, 0, 0, 0, 0, 1, 1, 2, 2, 4, 4, 6, 7, 9, 11, 14, 16, 20, 24, 28, 34, 40, 47, 55, 65, 75, 88, 102, 118, 136, 158, 180, 208, 238, 272, 311, 355, 403, 459, 521, 590, 668, 756, 852, 962, 1084, 1218, 1370, 1538, 1724, 1932, 2163, 2417, 2701, 3015, 3361, 3745, 4170, 4636, 5154, 5724
OFFSET
0,8
COMMENTS
Number of partitions of n into at least two distinct parts >= 2.
FORMULA
G.f.: x - 1/(1 - x) + Product_{k>=2} (1 + x^k).
a(n) = A025147(n) - 1 for n > 1.
EXAMPLE
a(9) = 4 because we have [7, 2], [6, 3], [5, 4] and [4, 3, 2].
MATHEMATICA
nmax = 60; CoefficientList[Series[Sum[x^(k (k + 3)/2)/Product[(1 - x^j), {j, 1, k}], {k, 2, nmax}], {x, 0, nmax}], x]
nmax = 60; CoefficientList[Series[x - 1/(1 - x) + 1/((1 + x) QPochhammer[x, x^2]), {x, 0, nmax}], x]
Join[{0, 0}, Table[-1 + Sum[(-1)^(n - k) PartitionsQ[k], {k, 0, n}], {n, 2, 60}]]
CROSSREFS
KEYWORD
nonn
AUTHOR
_Ilya Gutkovskiy_, Aug 13 2018
STATUS
approved
Square array A(n,k), n >= 0, k >= 0, read by antidiagonals, where column k is the expansion of Product_{j>=1} 1/(1 - j*x^j)^k.
+10
33
1, 1, 0, 1, 1, 0, 1, 2, 3, 0, 1, 3, 7, 6, 0, 1, 4, 12, 18, 14, 0, 1, 5, 18, 37, 49, 25, 0, 1, 6, 25, 64, 114, 114, 56, 0, 1, 7, 33, 100, 219, 312, 282, 97, 0, 1, 8, 42, 146, 375, 676, 855, 624, 198, 0, 1, 9, 52, 203, 594, 1276, 2030, 2178, 1422, 354, 0, 1, 10, 63, 272, 889, 2196, 4155, 5736, 5496, 3058, 672, 0
OFFSET
0,8
LINKS
FORMULA
G.f. of column k: Product_{j>=1} 1/(1 - j*x^j)^k.
A(0,k) = 1; A(n,k) = (k/n) * Sum_{j=1..n} A078308(j) * A(n-j,k). - Seiichi Manyama, Aug 16 2023
EXAMPLE
G.f. of column k: A_k(x) = 1 + k*x + (1/2)*k*(k + 5)*x^2 + (1/6)*k*(k^2 + 15*k + 20)*x^3 + (1/24)*k*(k^3 + 30*k^2 + 155*k + 150)*x^4 + (1/120)*k*(k^4 + 50*k^3 + 575*k^2 + 1750*k + 624)*x^5 + ...
Square array begins:
1, 1, 1, 1, 1, 1, ...
0, 1, 2, 3, 4, 5, ...
0, 3, 7, 12, 18, 25, ...
0, 6, 18, 37, 64, 100, ...
0, 14, 49, 114, 219, 375, ...
0, 25, 114, 312, 676, 1276, ...
MATHEMATICA
Table[Function[k, SeriesCoefficient[Product[1/(1 - i x^i)^k, {i, 1, n}], {x, 0, n}]][j - n], {j, 0, 11}, {n, 0, j}] // Flatten
PROG
(PARI) first(n, k) = my(res = matrix(n, k)); for(u=1, k, my(col = Vec(prod(j=1, n, 1/(1 - j*x^j)^(u-1)) + O(x^n))); for(v=1, n, res[v, u] = col[v])); res \\ Iain Fox, Dec 28 2017
KEYWORD
nonn,tabl
AUTHOR
_Ilya Gutkovskiy_, Dec 28 2017
STATUS
approved
Expansion of (-1 + Product_{k>=1} 1 / (1 + (-x)^k))^4.
+10
32
1, 0, 4, 4, 10, 16, 26, 44, 63, 100, 144, 212, 297, 420, 584, 796, 1081, 1452, 1940, 2556, 3355, 4372, 5668, 7288, 9327, 11892, 15076, 19012, 23884, 29904, 37276, 46284, 57276, 70680, 86918, 106528, 130220, 158784, 193054, 234076, 283178, 341824, 411616, 494512, 592933
OFFSET
4,3
LINKS
FORMULA
G.f.: (-1 + Product_{k>=1} (1 + x^(2*k - 1)))^4.
a(n) ~ A112160(n). - Vaclav Kotesovec, Feb 20 2021
MAPLE
g:= proc(n) option remember; `if`(n=0, 1, add(add([0, d, -d, d]
[1+irem(d, 4)], d=numtheory[divisors](j))*g(n-j), j=1..n)/n)
end:
b:= proc(n, k) option remember; `if`(k<2, `if`(n=0, 1-k, g(n)),
(q-> add(b(j, q)*b(n-j, k-q), j=0..n))(iquo(k, 2)))
end:
a:= n-> b(n, 4):
seq(a(n), n=4..48); # Alois P. Heinz, Feb 07 2021
MATHEMATICA
nmax = 48; CoefficientList[Series[(-1 + Product[1/(1 + (-x)^k), {k, 1, nmax}])^4, {x, 0, nmax}], x] // Drop[#, 4] &
KEYWORD
nonn
AUTHOR
_Ilya Gutkovskiy_, Feb 07 2021
STATUS
approved

Search completed in 1.118 seconds