[seqfan] PrimeQ and Gaussian integers

[israel at math.ubc.ca]

Today I submitted a correction to A058770, but I suspect the issue might 
affect other sequences as well. A058770 is "Numbers n such that n * (1+i)^n 
+ 1 is a Gaussian prime". The current Data section has 1, 2, 3, 5, 9, 19, 
20, 29, 30, 68, 142, 143, 150, 159, 160, 198, 468, 782, 858, 1100, 1137, 
3337, 3638, 3909, 4845, 16895, 21768, 30349, 42692, 48470, 65208, but I 
have removed 160, 21768, and 65208. There is Mathematica code:

 Do[ If[ PrimeQ[ n * (1 + I)^n + 1], Print[n] ], {n, 1, 4000} ]

Although I'm not a Mathematica expert, I think the problem is this: PrimeQ, 
when given an argument that is an ordinary (i.e. rational) integer, will 
test whether that argument is a (rational) prime; when given a Gaussian 
integer that is not an ordinary integer, it will test for a Gaussian prime. 
The trouble here is that for the three values n = 160, 21768, and 65208, n 
* (1+i)^n + 1 is a (rational) prime == 1 (mod 4), and therefore by Fermat 
is not a Gaussian prime. The point is that PrimeQ does not see the I in its 
input, only the result of expanding out the (1 + I)^n, which is a rational 
integer when n is even, and so it only tests for a rational prime. I think 
the fix would be to insert the option GaussianIntegers -> true in the call 
to PrimeQ.

Cheers,
Robert