OFFSET
1,2
COMMENTS
For a sequence A with terms a(1), a(2), a(3).... , let R(0) = 1, and for k >= 1 let R(k) = rad(a(1)*a(2)*...*a(k)). Define the Rad-transform of A to be R(n)/R(n-1); n >= 1, where rad is A007947. This sequence is the Rad transform of the squarefree numbers, A = A005117; see Example.
The sequence consists of only 1's and primes.
EXAMPLE
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 *)
PROG
(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
CROSSREFS
KEYWORD
nonn
AUTHOR
David James Sycamore, Apr 19 2024
EXTENSIONS
More terms from Michel Marcus, Apr 19 2024
STATUS
approved