|
| MATHEMATICA
|
max = 38; m[k_, n_] := (s = 1; Do[ If[ Mod[s, n] == 0, Return[e], s = s + k^e ] , {e, 1, max}]); c[k_, t_, n_] := Sum[ m[k, n/GCD[n, u*(k-1) + t]]^(-1), {u, 0, n-1}]; (* b = A002729 *) b[n_] := b[n] = (s = 0; Do[ If[ GCD[n, k] == 1 , s = s + 2^c[k, t, n]] , {k, 1, n-1}, {t, 0, n-1}]; s / (n*EulerPhi[n]) ); b[0] = 1; b[1] = 2; a[n_] := Sum[ MoebiusMu[n/d]*b[d], {d, Divisors[n]}]; Table[a[n], {n, 1, max}] (* From }] (* _Jean-François Alcover, _, Dec 06 2011, after Maple *)
|
