OFFSET
1,2
COMMENTS
Let d(n) and sigma(n) be number and sum of unitary divisors of n; then unitary harmonic mean of unitary divisors is H(n)=n*d(n)/sigma(n).
Each term appears a finite number of times in the sequence (Hagis and Lord, 1975). - Amiram Eldar, Mar 10 2023
REFERENCES
N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).
LINKS
Donovan Johnson, Table of n, a(n) for n = 1..290
Peter Hagis, Jr. and Graham Lord, Unitary harmonic numbers, Proc. Amer. Math. Soc., 51 (1975), 1-7.
Peter Hagis, Jr. and Graham Lord, Unitary harmonic numbers, Proc. Amer. Math. Soc., 51 (1975), 1-7. (Annotated scanned copy)
FORMULA
MAPLE
A034444 := proc(n) 2^nops(ifactors(n)[2]) ; end: A034448 := proc(n) local ans, i, ifs ; ans :=1 ; ifs := ifactors(n)[2] ; for i from 1 to nops(ifs) do ans := ans*(1+ifs[i][1]^ifs[i][2]) ; od ; RETURN(ans) ; end: A006086 := proc(n) n*A034444(n)/A034448(n) ; end: for n from 1 to 5000000 do uhn := A006086(n) : if type(uhn, 'integer') then printf("%d, ", uhn) ; fi ; od : # R. J. Mathar, Jun 06 2007
MATHEMATICA
ud[n_] := 2^PrimeNu[n]; usigma[n_] := Sum[ If[ GCD[d, n/d] == 1, d, 0], {d, Divisors[n]}]; a[n_] := n*ud[n]/usigma[n]; a[1] = 1; Reap[ Do[ If[ IntegerQ[h = a[n]], Print[h]; Sow[h]], {n, 1, 10^7}]][[2, 1]] (* Jean-François Alcover, May 16 2013 *)
uh[n_] := n * Times @@ (2/(1 + Power @@@ FactorInteger[n])); uh[1] = 1; Select[Array[uh, 10^6], IntegerQ] (* Amiram Eldar, Mar 10 2023 *)
PROG
(PARI) {ud(n)=2^omega(n)} {sud(n) = sumdiv(n, d, if(gcd(d, n/d)==1, d))} {H(n)=n*ud(n)/sud(n)} for(n=1, 10000000, if(((n*ud(n))%sud(n))==0, print1(H(n)", "))) \\ Herman Jamke (hermanjamke(AT)fastmail.fm), Mar 02 2008
(PARI) uhmean(n) = {my(f = factor(n)); n*prod(i=1, #f~, 2/(1+f[i, 1]^f[i, 2])); };
lista(kmax) = {my(uh); for(k = 1, kmax, uh = uhmean(k); if(denominator(uh) == 1, print1(uh, ", "))); } \\ Amiram Eldar, Mar 10 2023
(Haskell)
import Data.Ratio ((%), numerator, denominator)
a006087 n = a006087_list !! (n-1)
a006087_list = map numerator $ filter ((== 1) . denominator) $
map uhm [1..] where uhm n = (n * a034444 n) % (a034448 n)
-- Reinhard Zumkeller, Mar 17 2012
CROSSREFS
KEYWORD
nonn,nice
AUTHOR
EXTENSIONS
More terms from R. J. Mathar, Jun 06 2007
More terms from Herman Jamke (hermanjamke(AT)fastmail.fm), Mar 02 2008
STATUS
approved