A333613 a(1) = 1; thereafter a(n) = Sum_{k = 1..n} a(k/gcd(n,k)).

1, 2, 4, 7, 15, 21, 51, 78, 158, 230, 568, 661, 1797, 2595, 5117, 7789, 19095, 21702, 59892, 81801, 171329, 258028, 630942, 713093, 1887828, 2776798, 5727675, 8335692, 20702970, 21420664, 62826604, 92041835, 189376593, 281410640, 656577018, 742729123, 2087788417

Offset

1,2

Links

Vaclav Kotesovec, Table of n, a(n) for n = 1..3930

Vaclav Kotesovec, Plot of a(n)^(1/n) for n = 1..10000

Formula

a(1) = 1; a(n) = Sum_{k = 1..n} a(lcm(n, k)/n).

a(1) = 1; a(n) = Sum_{d|n} Sum_{k = 1..d, gcd(d, k) = 1} a(k).

Maple

A333613:= proc(n) option remember;
if n<3 then n;
else add( A333613(lcm(n, j)/n), j = 1..n);
end if; end proc;
seq(A333613(n), n=1..40); # G. C. Greubel, Mar 08 2021

Mathematica

a[1] = 1; a[n_] := a[n] = Sum[a[k/GCD[n, k]], {k, n}]; Table[a[n], {n, 37}]
a[1] = 1; a[n_] := a[n] = Sum[Sum[If[GCD[k, d] == 1, a[k], 0], {k, d}], {d, Divisors[n]}]; Table[a[n], {n, 37}]

Prog

(Sage)
@CachedFunction
def A333613(n): return 1 if n==1 else sum( A333613(lcm(n, j)/n) for j in (1..n) )
[A333613(n) for n in (1..40)] # G. C. Greubel, Mar 08 2021

Crossrefs

Cf. A006874, A045545, A057661, A332791.

Sequence in context: A152477 A019998 A002956 * A171276 A027167 A259090

Adjacent sequences: A333610 A333611 A333612 * A333614 A333615 A333616

Keyword

nonn

Author

Ilya Gutkovskiy, Mar 28 2020

Status

approved