primes1mod4(n)=
{
p=n;s=listcreate(100);A=(2*n)^2+1;listput(s,n);\
if(n%4==1&&isprime(n)==1,\
 for(k=1,4,v=factor(A);\
  for(i=1,omega(A),P=v[i,1];listput(s,P));\
 p=p*prod(i=1,omega(A),v[i,1]);A=(2*p)^2+1);\
for(j=1,length(s),print(s[j])),\
print("n must be a prime congruent to 1 mod 4"))
}

addhelp(primes1mod4 , "primes1mod4(n) : Generates a list of primes congruent to 1 mod 4 using the procedure described in the book.")

printp("primes1mod4(n) : Generates a list of primes congruent to 1 mod 4 using the procedure described in the book.")