The output of the following GAP code is the code to be used in Pari/GP to find $\mathcal P_k$ for a given $k\geq 3$. Note that one can use the same idea to check if a given given prime $p$ is in $\mathcal P_k$ for a given $k$ without trying to calculate explicitly the set $\mathcal P_k$. The program is here. However, to confirm that $p\not\in\mathcal P_k$ with $k$ large can take a long time since one still needs to go through all possible $a\cdot d-b\cdot c$ (see below). The GAP used is of version 4.4.12 and Pari/GP is of version 2.3.5. Note that "factor(L)" does not seem to work for Pari/GP version 2.4.3. This command is used to factor every number contained in the list L.
## 1. Run
##
## gap < pk.gap > output_filename
##
## so that output_filename will capture the output.
##
## 2. After the program finishes, edit output_filename to remove
## the gap logo, and "gap> > >" (at the beginning of the
## output file, a few lines down, and the end of the file.)
##
## 3. Run
##
## gp < output_filename > pk_filename
##
## Thus, the final list of primes will be held in pk_filename
##
## Define a function in gp to collect the primes after factorization
##
Print( "\n\n",
"\\\\A function to collect primes\n",
"collect( W ) = {\n",
"U = [];\n",
"for( i=1, length(L), U = setunion( U, Vec(W[i][,1]) ) );\n",
"U = eval(U);\n",
"}\n\n"
);
## Tell gap how many loops to compute.
zeta:= 0;
i:= 0;
total_time:= 0;
for k in [4..60] do
Print( "\n\n" );
Print( "\n{\nL = " );
ini_time:= Runtime();;
## zeta is a primitive root of $k$-th unity
## which generate a subgroup $G$ of order $k$ of the multiplicative
## group of the field Q[zeta].
zeta:= E(k);
## C hold the elements $\{zeta^i-1\mid i =1,2,\dots, k-1\}$ to be
## used in computing $|\Norm( (a-1)(c-1)-(b-1)(d-1) )|$
## where $a,b,c,d$ are from the group $G$.
C:= List([1..k-1], i-> zeta^i-zeta^0);
## W will hold the list of integers $|\Norm( (a-1)(c-1)-(b-1)(d-1) )|$
## where $a,b,c,d$ are elements from $G$ with $(a,b)\not=(c,d)$ and
## $(a,d)\not=(b,c)$.
W:= [];
for i in [1..k-2] do
if k mod i = 0 then
a:= C[i];
for j in [i+1..k-1] do
b:= C[j];
for s in [j..k-1] do
c:= C[s];
for t in [i..k-1] do
d:= C[t];
if (i <> j) and (i <> s) and (t <> j) and (t <> s) then
if not(k mod t = 0) or (t >= i) then
n:= AbsoluteValue(Norm( a*d-b*c ));
if n > 1 then
AddSet( W, n );
fi;
fi;
fi;
od;
od;
od;
fi;
od;
## Produce the commands to be used in GP for factorizing
## the integers held in W
Print( W, ";\n}\n" );
Print( "collect( factor( L ) );\n",
"print( \"k = ",
k,
"\" );\n",
"print( vecsort( U ) );\n\n",
"print( \" \" );\n"
);
### Print the time used for computing the list, as comments in GP,
### every time a loop is completed.
end_time:= Runtime();;
Print( "\\\\ Runtime this loop: ", HMSMSec( end_time - ini_time ), "\n\n");
total_time:= total_time + (end_time - ini_time);
ini_time:= end_time;;
od;
Print( "\\\\ Total runtime: ", HMSMSec( total_time ), "\n\n\n" );
quit;