AGAPProgram

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;