OFFSET
1,2
COMMENTS
From Robert Israel, Dec 29 2022: (Start)
If n is prime, a(n) = n.
If n is an odd prime, a(n^2) = n.
If p and q are distinct primes with p | 2*q-1, then a(p*q) = q, except in the case of 3*5, where both 3 | 2*5-1 and 5 | 2*3-1, and a(3*5) = 1.
For semiprimes p*q where p <> q, p does not divide 2*q-1 and q does not divide 2*p-1, a(p*q) = p*q. (End)
LINKS
Robert Israel, Table of n, a(n) for n = 1..10000
FORMULA
a(n) = denominator of A018804(n)/n.
a(n) divides n, so in particular a(n) <= n. - Charles R Greathouse IV, Sep 08 2022
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 which has denominator a(3)=3.
MAPLE
f:= proc(n) local i; denom(add(igcd(i, n), i=1..n)/n) end proc:
map(f, [$1..100]); # Robert Israel, Dec 29 2022
MATHEMATICA
Table[Denominator[Sum[GCD[I, j], {j, 1, I}]/I], {I, 100}]
f[p_, e_] := e*(p - 1)/p + 1; a[n_] := Denominator[Times @@ f @@@ FactorInteger[n]]; Array[a, 100] (* Amiram Eldar, Apr 28 2023 *)
PROG
(Haskell)
map denominator (map (\i -> sum (map (\j -> gcd i j) [1..i]) % i) [1..])
(PARI) a(n) = denominator(sum(i=1, n, gcd(i, n))/n); \\ Michel Marcus, Aug 08 2022
(PARI) a(n, f=factor(n))=n/gcd(prod(i=1, #f~, (f[i, 2]*(f[i, 1]-1)/f[i, 1] + 1)*f[i, 1]^f[i, 2]), n) \\ Charles R Greathouse IV, Sep 08 2022
(Python)
from math import gcd, prod
from sympy import factorint
def A356473(n): return (p:=prod(f:=factorint(n)))//gcd(p, prod((p-1)*e+p for p, e in f.items())) # Chai Wah Wu, Sep 08 2022
CROSSREFS
KEYWORD
AUTHOR
Matthias Kaak, Aug 08 2022
STATUS
approved