OFFSET
1,2
COMMENTS
The sequence of bi-unitary perfect numbers obeying a(n) = 2*n consists of only 6, 60, 90 [Wall].
Row sum of row n of the irregular table of the bi-unitary divisors, A222266.
LINKS
Reinhard Zumkeller, Table of n, a(n) for n = 1..10000
Krishnaswami Alladi, On arithmetic functions and divisors of higher order, J. Austral. Math. Soc. 23 (series A) (1977), 9-27.
József Sándor and Borislav Crstici, Perfect numbers: Old and new issues; perspectives, in Handbook of number theory, II, p. 45.
László Tóth, On the bi-unitary analogues of Euler's arithmetical function and the gcd-sum function J. Int. Seq. 12 (2009), Article 09.5.2.
Charles R. Wall, Bi-unitary perfect numbers, Proc. Am. Math. Soc. 33 (1) (1972), 39-42.
Eric Weisstein's World of Mathematics, Biunitary Divisor.
Tomohiro Yamada, 2 and 9 are the only biunitary superperfect numbers, arXiv:1705.00189 [math.NT], 2017.
FORMULA
Multiplicative with a(p^e) = (p^(e+1)-1)/(p-1) if e is odd, a(p^e) = (p^(e+1)-1)/(p-1) -p^(e/2) if e is even.
Dirichlet g.f.: zeta(s-1) * zeta(s) * zeta(2*s-1) * Product_{p prime} (1 - 2/p^(2*s-1) + 1/p^(3*s-2) + 1/p^(3*s-1) - 1/p^(4*s-2)). - Amiram Eldar, Aug 28 2023
EXAMPLE
The divisors of n=16 are d=1, 2, 4, 8 and 16. The greatest common unitary divisor of (1,16) is 1, of (2,8) is 1, of (4,4) is 4, of (8,2) is 1, of (16,1) is 1 (see A165430). So 1, 2, 8 and 16 are bi-unitary divisors of 16, which sum to a(16) = 1 + 2 + 8 + 16 = 27.
MAPLE
MATHEMATICA
f[n_] := Select[Divisors[n], Function[d, CoprimeQ[d, n/d]]]; Table[DivisorSum[n, # &, Last@ Intersection[f@ #, f[n/#]] == 1 &], {n, 76}] (* Michael De Vlieger, May 07 2017 *)
a[n_] := If[n==1, 1, Product[{p, e} = pe; If[OddQ[e], (p^(e+1)-1)/(p-1), ((p^(e+1)-1)/(p-1)-p^(e/2))], {pe, FactorInteger[n]}]]; Array[a, 80] (* Jean-François Alcover, Sep 22 2018 *)
PROG
(Haskell)
a188999 n = product $ zipWith f (a027748_row n) (a124010_row n) where
f p e = (p ^ (e + 1) - 1) `div` (p - 1) - (1 - m) * p ^ e' where
(e', m) = divMod e 2
-- Reinhard Zumkeller, Mar 04 2013
(PARI) udivs(n) = {my(d = divisors(n)); select(x->(gcd(x, n/x)==1), d); }
gcud(n, m) = vecmax(setintersect(udivs(n), udivs(m)));
biudivs(n) = select(x->(gcud(x, n/x)==1), divisors(n));
a(n) = vecsum(biudivs(n)); \\ Michel Marcus, May 07 2017
(PARI) a(n) = {f = factor(n); for (i=1, #f~, p = f[i, 1]; e = f[i, 2]; f[i, 1] = if (e % 2, (p^(e+1)-1)/(p-1), (p^(e+1)-1)/(p-1) -p^(e/2)); f[i, 2] = 1; ); factorback(f); } \\ Michel Marcus, Nov 09 2017
(Python)
from math import prod
from sympy import factorint
def A188999(n): return prod((p**(e+1)-1)//(p-1)-(0 if e&1 else p**(e>>1)) for p, e in factorint(n).items()) # Chai Wah Wu, Dec 28 2024
CROSSREFS
KEYWORD
mult,nonn,easy,changed
AUTHOR
R. J. Mathar, Apr 15 2011
STATUS
approved