%I #9 Mar 24 2023 11:39:06
%S 1,6,60,270,420,630,672,2970,5460,8190,10080,22848,30240,99792,136500,
%T 172900,204750,208656,245700,249480,312480,332640,342720,385560,
%U 491400,695520,708288,791700,819000,861840,1028160,1037400,1187550,1228500,1421280,1528800,1571328
%N Bi-unitary arithmetic numbers k whose mean bi-unitary divisor is a bi-unitary divisor of k.
%C Also, bi-unitary harmonic numbers k whose harmonic mean of the bi-unitary divisors of k is a bi-unitary divisor of k.
%H Amiram Eldar, <a href="/A361787/b361787.txt">Table of n, a(n) for n = 1..236</a>
%e 6 is a term since the arithmetic mean of its bi-unitary divisors, {1, 2, 3, 6}, is 3, and 3 is also a bi-unitary divisor of 6.
%e 60 is a term since the arithmetic mean of its bi-unitary divisors, {1, 3, 4, 5, 12, 15, 20, 60}, is 15, and 15 is also a bi-unitary divisor of 60.
%t biudivQ[f_, d_] := AllTrue[f, OddQ[Last[#]] || IntegerExponent[d, First[#]] != Last[#]/2 &]; biuDivs[n_] := Module[{d = Divisors[n], f = FactorInteger[n]}, Select[d, biudivQ[f, #] &]]; Select[Range[10^5], IntegerQ[(r = Mean[(i = biuDivs[#])])] && MemberQ[i, r] &]
%o (PARI) isbdiv(f, d) = {for (i=1, #f~, if(f[i, 2]%2 == 0 && valuation(d, f[i, 1]) == f[i, 2]/2, return(0))); 1;}
%o is(n) = {my(f = factor(n), r, p, e); r = prod(i=1, #f~, p = f[i, 1]; e = f[i, 2]; if(e%2, (p^(e+1)-1)/((e + 1)*(p-1)), ((p^(e+1)-1)/(p-1)-p^(e/2))/e)); denominator(r) == 1 && n%r==0 && isbdiv(f, r); }
%Y Subsequence of A286325 and A361786.
%Y Similar sequence: A007340, A353039, A361387.
%Y Cf. A188999, A222266, A286324.
%K nonn
%O 1,2
%A _Amiram Eldar_, Mar 24 2023