%I #44 Jan 19 2025 15:24:00

%S 1,2,3,4,5,6,6,7,8,10,9,12,10,12,15,11,12,16,13,20,18,18,14,21,15,20,

%T 16,24,17,30,18,19,27,24,30,32,20,26,30,35,21,36,22,36,40,28,23,33,24,

%U 30,36,40,25,32,45,42,39,34,26,60,27,36,48,28,50,54,29,48,42,60,30,56,31

%N In canonical prime factorization of n replace p^e with its index in A000961.

%C The definition of the sequence has been corrected, given that it uses A095874, the indices in the list A000961 of "powers of primes" starting with A000961(1) = 1, rather than A322981, index of p^e in the list of prime powers A246655, as written in the original definition. See A333235 for the variant of this sequence which uses A322981 and A246655 instead, maybe the more natural choice given that the factorization of integers consists of prime powers > 1. - _M. F. Hasler_, Jun 15 2021

%H Ivan Neretin, <a href="/A097621/b097621.txt">Table of n, a(n) for n = 1..10000</a>

%F Multiplicative with: p^e -> A095874(p^e), p prime.

%F a(A000961(n)) = n; a(a(n)) = A097622(n); a(a(a(n))) = A097623(n);

%F a(n) <= n; a(n) = n iff 60 mod n = 0: a(A018266(n)) = A018266(n);

%F a(A097624(n)) = n and a(m) <> n for n < A097624(n).

%e n=600 = 2^3 * 3 * 5^2 -> A095874(8)*A095874(3)*A095874(25) = 7 * 3 * 15 = 315.

%p N:= 1000: # to get a(1) to a(N)

%p Primes:= select(isprime,[2, seq(2*i+1,i=1..(N-1)/2)]):

%p PP:= sort([1,seq(seq(p^k, k=1..floor(log[p](N))),p=Primes)]):

%p for n from 1 to nops(PP) do B[PP[n]]:= n od:

%p seq(mul(B[f[1]^f[2]],f=ifactors(n)[2]),n=1..N); # _Robert Israel_, Sep 02 2015

%t pp = Select[Range@100, Length[FactorInteger[#]] == 1 &]; a = Table[Times @@ (Position[pp, #][[1, 1]] & /@ (#[[1]]^#[[2]] & /@ FactorInteger[n])), {n, 73}] (* _Ivan Neretin_, Sep 02 2015 *)

%o (PARI) f(n) = if(isprimepower(n), sum(i=1, logint(n, 2), primepi(sqrtnint(n, i)))+1, n==1); \\ A095874

%o a(n) = my(fr=factor(n)); for (k=1, #fr~, fr[k,1] = f(fr[k,1]^fr[k,2]); fr[k,2] = 1); factorback(fr); \\ _Michel Marcus_, May 29 2021

%o A097621(n)=vecprod([A095874(f[1]^f[2])|f<-factor(n)~]) \\ _M. F. Hasler_, Jun 15 2021

%o (Python)

%o from math import prod

%o from sympy import primepi, integer_nthroot, factorint

%o def A097621(n): return prod(1+int(primepi(m:=p**e)+sum(primepi(integer_nthroot(m,k)[0]) for k in range(2,m.bit_length()))) for p,e in factorint(n).items()) # _Chai Wah Wu_, Jan 19 2025

%Y Cf. A000961 (powers of primes), A246655 (prime powers), A003963, A018266, A095874 (index of n = p^e in A000961).

%Y Cf. A097622, A097623, A097624.

%Y Cf. A322981 (index of n = p^e in A246655), A333235 (variant of this sequence).

%K nonn,mult

%O 1,2

%A _Reinhard Zumkeller_, Aug 17 2004

%E Definition corrected by _M. F. Hasler_, Jun 16 2021

%E Example corrected by _Ray Chandler_, Jun 30 2021