OFFSET
1,2
COMMENTS
A divisor of n is called infinitary if it is a product of divisors of the form p^{y_a 2^a}, where p^y is a prime power dividing n and sum_a y_a 2^a is the binary representation of y.
Multiplicative: If e = Sum_{k >= 0} d_k 2^k (binary representation of e), then a(p^e) = Product_{k >= 0} (p^(2^k*{d_k+1}) - 1)/(p^(2^k) - 1). - Christian G. Bower and Mitch Harris, May 20 2005 [This means there is a factor p^2^k + 1 if d_k = 1, else the factor is 1. - M. F. Hasler, Oct 20 2022]
This sequence is an infinitary analog of the Dedekind psi function A001615. Indeed, a(n) = Product_{q in Q_n}(q+1) = n*Product_{q in Q_n} (1+1/q), where {q} are terms of A050376 and Q_n is the set of distinct q's whose product is n. - Vladimir Shevelev, Apr 01 2014
LINKS
Amiram Eldar, Table of n, a(n) for n = 1..10000 (terms 1..7417 from R. J. Mathar)
Graeme L. Cohen, On an integer's infinitary divisors, Math. Comp. 54 (189) (1990) 395-411.
Steven R. Finch, Unitarism and Infinitarism, February 25, 2004. [Cached copy, with permission of the author]
J. O. M. Pedersen, Tables of Aliquot Cycles [Broken link]
J. O. M. Pedersen, Tables of Aliquot Cycles [Via Internet Archive Wayback-Machine]
J. O. M. Pedersen, Tables of Aliquot Cycles [Cached copy, pdf file only]
Tomohiro Yamada, Infinitary superperfect numbers, arXiv:1705.10933 [math.NT], 2017.
FORMULA
Let n = Product(q_i) where {q_i} is a set of distinct terms of A050376. Then a(n) = Product(q_i + 1). - Vladimir Shevelev, Feb 19 2011
If n is squarefree, then a(n) = A001615(n). - Vladimir Shevelev, Apr 01 2014
a(n) = Sum_{k>=1} A077609(n,k). - R. J. Mathar, Oct 04 2017
a(n) = A126168(n)+n. - R. J. Mathar, Oct 05 2017
Multiplicative with a(p^e) = Product{k >= 0, e_k = 1} p^2^k + 1, where e = Sum e_k 2^k, i.e., e_k is bit k of e. - M. F. Hasler, Oct 20 2022
a(n) = iphi(n^2)/iphi(n), where iphi(n) = A091732(n). - Amiram Eldar, Sep 21 2024
EXAMPLE
If n = 8: 8 = 2^3 = 2^"11" (writing 3 in binary) so the infinitary divisors are 2^"00" = 1, 2^"01" = 2, 2^"10" = 4 and 2^"11" = 8; so a(8) = 1+2+4+8 = 15.
n = 90 = 2*5*9, where 2, 5, 9 are in A050376; so a(n) = 3*6*10 = 180. - Vladimir Shevelev, Feb 19 2011
MAPLE
isidiv := proc(d, n)
local n2, d2, p, j;
if n mod d <> 0 then
return false;
end if;
for p in numtheory[factorset](n) do
padic[ordp](n, p) ;
n2 := convert(%, base, 2) ;
padic[ordp](d, p) ;
d2 := convert(%, base, 2) ;
for j from 1 to nops(d2) do
if op(j, n2) = 0 and op(j, d2) <> 0 then
return false;
end if;
end do:
end do;
return true;
end proc:
idivisors := proc(n)
local a, d;
a := {} ;
for d in numtheory[divisors](n) do
if isidiv(d, n) then
a := a union {d} ;
end if;
end do:
a ;
end proc:
A049417 := proc(n)
local d;
add(d, d=idivisors(n)) ;
end proc:
seq(A049417(n), n=1..100) ; # R. J. Mathar, Feb 19 2011
MATHEMATICA
bitty[k_] := Union[Flatten[Outer[Plus, Sequence @@ ({0, #1} & ) /@ Union[2^Range[0, Floor[Log[2, k]]]*Reverse[IntegerDigits[k, 2]]]]]]; Table[Plus@@((Times @@ (First[it]^(#1 /. z -> List)) & ) /@ Flatten[Outer[z, Sequence @@ bitty /@ Last[it = Transpose[FactorInteger[k]]], 1]]), {k, 2, 120}]
(* Second program: *)
a[n_] := If[n == 1, 1, Sort @ Flatten @ Outer[ Times, Sequence @@ (FactorInteger[n] /. {p_, m_Integer} :> p^Select[Range[0, m], BitOr[m, #] == m &])]] // Total;
Array[a, 100] (* Jean-François Alcover, Mar 23 2020, after Paul Abbott in A077609 *)
PROG
(PARI) A049417(n) = {my(b, f=factorint(n)); prod(k=1, #f[, 2], b = binary(f[k, 2]); prod(j=1, #b, if(b[j], 1+f[k, 1]^(2^(#b-j)), 1)))} \\ Andrew Lelechenko, Apr 22 2014
(PARI) isigma(n)=vecprod([vecprod([f[1]^2^k+1|k<-[0..exponent(f[2])], bittest(f[2], k)])|f<-factor(n)~]) \\ M. F. Hasler, Oct 20 2022
(Haskell)
a049417 1 = 1
a049417 n = product $ zipWith f (a027748_row n) (a124010_row n) where
f p e = product $ zipWith div
(map (subtract 1 . (p ^)) $
zipWith (*) a000079_list $ map (+ 1) $ a030308_row e)
(map (subtract 1 . (p ^)) a000079_list)
-- Reinhard Zumkeller, Sep 18 2015
(Python)
from math import prod
from sympy import factorint
def A049417(n): return prod(p**(1<<i)+1 for p, e in factorint(n).items() for i, j in enumerate(bin(e)[-1:1:-1]) if j=='1') # Chai Wah Wu, Jul 11 2024
CROSSREFS
KEYWORD
nonn,mult
AUTHOR
Yasutoshi Kohmoto, Dec 11 1999
EXTENSIONS
More terms from Wouter Meeussen, Sep 02 2001
STATUS
approved