(PARI)
is(n) = {
my(p = gcd(n, eulerphi(n)), f, g);
if (isprime(p), return(n % p^2 == 0 && isprime(gcd(p+1, n))));
if (omega(p) != 2 || !issquarefree(n), return(0));
f = factor(n); g = factor(p);
1 == g[2, 1] % g[1, 1] &&
1 == sum(k=1, matsize(f)[1], f[k, 1] % g[1, 1] == 1) &&
1 == sum(k=1, matsize(f)[1], f[k, 1] % g[2, 1] == 1);
};
seq(N) = {
my(a = vector(N), k=0, n=1);
while(k < N, if(is(n), a[k++]=n); n++); a;
};