%I #17 Apr 25 2023 03:37:23
%S 1,1,4,6,17,10,41,28,54,40,122,44,183,82,130,140,340,102,436,156,304,
%T 254,673,184,665,360,604,388,1143,232,1335,584,854,656,1138,484,1997,
%U 842,1238,768,2486,508,2762,1044,1382,1314,3339,816,3139,1160,2276,1600
%N a(n) = Sum_{1<=k<=n, gcd(k,n)=1} (Sum_{i=1..k} gcd(i, k)).
%C a(n) = Sum_{1<=k<=n, gcd(k,n)=1} A018804(k) where A018804(n) = Sum_{k=1..n} gcd(k,n).
%H Amiram Eldar, <a href="/A127416/b127416.txt">Table of n, a(n) for n = 1..10000</a>
%F M * V where M = A054521 is an infinite lower triangular matrix and V = A018804 is a vector: (1, 3, 5, 8, 9, 15, 13, ...).
%e a(6) = 10 since the relative primes of 6 are 1 and 5, A018804(1) + A018804(5) = 1 + 9 = 10.
%t f[p_, e_] := (e*(p - 1)/p + 1)*p^e; pil[n_] := Times @@ (f @@@ FactorInteger[n]); a[n_] := Sum[If[GCD[n, k] == 1, pil[k], 0], {k, 1, n}]; Array[a, 100] (* _Amiram Eldar_, Jul 19 2019 *)
%o (PARI) f(n) = sum(i=1, n, gcd(n, i)); \\ A018804
%o a(n) = sum(k=1, n, if (gcd(k, n) == 1, f(k))); \\ _Michel Marcus_, Jul 19 2019
%Y Cf. A018804, A054521.
%K nonn
%O 1,3
%A _Gary W. Adamson_, Jan 13 2007
%E More terms from _Amiram Eldar_, Jul 19 2019