OFFSET
1,2
LINKS
Amiram Eldar, Table of n, a(n) for n = 1..10000
FORMULA
a(n) = numerator(A018804(n)/n).
a(n) << n^(1+e) for any e > 0. a(n) > 1 for all n > 1. - Charles R Greathouse IV, Sep 08 2022
Sum_{k=1..n} a(k)/A356473(k) ~ (n/zeta(2)) * (log(n) + 2*gamma - 1 - zeta'(2)/zeta(2)), where gamma is Euler's constant (A001620). - Amiram Eldar, Dec 25 2024
EXAMPLE
For n = 3, the average of the gcd's is (gcd(1,3) + gcd(2,3) + gcd(3,3))/3 = (1 + 1 + 3)/3 = 5/3 and its numerator is a(3)=5.
MATHEMATICA
Table[Numerator[Sum[GCD[I, j], {j, 1, I}]/I], {I, 100}]
f[p_, e_] := e*(p - 1)/p + 1; a[n_] := Numerator[Times @@ f @@@ FactorInteger[n]]; Array[a, 100] (* Amiram Eldar, Apr 28 2023 *)
PROG
(Haskell) map numerator (map (\i -> sum (map (\j -> gcd i j) [1..i]) % i) [1..])
(PARI) a(n) = numerator(sum(i=1, n, gcd(i, n))/n); \\ Michel Marcus, Aug 08 2022
(PARI) a(n, f=factor(n))=my(k=prod(i=1, #f~, (f[i, 2]*(f[i, 1]-1)/f[i, 1] + 1)*f[i, 1]^f[i, 2])); k/gcd(k, n) \\ Charles R Greathouse IV, Sep 08 2022
(Python)
from math import prod, gcd
from sympy import factorint
def A356472(n):
f = factorint(n)
return (m:=prod((p-1)*e+p for p, e in f.items()))//gcd(prod(f), m) # Chai Wah Wu, Sep 08 2022
CROSSREFS
KEYWORD
easy,frac,nonn
AUTHOR
Matthias Kaak, Aug 08 2022
STATUS
approved
