OFFSET
0,3
LINKS
Alois P. Heinz, Table of n, a(n) for n = 0..2000
Eric Weisstein's World of Mathematics, Unitary Divisor
FORMULA
a(n) = [x^n] 1/(1 - Sum_{d|n, gcd(d, n/d) = 1} x^d).
a(n) = 2 if n is a prime power (A246655).
EXAMPLE
a(8) = 2 because 8 has 4 divisors {1, 2, 4, 8} among which 2 are unitary divisors {1, 8} therefore we have [8] and [1, 1, 1, 1, 1, 1, 1, 1].
MAPLE
a:= proc(n) option remember; local b, l; l, b:=
select(x-> igcd(x, n/x)=1, numtheory[divisors](n)),
proc(m) option remember; `if`(m=0, 1,
add(`if`(j>m, 0, b(m-j)), j=l))
end; b(n)
end:
seq(a(n), n=0..60); # Alois P. Heinz, Aug 01 2017
MATHEMATICA
Join[{1}, Table[d = Divisors[n]; Coefficient[Series[1/(1 - Sum[Boole[GCD[n/d[[k]], d[[k]]] == 1] x^d[[k]], {k, Length[d]}]), {x, 0, n}], x, n], {n, 1, 53}]]
PROG
(Python)
from sympy import divisors, gcd
from sympy.core.cache import cacheit
@cacheit
def a(n):
l=[x for x in divisors(n) if gcd(x, n//x)==1]
@cacheit
def b(m): return 1 if m==0 else sum(b(m - j) for j in l if j <= m)
return b(n)
print([a(n) for n in range(61)]) # Indranil Ghosh, Aug 01 2017, after Maple code
CROSSREFS
KEYWORD
nonn
AUTHOR
Ilya Gutkovskiy, Aug 01 2017
STATUS
approved