cp's OEIS Frontend

A024621 Positions of nonprimes among the powers of primes (A000961).

%I A024621 #12 Nov 05 2024 20:39:47
%S A024621 1,4,7,8,11,15,16,19,24,28,33,42,43,45,54,69,71,79,87,92,118,121,137,
%T A024621 153,171,188,199,245,248,293,314,341,360,361,364,393,446,483,526,559,
%U A024621 605,651,699,719,750,857,894,930,935,1050,1079,1215,1305,1348,1436,1479,1514
%N A024621 Positions of nonprimes among the powers of primes (A000961).
%o A024621 (PARI) lista(nn) = {vec = vector(nn, i, i); powpr = select((i->((omega(i)==1) || (i==1))), vec); for (i = 1, #powpr, if (! isprime(powpr[i]), print1(i, ", ")););} \\ _Michel Marcus_, Oct 02 2013
%o A024621 (Python)
%o A024621 from sympy import primepi, integer_nthroot
%o A024621 def A024621(n):
%o A024621     def bisection(f,kmin=0,kmax=1):
%o A024621         while f(kmax) > kmax: kmax <<= 1
%o A024621         while kmax-kmin > 1:
%o A024621             kmid = kmax+kmin>>1
%o A024621             if f(kmid) <= kmid:
%o A024621                 kmax = kmid
%o A024621             else:
%o A024621                 kmin = kmid
%o A024621         return kmax
%o A024621     def f(x): return int(n+x-1-sum(primepi(integer_nthroot(x,k)[0]) for k in range(2,x.bit_length())))
%o A024621     return n+primepi(bisection(f,n,n)) # _Chai Wah Wu_, Nov 05 2024
%Y A024621 Complement of A024620.
%K A024621 nonn
%O A024621 1,2
%A A024621 _Clark Kimberling_