%I #33 Feb 15 2022 02:46:39
%S 1,1,4,7,20,16,56,50,93,80,220,110,364,224,340,372,816,354,1140,580,
%T 966,880,2024,820,2200,1456,2304,1666,4060,1240,4960,2856,3850,3264,
%U 5180,2706,8436,4560,6396,4440,11480,3612,13244,6710,8400,8096,17296,6344,17297,8600
%N a(n) = Sum_{1<=k<=n, gcd(k,n)=1}, A000217(k).
%C From _Wolfdieter Lang_, Jun 14 2011: (Start)
%C Such sums are over a reduced residue system modulo n. See the Apostol reference, p. 133, for the definition or the wikipedia link given under A189918.
%C This sum over triangular numbers can be found using the results given in exercise 16 of the Apostol reference on p. 48, together with the definition of phi_1(n) and phi_2(n) from the exercise 15.
%C The result for n >= 2 coincides with the formula given below, using Product_{p|n} (1 - p) = mu(rad(n))*rad(n)*phi(n)/n, with the definitions given there.
%C (End)
%D T. Apostol, Introduction to Analytic Number Theory, Springer, 1986.
%H Seiichi Manyama, <a href="/A127415/b127415.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 = A000217: (1, 3, 6, 10, ...).
%F From _Wolfdieter Lang_, May 17 2011: (Start)
%F a(n) = (n/(3!*2))*((2*n+3)*n + mu(rad(n))*rad(n))*(phi(n)/n), n >= 2, with rad(n) the squarefree kernel of n (the largest squarefree number dividing n, see A007947), the Moebius function mu(n)=A008683(n), and the Euler totient function phi(n)= A000010(n).
%F Note that phi(n)/n = A076512(n)/A109395(n) = phi(rad(n))/rad(n).
%F Proof via inclusion-exclusion.
%F (End)
%e a(6) = 16 since the relative primes of 6 are 1 and 5 and (1 + 15) = 16.
%e a(6) = (6/(3!*2))*(15*6 + 1*6)*(1/2)*(2/3)= 16.
%t rad[n_] := Times @@ (FactorInteger[n][[ All, 1]]); a[n_] := (n/(3!*2))*((2*n+3)*n + MoebiusMu[ rad[n]]*rad[n])*(EulerPhi[n] / n); a[1] = 1; Table[ a[n], {n, 1, 33}] (* _Jean-François Alcover_, Oct 03 2011 *)
%o (PARI) a(n)=if(n<3,return(1));my(s=factor(n)[,1]); s=prod(i=1,#s,s[i]); (n/12)*((2*n+3)*n + moebius(s)*s)*(eulerphi(n)/n) \\ _Charles R Greathouse IV_, May 17 2011
%o (PARI) a(n) = sum(k=1, n, if (gcd(n,k)==1, k*(k+1)/2)); \\ _Michel Marcus_, Feb 01 2016
%Y Cf. A000217, A023896, A076512/A109395, A189918.
%K nonn,easy
%O 1,3
%A _Gary W. Adamson_, Jan 13 2007
%E More terms and formula from _Wolfdieter Lang_, May 17 2011
%E More terms from _Michel Marcus_, Feb 01 2016