%I #21 Aug 31 2024 19:19:55

%S 1,4,9,25,12,49,20,121,169,28,45,289,361,63,44,529,52,841,60,961,99,

%T 68,175,1369,76,117,1681,84,1849,92,2209,153,2809,275,171,116,3481,

%U 3721,124,325,132,4489,207,140,5041,5329,148,539,156,6241,164,6889,425,172

%N a(n) = A073481(n)*A005117(n).

%C Define f(x) to be lpf(k)*k, where lpf(k) = A020639(k). This sequence contains the mappings of f(x) across squarefree numbers A005117.

%C a(1) = 1 by definition. 1 is the empty product and has no least prime factor.

%C Define sequence R_k = { m : rad(m) | k }, where rad(n) = A007947(n) and k > 1 is squarefree. Then the sequence k*R_k contains all numbers divisible by squarefree k that are also not divisible by any prime q coprime to k.

%C Plainly, k is the first term in the sequence k*R_k, because 1 is the first term in R_k. Hence a(n) is the second term in k*R_k for n > 1, since lpf(k) is the second term in R_k.

%H Michael De Vlieger, <a href="/A366786/b366786.txt">Table of n, a(n) for n = 1..10000</a>

%F a(n) = A065642(A005117(n)), n > 1.

%F a(n) = A285109(A005117(n)).

%F a(n) = A020639(A005117(n))*A005117(n).

%F For prime p, a(p) = p^2.

%F For composite squarefree k, a(k) = (p^2 * m) such that (p^2 * m) is in A364996.

%F Permutation of the union of {1}, A001248, and A366825.

%e Let b(n) = A005117(n).

%e a(2) = 4 = b(2)*lpf(b(2)) = 2*lpf(2) = 2*2. In {2*A000079}, 4 is the second term.

%e a(5) = 12 = b(5)*lpf(b(5)) = 6*lpf(6) = 6*2. In {6*A003586}, 12 is the second term..

%e a(11) = 45 = b(11)*lpf(b(11)) = 15*lpf(15) = 15*3. In {15*A003593}, 45 is the second term, etc.

%t nn = 120; s = Select[Range[nn], SquareFreeQ];

%t Array[#*FactorInteger[#][[1, 1]] &[s[[#]]] &, Length[s]]

%o (PARI) apply(x->(if (x==1,1, x*vecmin(factor(x)[,1]))), select(issquarefree, [1..150])) \\ _Michel Marcus_, Dec 17 2023

%o (Python)

%o from math import isqrt

%o from sympy import mobius, primefactors

%o def A366786(n):

%o def f(x): return n+x-sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1))

%o def bisection(f,kmin=0,kmax=1):

%o while f(kmax) > kmax: kmax <<= 1

%o while kmax-kmin > 1:

%o kmid = kmax+kmin>>1

%o if f(kmid) <= kmid:

%o kmax = kmid

%o else:

%o kmin = kmid

%o return kmax

%o return (m:=bisection(f))*min(primefactors(m),default=1) # _Chai Wah Wu_, Aug 31 2024

%Y Cf. A000040, A001248, A005117, A020639, A065642, A073481, A120944, A285109, A364996, A366807, A366825.

%K nonn,easy

%O 1,2

%A _Michael De Vlieger_, Dec 16 2023