A006087 Unitary harmonic means H(n) of the unitary harmonic numbers (A006086).
1, 2, 3, 4, 4, 7, 7, 6, 9, 13, 10, 13, 10, 7, 11, 15, 10, 15, 9, 12, 7, 17, 12, 18, 16, 14, 19, 20, 19, 12, 15, 20, 10, 20, 18, 22, 19, 13, 12, 13, 17, 29, 18, 33, 20, 23, 29, 34, 23, 22, 31, 38, 24, 23, 38, 33, 37, 40, 19, 38, 24, 37, 29, 40, 22, 34, 24, 33
Offset: 1
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)
Programs
-
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
-
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 *)
-
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
Formula
Extensions
More terms from R. J. Mathar, Jun 06 2007
More terms from Herman Jamke (hermanjamke(AT)fastmail.fm), Mar 02 2008
Comments