A246716 - OEIS

%I #85 Sep 22 2025 16:01:18

%S 1,2,3,4,5,7,8,9,11,12,13,16,17,18,19,20,23,24,25,27,28,29,30,31,32,

%T 36,37,40,41,42,43,44,45,47,48,49,50,52,53,54,56,59,60,61,63,64,66,67,

%U 68,70,71,72,73,75,76,78,79,80,81,83,84,88,89,90,92,96,97,98,99,100

%N Positive numbers that are not the product of (exactly) two distinct primes.

%C Non-disjoint union of A100959 and A000961. Disjoint union of A100959 and A001248.

%C Complement of A006881, then inheriting the "opposite" of the properties of A006881.

%C a(n+1) - a(n) <= 4 (gap upper bound) - (that is because among four consecutive integers there is always a multiple of 4, then there is a number which is not the product of two distinct primes). E.g., a(26)-a(25) = a(62)-a(61) = 4. Is it true that for any k <= 4 there are infinitely many numbers n such that a(n+1) - a(n) = k?

%C If r = A006881(n+1) - A006881(n) - 1 > 1, it indicates that there are r terms of (a(j)) starting with j = A006881(n) - n + 1 which are consecutive integers. E.g., A006881(8) - A006881(7) - 1 = 6, so there are 6 consecutive terms in (a(j)), starting with j = A006881(7) - 7 + 1 = 20.

%H Giuseppe Coppoletta, <a href="/A246716/b246716.txt">Table of n, a(n) for n = 1..10000</a>

%e 7 is a term because 7 is prime, so it has only one prime divisor.

%e 8 and 9 are terms because neither of them has two distinct prime divisors.

%e 30 is a term because it is the product of three primes.

%e But 35 is not a term because it is the product of two distinct primes.

%p filter:= n -> map(t -> t[2],ifactors(n)[2]) <> [1,1]:

%p select(filter, [$1..1000]); # _Robert Israel_, Nov 02 2014

%t Select[Range[125], Not[PrimeOmega[#] == PrimeNu[#] == 2] &] (* _Alonso del Arte_, Nov 03 2014 *)

%o (PARI) isok(n) = (omega(n)!=2) || (bigomega(n) != 2); \\ _Michel Marcus_, Nov 01 2014

%o (Magma) [n: n in [1..100] | #PrimeDivisors(n) ne 2 or &*[t[2]: t in Factorization(n)] ne 1]; // _Bruno Berselli_, Nov 12 2014

%o (SageMath)

%o def A246716_list(n) :

%o     R = []

%o     for i in (1..n) :

%o         d = prime_divisors(i)

%o         if len(d) != 2 or d[0]*d[1] != i : R.append(i)

%o     return R

%o A246716_list(100)

%o (SageMath) [n for n in (1..100) if sloane.A001221(n)!=2 or sloane.A001222(n)!=2] # _Giuseppe Coppoletta_, Jan 19 2015

%o (Python)

%o from math import isqrt

%o from sympy import primepi, primerange

%o def A246716(n):

%o     def f(x): return int(n-(t:=primepi(s:=isqrt(x)))-(t*(t-1)>>1)+sum(primepi(x//k) for k in primerange(1, s+1)))

%o     def bisection(f,kmin=0,kmax=1):

%o         while f(kmax) > kmax: kmax <<= 1

%o         while kmax-kmin > 1:

%o             kmid = kmax+kmin>>1

%o             if f(kmid) <= kmid:

%o                 kmax = kmid

%o             else:

%o                 kmin = kmid

%o         return kmax

%o     return bisection(f) # _Chai Wah Wu_, Aug 30 2024

%Y Cf. A001358, A100959, A006881, A007774, A001221, A001222, A007304, A000977, A001248, A000961.

%K nonn,easy

%O 1,2

%A _Giuseppe Coppoletta_, Nov 01 2014