%I #64 Aug 05 2024 16:03:48
%S 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,17,18,19,20,21,22,23,24,25,26,27,
%T 28,29,30,31,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,49,50,51,52,
%U 53,54,55,56,57,58,59,60,61,62,63,65,66,67,68,69,70,71,72,73,74,75,76
%N Biquadratefree numbers: numbers that are not divisible by any 4th power greater than 1.
%C Differs from A023809 at entries 0, 81, 162, 225, 226, etc. - _R. J. Mathar_, Oct 18 2008
%C Density is 1/zeta(4) = A215267 = 0.923938.... - _Charles R Greathouse IV_, Sep 02 2015
%C The Schnirelmann density of the biquadratefree numbers is 145/157 (Orr, 1969). - _Amiram Eldar_, Mar 12 2021
%C This sequence has arbitrarily large gaps and hence is not a Beatty sequence. - _Charles R Greathouse IV_, Jan 27 2022
%H Amiram Eldar, <a href="/A046100/b046100.txt">Table of n, a(n) for n = 1..10000</a> (terms 1..1000 from T. D. Noe)
%H Richard C. Orr, <a href="https://doi.org/10.1112/jlms/s1-44.1.313">On the Schnirelmann density of the sequence of k-free integers</a>, Journal of the London Mathematical Society, Vol. 1, No. 1 (1969), pp. 313-319.
%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/Biquadratefree.html">Biquadratefree</a>.
%F A051903(a(n)) < 4. - _Reinhard Zumkeller_, Sep 03 2015
%F Sum_{n>=1} 1/a(n)^s = zeta(s)/zeta(4*s), for s > 1. - _Amiram Eldar_, Dec 27 2022
%p A046100 := proc(n)
%p option remember;
%p local a,p,is4free;
%p if n = 1 then
%p return 1;
%p else
%p for a from procname(n-1)+1 do
%p is4free := true ;
%p for p in ifactors(a)[2] do
%p if op(2,p) >= 4 then
%p is4free := false;
%p break;
%p end if;
%p end do:
%p if is4free then
%p return a;
%p end if;
%p end do:
%p end if;
%p end proc: # _R. J. Mathar_, Aug 08 2012
%t lst={};Do[a=0;Do[If[FactorInteger[m][[n, 2]]>4, a=1], {n, Length[FactorInteger[m]]}];If[a!=1, AppendTo[lst, m]], {m, 5!}];lst (* _Vladimir Joseph Stephan Orlovsky_, Sep 27 2008 *)
%t Select[Range[100],Max[FactorInteger[#][[;;,2]]]<4&] (* _Harvey P. Dale_, Jul 13 2023 *)
%o (PARI) is(n)=n==1 || vecmax(factor(n)[,2])<4 \\ _Charles R Greathouse IV_, Jun 16 2012
%o (Sage)
%o def is_biquadratefree(n):
%o return all(c[1] < 4 for c in n.factor())
%o def A046100_list(n): return [i for i in (1..n) if is_biquadratefree(i)]
%o A046100_list(76) # _Peter Luschny_, Aug 08 2012
%o (Haskell)
%o a046100 n = a046100_list !! (n-1)
%o a046100_list = filter ((< 4) . a051903) [1..]
%o -- _Reinhard Zumkeller_, Sep 03 2015
%o (Python)
%o from sympy import mobius, integer_nthroot
%o def A046100(n):
%o def f(x): return n+x-sum(mobius(k)*(x//k**4) for k in range(1, integer_nthroot(x,4)[0]+1))
%o m, k = n, f(n)
%o while m != k:
%o m, k = k, f(k)
%o return m # _Chai Wah Wu_, Aug 05 2024
%Y Cf. A046101, A005117 (2-free), A004709 (3-free).
%Y Cf. A023809, A051903, A215267.
%Y Subsequence of A209061.
%K nonn
%O 1,2
%A _Eric W. Weisstein_
%E Name edited by _Amiram Eldar_, Jul 29 2024