A122410 a(n) is the sum of j's for those k's, 1 <= k <= n, where gcd(k,n) = p^j, p = prime.
0, 1, 1, 3, 1, 3, 1, 7, 4, 5, 1, 8, 1, 7, 6, 15, 1, 10, 1, 14, 8, 11, 1, 18, 6, 13, 13, 20, 1, 14, 1, 31, 12, 17, 10, 26, 1, 19, 14, 32, 1, 20, 1, 32, 22, 23, 1, 38, 8, 26, 18, 38, 1, 31, 14, 46, 20, 29, 1, 36, 1, 31, 30, 63, 16, 32, 1, 50, 24, 34, 1, 58, 1, 37, 32, 56, 16, 38, 1, 68, 40
Offset: 1
Keywords
Examples
The positive integers k, k <= 12, where gcd(k,12) = a power of a prime, are 1, 2, 3, 4, 8, 9 and 10; gcd(1,12) = p^0, gcd(2,12) = 2^1, gcd(3,12) = 3^1, gcd(4,12) = 2^2, gcd(8,12) = 2^2, gcd(9,12) = 3^1 and gcd(10,12) = 2^1. The sum of the exponents raising the primes is 0+1+1+2+2+1+1 = 8. So a(12) = 8.
Links
- Antti Karttunen, Table of n, a(n) for n = 1..65537
Programs
-
Mathematica
f[n_] := Plus @@ Last /@ Flatten[Select[FactorInteger[GCD[Range[n], n]], Length[ # ] == 1 &], 1]; Table[f[n], {n, 80}] (* Ray Chandler, Sep 06 2006 *) a[n_] := Module[{f = FactorInteger[n], p, e}, p = f[[;; , 1]]; e = f[[;; , 2]]; n * Times @@ (1-1/p) * Total[p*(1-1/p^e)/(p - 1)^2]]; a[1] = 0; Array[a, 100] (* Amiram Eldar, Jun 21 2025 *) -
PARI
A122410(n) = sum(k=1,n,isprimepower(gcd(n,k))); \\ Antti Karttunen, Feb 25 2018
-
PARI
a(n) = {my(f = factor(n), p = f[,1], e = f[,2]); eulerphi(f) * sum(i = 1, #p, p[i] * (1 - 1/p[i]^e[i]) / (p[i] - 1)^2);} \\ Amiram Eldar, Jun 21 2025
Formula
From Ridouane Oudra, Jun 06 2025: (Start)
a(p^m) = (p^m-1)/(p-1), for p prime and m >= 0.
a(n) = Sum_{p|n, p prime} phi(n/p), for n a squarefree.
More generally, for all n we have:
a(n) = Sum_{d|n, d is a prime power} A100995(d)*phi(n/d).
a(n) = Sum_{p|n, p prime} ((p^v(n,p)-1)/(p-1))*phi(n/p^v(n,p)), where p^v(n,p) is the highest power of p dividing n.
a(n) = phi(n) * Sum_{p|n, p prime} p*(1-p^(-v(n,p)))/((1-p)^2). (End)
Sum_{k=1..n} a(k) ~ c * zeta(2) * n^2 / 2, where c = Sum_{p prime} (p^2/(p^2-1)^2) = 0.68073222355480674093... . - Amiram Eldar, Jun 21 2025
Extensions
Extended by Ray Chandler, Sep 06 2006