%I #32 Aug 05 2024 06:16:14
%S 1,5,12,24,37,61,80,112,145,189,220,288,325,389,464,544,593,701,756,
%T 880,989,1093,1160,1336,1441,1565,1700,1880,1965,2205,2296,2488,2665,
%U 2829,3028,3328,3437,3621,3832,4152,4273,4621,4748,5040,5373,5597,5736,6168
%N Sum of gcd(x, y) for 1 <= x, y <= n.
%C a(n) is also the entrywise 1-norm of the n X n GCD matrix.
%H Charles R Greathouse IV, <a href="/A018806/b018806.txt">Table of n, a(n) for n = 1..10000</a>
%F Sum_{k=1..n} phi(k)*(floor(n/k))^2. - _Vladeta Jovovic_, Nov 10 2002
%F a(n) ~ kn^2 log n, with k = 6/Pi^2. - _Charles R Greathouse IV_, Jun 21 2013
%F G.f.: Sum_{k >= 1} phi(k)*x^k*(1+x^k)/((1-x^k)^2*(1-x)). - _Robert Israel_, Jan 14 2015
%p N:= 1000 # to get a(1) to a(N)
%p g:= add(numtheory:-phi(k)*x^k*(1+x^k)/((1-x^k)^2*(1-x)),k=1..N):
%p S:= series(g, x, N+1):
%p seq(coeff(S,x,j), j=1..N); # _Robert Israel_, Jan 14 2015
%t Table[nn = n;Total[Level[Table[Table[GCD[i, j], {i, 1, nn}], {j, 1, nn}], {2}]], {n, 1, 48}] (* _Geoffrey Critzer_, Jan 14 2015 *)
%o (PARI) a(n)=2*sum(i=1,n,sum(j=1,i-1,gcd(i,j)))+n*(n+1)/2 \\ _Charles R Greathouse IV_, Jun 21 2013
%o (PARI) a(n)=sum(k=1,n,eulerphi(k)*(n\k)^2) \\ _Charles R Greathouse IV_, Jun 21 2013
%o (Python)
%o from sympy import totient
%o def A018806(n): return sum(totient(k)*(n//k)**2 for k in range(1,n+1)) # _Chai Wah Wu_, Aug 05 2024
%Y Cf. A000010, A018805, A064951, A344479.
%K nonn
%O 1,2
%A _David W. Wilson_