|
%I M0452 #43 Mar 10 2023 09:09:55
%S 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,
%T 20,19,12,15,20,10,20,18,22,19,13,12,13,17,29,18,33,20,23,29,34,23,22,
%U 31,38,24,23,38,33,37,40,19,38,24,37,29,40,22,34,24,33
%N Unitary harmonic means H(n) of the unitary harmonic numbers (A006086).
%C 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).
%C Each term appears a finite number of times in the sequence (Hagis and Lord, 1975). - _Amiram Eldar_, Mar 10 2023
%D N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).
%H Donovan Johnson, <a href="/A006087/b006087.txt">Table of n, a(n) for n = 1..290</a>
%H Peter Hagis, Jr. and Graham Lord, <a href="https://doi.org/10.1090/S0002-9939-1975-0369231-9">Unitary harmonic numbers</a>, Proc. Amer. Math. Soc., 51 (1975), 1-7.
%H Peter Hagis, Jr. and Graham Lord, <a href="/A006086/a006086.pdf">Unitary harmonic numbers</a>, Proc. Amer. Math. Soc., 51 (1975), 1-7. (Annotated scanned copy)
%F a(n) = A103339(A006086(n)). - _Reinhard Zumkeller_, Mar 17 2012
%p 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
%t 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 *)
%t uh[n_] := n * Times @@ (2/(1 + Power @@@ FactorInteger[n])); uh[1] = 1; Select[Array[uh, 10^6], IntegerQ] (* _Amiram Eldar_, Mar 10 2023 *)
%o (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
%o (PARI) uhmean(n) = {my(f = factor(n)); n*prod(i=1, #f~, 2/(1+f[i, 1]^f[i, 2])); };
%o lista(kmax) = {my(uh); for(k = 1, kmax, uh = uhmean(k); if(denominator(uh) == 1, print1(uh, ", ")));} \\ _Amiram Eldar_, Mar 10 2023
%o (Haskell)
%o import Data.Ratio ((%), numerator, denominator)
%o a006087 n = a006087_list !! (n-1)
%o a006087_list = map numerator $ filter ((== 1) . denominator) $
%o map uhm [1..] where uhm n = (n * a034444 n) % (a034448 n)
%o -- _Reinhard Zumkeller_, Mar 17 2012
%Y Cf. A006086, A034444, A034448, A077610.
%K nonn,nice
%O 1,2
%A _N. J. A. Sloane_
%E More terms from _R. J. Mathar_, Jun 06 2007
%E More terms from Herman Jamke (hermanjamke(AT)fastmail.fm), Mar 02 2008
|
