(PARI) isok(n) = my(f = factor(n^2-1)); #f~ == primepi(f[#f~, 1]); \\
Michel Marcus, Jul 02 2016
(Python)
from sympy import prime, primefactors
def is_ok(m):
pf = primefactors(m*m - 1)
return prime(len(pf)) == pf[-1]
print([ m for m in range(3, 1270000, 2) if is_ok(m)])
(Python)
from itertools import islice
from heapq import heappop, heappush
from sympy import integer_nthroot, nextprime
def
A194099_gen(): # generator of terms
h, hset = [(1, (1, ))], {1}
while True:
m, ps = heappop(h)
a, b = integer_nthroot(m+1, 2)
if b:
yield a
for p in ps:
mp = m*p
if mp not in hset:
heappush(h, (mp, ps))
hset.add(mp)
q = nextprime(max(ps, default=1))
mp = m*q
if mp not in hset:
heappush(h, (mp, (ps+(q, ))))
hset.add(mp)
