OFFSET
1,2
COMMENTS
From Marius A. Burtea, Dec 26 2021: (Start)
If p is a prime number in A005384 (Sophie Germain prime), then it is a term. Indeed, p = p*1 and p + p + 1 = 2*p + 1 is prime.
All terms (m >= 2) are squarefree numbers (A005117). Indeed, if m = p^2*k, p >= 2, k >= 1, then m = p*(p*k) and p^2*k + p + p*k = p*(p*k + 1 + k ) is not prime. (End).
EXAMPLE
1 is in the sequence because 1 = 1*1 and 1*1 + 1 + 1 = 3 is prime;
30 is in the sequence because A038548(30) = 4 has 4 factorizations:
30 = 1*30 = 2*15 = 3*10 = 5*6 and
30 + 1 + 30 = 61 is prime;
30 + 2 + 15 = 47 is prime;
30 + 3 + 10 = 43 is prime;
30 + 5 + 6 = 41 is prime.
MAPLE
A350230 := proc(n)
local a, b ;
for a in numtheory[divisors](n) do
b := n/a ;
if not isprime(a*b+b+a) then
return false;
end if;
end do:
true ;
end proc:
for n from 1 to 500 do
if isA350230(n) then
printf("%d, ", n) ;
end if;
end do: # R. J. Mathar, Jan 24 2022
MATHEMATICA
t={}; Do[ds=Divisors[n]; If[EvenQ[Length[ds]], ok=True; k=1; While[k<=Length[ds]/2&&(ok=PrimeQ[ds[[k]]*ds[[-k]]+ds[[k]]+ds[[-k]]]), k++]; If[ok, AppendTo[t, n]]], {n, 2, 4000}]; t
PROG
(Python)
from sympy import divisors, isprime
def ok(n):
divs = divisors(n)
if n == 0: return False
return all(isprime(a*b+a+b) for a, b in ((d, n//d) for d in divs))
print([k for k in range(427) if ok(k)]) # Michael S. Branicky, Dec 21 2021
(PARI) isok(m) = sumdiv(m, d, isprime(m+d+m/d)) == numdiv(m); \\ Michel Marcus, Dec 25 2021
(Magma) [n:n in [1..450]|forall{d: d in Divisors(n)| IsPrime(n+d+n div d)}]; // Marius A. Burtea, Dec 26 2021
CROSSREFS
KEYWORD
nonn
AUTHOR
Michel Lagneau, Dec 21 2021
STATUS
approved