## On the Density of Some Exceptional Primes

```As discussed in the note On Case 1 of Fermat's Last Theorem, the
congruence

(x+1)^p - (x)^p  -  1  =  0    (mod p^2)             (1)

must have a non-trivial solution x coprime to p in order for there to
exist an integer solution of Fermat's equation x^p + y^p + z^p = 0
with xyz coprime to p.  (Incidentally, the numbers 1093 and 3511 are
the only known primes such that x=1 is a root of congruence (1)).
However, this congruence has at least one non-trivial root for every
prime p of the form 3k+1.  (Here non-trivial means other than x = 0
or -1 modulo p).  On the other hand, there are no non-trivial
solutions at all for most primes of the form 3k-1.  The exceptions
less than 7000 are listed below

59   83  179  227  419  443  701  857  887  911  929  971  977
1091 1109 1193 1217 1223 1259 1283 1289 1439 1487 1493 1613 1637
1811 1847 1901 1997 2003 2087 2243 2423 2477 2579 2591 2729 2777
2969 3089 3137 3191 3203 3251 3467 3527 3533 3881 3917 3929 4001
4079 4091 4643 4649 4691 4889 4943 5189 5303 5351 5711 5843 5849
5903 6257 6359 6389 6551 6569 6581 6947 6959

The total numbers of primes of the form 3k-1, and the numbers of
these "exceptional" primes, up to various levels are listed below

Tot  Excpt  dens

1000   86   13   0.151
2000  153   30   0.196
3000  221   40   0.181
4000  277   51   0.184
5000  337   59   0.175
6000  397   66   0.166
7000  454   74   0.162

Clearly if x is a root of (1) for a given prime p, then so is x+kp
for every k, so we can represent all the non-trivial roots by just
those in the range 1 to p-2.  Furthermore, if x is a root, then
-(x+1) is also a root, in view of

(-(x+1)+1)^p - (-(x+1))^p  =  (-x)^p - (-(x+1))^p  (mod p^2)

Therefore, we can represent all the non-trivial roots by the set of
roots in the range from 1 to (p-1)/2.

The table below shows the number of primes of the form 3k-1
less than 160000, along with the number exceptional primes in
this range.

# of primes     number of
range           the form      exceptional
3k-1          primes         ratio
--------------   -----------    -----------    ---------
0 -  10000       616            96          0.15584
10000 -  20000       520            83          0.15961
20000 -  30000       497            83          0.16700
30000 -  40000       479            72          0.15031
40000 -  50000       463            66          0.14254
50000 -  60000       466            77          0.16523
60000 -  70000       449            74          0.16481
70000 -  80000       448            65          0.14508
80000 -  90000       436            71          0.16284
90000 - 100000       432            78          0.18055
100000 - 110000       427            66          0.15456
110000 - 120000       432            67          0.15509
120000 - 130000       440            73          0.16590
130000 - 140000       427            63          0.14754
140000 - 150000       408            64          0.15686
150000 - 160000       421            59          0.14014
-----         -----         --------
0 - 160000      7361          1157          0.15717

It certainly appears that the density of exceptional primes among
all primes of the form 3k-1 is fairly constant, independent of the
size of the primes involved.  It's also worth noting that whenever
a prime of this form has (non-trivial) solutions, the number of
solutions in the reduced range 1 to (p-1)/2 is a multiple of
three (as explained below).  The distribution of the number of
solutions for the primes less than 160000 is shown below

# of primes     # of solutions
p = 3k-1        solutions
-----------     --------------
6204                0
1055                3
98                6
4                9

The density seems to drop roughly by factors of 1/6, 1/11, and 1/24
for each additional set of three roots, suggesting that primes with
12 roots would have a density of roughly 1/48 times the density of
primes with 9 roots.

The reason that these roots occur in sets of 3 is easier to see if we
consider the whole range of residues from 1 to p-2.  Over this range
the solutions of

(x+1)^p - x^p - 1  =  0   (mod p^2)

occur in sets of six related values, because if x is a solution, then
another solution is given by replacing x with either -(1+x) as we saw
previously, or with -(1 + 1/x), where the inverse is taken modulo p.
The former has a period of 2, whereas the latter has a period of 3,
and together these substitutions generate a group of order 6 with
the members

a(x) = x             b(x) = 1/x            c(x) =  -(1+x)

d(x) = -1/(1+x)      e(x) = -x/(1+x)       f(x) = -(1+x)/x

This is isomorphic to the group of symmetries of an equilateral
triangle.  The operation with period 2 is like flipping the triangle
over in the plane, and the operation with period 3 is like rotating
the triangle through 120 degrees.  In general, these operations can
be used to partition the residues modulo p into equivalence classes
of 6 elements that are related according to these operations.  For
example, with p=59 we get the following sets

a   b   c   d   e   f
--  --  --  --  --  --
1   1  57  29  29  57
2  30  56  39  19  28
3  20  55  44  14  38  *
4  15  54  47  11  43  *
5  12  53  49   9  46
6  10  52  42  16  48
7  17  51  22  36  41
8  37  50  13  45  21
18  23  40  31  27  35
24  32  34  33  25  26

The two sets marked with asterisks are the values that satisfy (1)
in the range from 1 to p-2.  Obviously for each set we have the
relations
a + c  =  d + e  =  b + f  =  -1
(mod p)
ab  =  cd  =  ef  =  1

Notice that for primes congruent to +1 (mod 3) there is always a
degenerate set corresponding to the algebraic factor (1 + x + x^2)^2
of equation (1) when p is of this form.  (For more on this, see the
note Sums of Powers in Terms of Symmetric Functions.  For example,
with p=31 we have the sets

a   b   c   d   e   f
--  --  --  --  --  --
1   1  29  15  15  29
2  16  28  10  20  14
3  21  27  23   7   9
4   8  26   6  24  22
5  25  25   5  25   5
11  17  19  18  12  13

Here we see that 5 and 25 are roots of (1+x+x^2)^2 = 0 (mod 31)^2.
Also, the set containing 1 always consists of just three distinct
elements, +1, -2, and (p-1)/2.  So, if p is of the form 6k+1, the
residues from 1 to p-2 are partitioned into one set with 3 members,
one set with 2 members, and then the remaining

(p-2)-3-2  =  6k+1-2-3-2 = 6(k-1)

residues are partitioned into k-1 sets.  On the other hand, if p
is of the form 6k-1, the algebraic factor (1+x+x^2) occurs only
to the 1st power, so there is no degenerate set of 2 elements.
Thus the residues are patitioned into one set of 3 elements, and
then the remaining

(p-2)-3  =  6k-1-2-3  =  6(k-1)

residues are partitioned into k-1 sets.  We recall that only the
primes 1093 and 3511 have the degenerate set containing 1 as solutions
of (1).

Just to give one more example of an exceptional prime, with p=83 we
have the partition

a   b   c   d   e   f
--  --  --  --  --  --
1   1  81  41  41  81
2  42  80  55  27  40
3  28  79  62  20  54
4  21  78  33  49  61
5  50  77  69  13  32
6  14  76  71  11  68
7  12  75  31  51  70
8  52  74  46  36  30  *
9  37  73  58  24  45
10  25  72  15  67  57
16  26  66  39  43  56
17  44  65  23  59  38
18  60  64  48  34  22
19  35  63  29  53  47

Again the set marked with an asterisk represents the solutions of
congruence (1) in the range 1 to p-2.

Regarding the density of these primes, notice that the expression
(x+1)^p - x^p - 1  is always a multiple of p, so we can divide that
out, and we are left with a polynomial that wish to set equal to
zero modulo p (so that the full expression vanishes modulo p^2).
For any given value of x, if we assume the value of this polynomial
is randomly distributed over the residues mod p, we can say that
the chances are about 1/p that it will vanish.

Now, the number of independent x values that we can try is really
only (p-5)/6  (which is k-1 from above), because the solutions
cover the equivalence classes described above.  Thus, for any given
prime p we have (p-5)/6 attempts, each of which has a 1/p probability
of success, so overall the probability of at least one solution
(up to equivalence class) should be about

[(p-5)/6](1/p)  =  (1/6)(1 - 5/p)

Thus it isn't too surprising that the density is roughly 1/6 = 0.167.

It would be interesting to know if there is any other way of
characterizing these exceptional primes, besides evaluating (1).
```