|
%I #18 Sep 22 2021 07:41:47
%S 9,22,25,39,62,69,77,87,91,94,95,106,115,119,121,122,133,134,142,146,
%T 159,183,187,202,213,214,218,219,226,235,237,249,253,259,262,265,274,
%U 287,289,291,299,303,305,309,314,335,362,381,386,393,403,411,417,422,446,458,469,473,489,501,502,505
%N Semiprimes (A001358) p*q such that p*q+p+q is also a semiprime.
%H Robert Israel, <a href="/A330477/b330477.txt">Table of n, a(n) for n = 1..10000</a>
%e a(3) = 25 is a member because 25 = 5*5 and 25+5+5 = 5*7 is also a semiprime.
%p N:= 1000:
%p Primes:= select(isprime, [2,seq(i,i=3..N)]):
%p SP:= sort([seq(seq([p,q],q=select(t -> t >= p and p*t<=N, Primes)),p=Primes)],(a,b) -> a[1]*a[2]<b[1]*b[2]):
%p map(t -> t[1]*t[2], select(t -> numtheory:-bigomega(t[1]*t[2]+t[1]+t[2])=2, SP));
%t Select[Union@ Apply[Join, Table[Flatten@{p #, Sort[{p, #}]} & /@ Prime@ Range@ PrimePi@ Floor[Max[#]/p], {p, #}]] &@ Prime@ Range@ 97, PrimeOmega[Total@ #] == 2 &][[All, 1]] (* _Michael De Vlieger_, Dec 15 2019 *)
%o (PARI) issemi(n)=bigomega(n)==2
%o list(lim)=my(v=List()); forprime(p=2, sqrtint(lim\=1), forprime(q=p, lim\p, if(issemi(p*q+p+q), listput(v,p*q)))); Set(v) \\ _Charles R Greathouse IV_, Dec 16 2019
%o (Python)
%o from sympy import factorint
%o def is_semiprime(n): return sum(e for e in factorint(n).values()) == 2
%o def ok(n):
%o f = factorint(n, multiple=True)
%o if len(f) != 2: return False
%o p, q = f
%o return len(factorint(p*q + p + q, multiple=True)) == 2
%o print(list(filter(ok, range(506)))) # _Michael S. Branicky_, Sep 22 2021
%Y Cf. A001358.
%Y Contains A108570.
%K nonn
%O 1,1
%A _J. M. Bergot_ and _Robert Israel_, Dec 15 2019
|
