A085971 - OEIS

%I #12 Aug 20 2024 13:19:50

%S 2,3,5,6,7,10,11,12,13,14,15,17,18,19,20,21,22,23,24,26,28,29,30,31,

%T 33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,50,51,52,53,54,55,56,

%U 57,58,59,60,61,62,63,65,66,67,68,69,70,71,72,73,74,75,76,77,78

%N Union of primes and numbers that are not prime powers (A000040, A024619).

%C Complement of A025475;

%C A085972(n) = Max{k: a(k)<=n};

%C different from A007916 and A052485, as a(28)=36;

%C A085818(a(n)) = A000040(n).

%F a(n) = n + o(sqrt n). - _Charles R Greathouse IV_, Oct 19 2015

%t With[{nn=100},Union[Join[Prime[Range[PrimePi[nn]]],DeleteCases[Range[2,80], _?(PrimePowerQ[#]&)]]]] (* _Harvey P. Dale_, May 15 2019 *)

%o (PARI) is(n)=isprimepower(n)<2 && n>1 \\ _Charles R Greathouse IV_, Oct 19 2015

%o (Python)

%o from sympy import primepi, integer_nthroot

%o def A085971(n):

%o     def f(x): return int(n+sum(primepi(integer_nthroot(x,k)[0]) for k in range(2,x.bit_length())))

%o     kmin, kmax = 1,2

%o     while f(kmax) >= kmax:

%o         kmax <<= 1

%o     while True:

%o         kmid = kmax+kmin>>1

%o         if f(kmid) < kmid:

%o             kmax = kmid

%o         else:

%o             kmin = kmid

%o         if kmax-kmin <= 1:

%o             break

%o     return kmax # _Chai Wah Wu_, Aug 20 2024

%K nonn,easy

%O 1,1

%A _Reinhard Zumkeller_, Jul 06 2003