OFFSET
1,2
COMMENTS
A self-inverse permutation on the positive integers: a(a(n)) = n.
LINKS
Alois P. Heinz, Table of n, a(n) for n = 1..10000
FORMULA
Sum_{k=1..n} a(k) ~ c * n^3, where c = (1/3) * Product_{p prime} ((p-1)*(p^6 + q(p) +(p^3-1)*q(p)^2))/(p^7 - p*q(p)^2) = 0.3120270364..., where q(p) = nextprime(p) = A151800(p) if p has an odd index, and q(p) = prevprime(p) = A151799(p) otherwise. - Amiram Eldar, Sep 17 2023
EXAMPLE
n = 2^3 => a(n) = 3^4 = 81.
n = 3^2 => a(n) = 2^1 = 2.
MAPLE
a:= n-> mul(ithprime(i[1])^i[2], i=map(x->map(y->y-1+2*irem(y, 2),
[numtheory[pi](x[1]), x[2]]), ifactors(n)[2])):
seq(a(n), n=1..100); # Alois P. Heinz, Jul 17 2013
MATHEMATICA
a[n_] := If[n == 1, 1, Product[{p, e} = pe; Prime[BitXor[PrimePi[p] - 1, 1] + 1]^(BitXor[e - 1, 1] + 1), {pe, FactorInteger[n]}]];
Array[a, 100] (* Jean-François Alcover, May 31 2019, after Andrew Howroyd *)
PROG
(PARI) a(n)={my(f=factor(n)); prod(i=1, #f~, my(p=f[i, 1], e=f[i, 2]); prime( bitxor( primepi(p)-1, 1)+1)^(bitxor(e-1, 1)+1))} \\ Andrew Howroyd, Jul 23 2018
(Python)
primes = [2]*2
primes[1] = 3
def addPrime(k):
for p in primes:
if k%p==0: return
if p*p > k: break
primes.append(k)
for n in range(5, 1000000, 6):
addPrime(n)
addPrime(n+2)
for n in range(1, 99):
p = 1
j = n
i = 0
while j>1:
e = 0
while j % primes[i] == 0:
j /= primes[i]
e+=1
if e:
e = ((e-1)^1) + 1
p*= primes[i^1]**e
i += 1
print(str(p), end=', ')
CROSSREFS
KEYWORD
nonn,easy,mult
AUTHOR
Alex Ratushnyak, Jul 07 2013
EXTENSIONS
Keyword:mult added by Andrew Howroyd, Jul 23 2018
STATUS
approved