A372112 Rad-transform of the squarefree numbers A005117 (see Comments).
1, 2, 3, 5, 1, 7, 1, 11, 13, 1, 1, 17, 19, 1, 1, 23, 1, 29, 1, 31, 1, 1, 1, 37, 1, 1, 41, 1, 43, 1, 47, 1, 53, 1, 1, 1, 59, 61, 1, 1, 1, 67, 1, 1, 71, 73, 1, 1, 1, 79, 1, 83, 1, 1, 1, 89, 1, 1, 1, 1, 97, 101, 1, 103, 1, 1, 107, 109, 1, 1, 113, 1, 1, 1, 1, 1, 1, 127, 1, 1
Offset: 1
Keywords
Examples
For A005117, R(k) (k >= 0) is 1 U A128040. That is, 1,1,2,6,30,30,210,... from which: a(1) = 1/1 = 1, a(2) = 2/1 = 2, a(3) = 6/2 = 3, 4(4) = 30/6 = 5, a(5) = 30/30 = 1, and so on. Note that the first term of R(k) is 1, the empty product (product of the first 0 terms of A005117).
Links
- Michael De Vlieger, Table of n, a(n) for n = 1..10000
Programs
-
Mathematica
k = r = s = 1; f[x_] := f[x] = Times @@ FactorInteger[x][[All, 1]]; Reap[Do[While[! SquareFreeQ[k], k++]; s = f[s*f[k]]; Sow[s/r]; r = s; k++, {n, 120}] ][[-1, 1]] (* Michael De Vlieger, Apr 19 2024 *) -
PARI
rad(n) = factorback(factorint(n)[, 1]); \\ A007947 lista(nn) = my(v = select(issquarefree, [1..nn])); my(w = vector(#v, k, rad(prod(i=1, k, v[i])))); concat(1, vector(#w-1, k, w[k+1]/w[k])); \\ Michel Marcus, Apr 19 2024
-
Python
from math import isqrt from sympy import mobius, isprime def A372112(n): def f(x): return int(n-sum(mobius(k)*(x//k**2) for k in range(2, isqrt(x)+1))) m, k = n, f(n) while m != k: m, k = k, f(k) return m if isprime(m) else 1 # Chai Wah Wu, Dec 23 2024
Extensions
More terms from Michel Marcus, Apr 19 2024
Comments