Mersenne Primes and Lucas-Sequences

created 1998–Nov–28, updated 2016–Sep–20

Mersenne Primes are prime numbers of the form 2^{n}−1. For such numbers a necessary condition to be prime is: the exponent n must be a prime number, also. Therefore, Mersenne primes 2^{p_{k}} − 1 are indexed by p_{k}, respectively k.

k  p_{k} best fit error
1  2  4 0.5337
2  3  5 0.5509
3  5  8 0.6318
4  7  12 0.6087
5  13  17 0.7778
6  17  24 0.6999
7  19  35 0.5383
8  31  51 0.6043
9  61  75 0.8183
10  89  108 0.8215
11 107  157 0.6796
12 127  229 0.5551
13 521  333 1.5669
14 607  483 1.2562
15 1279  702 1.8213
16 2203  1021 2.1587
17 2281  1483 1.5380
18 3217  2155 1.4926
19 4253  3132 1.3578
20 4423  4552 0.9717
21 9689  6615 1.4647
22 9941  9613 1.0341
23 11213 13971 0.8026
24 19937 20303 0.9820
25 21701 29505 0.7355

k  p_{k}  best fit error
26  23209  42878 0.5413
27  44497  62313 0.7141
28  86243  90557 0.9524
29 110503  131602 0.8397
30 132049  191252 0.6904
31 216091  277937 0.7775
32 756839  403914 1.8738
33 859433  586990 1.4641
34 1257787  853046 1.4745
35 1398269 1239693 1.1279
36 2976221 1801589 1.6520
37 3021377 2618168 1.1540
38 6972593 3804866 1.8325
39 13466917 5529439 2.4355
40 20996011 8035684 2.6128
41 24036583 11677897 2.0583
42 25964951 16970960 1.5300
43 30402457 24663130 1.2327
44 32582657 35841813 0.9091
45 37156667 52087288 0.7134
46 42643801 75696103 0.5634
47 43112609 110005727 0.3919
48 57885161 159866353 0.3621
49 74207281 232326550 0.3194
50  … 337629680

Table 1 Exponents of known Mersenne Primes and their approximation.

Because the prime density around x is \frac{1}{ln x} we expect for a random chosen number 2^{n} − 1 to be prime with probability ≈ \frac{2}{n ln 2}. If these numbers were independently of each other we could sum up the probabilities and get the number of Mersenne primes up to a given exponent n to be \frac{2}{ln 2} ∑_{i=1}^{n} \frac{1}{i} ∼ 2 log_{2} n. But we must consider two effects which cause to replace the factor 2 by e^{γ} with γ = 0.57721566 …, the Euler-Mascheroni constant. Firstly the above mentioned restriction, a number 2^{n} −1 can only be prime if n itself is a prime number, and secondly, since Fermat we know that all factors of those numbers have the form 2jn + 1 with certain integers j. Thus, these finer probabilistic arguments assert an approximation formula for the exponent p_{k} of the form C2^{\frac{k}{e^{γ}}}. Making a best-quadratic fitting curve for the 45 values p_{5}, p_{6}, …, p_{49} — these p_{k} which equal 2, 3, 5 or 7 must be prime numbers because already (2p_{k} + 1)^{2} ≥ 2^{p_{k}} − 1 for k ≤ 4 — gives an approximation of ≈ 2.5784·1.4533^{k}, see table 1 and figure 1. Remark: 2^{\frac{1}{e^{γ}}} = 1.47576 . . .. The quotient of p_{k} to this estimation value is headed error in table 1. The correlation coefficient of this logarithmic fitting yields 0.9946 … ≈ 1 − \frac{1}{185}.

[figure: Mersenne-Primes; p_{k}; k]

Figure 1 Exponents of known Mersenne Primes with best fitting line.

According to George Woltmann’s GIMPS of 2016, Sep. 20, the next Mersenne Prime p_{50} is > 68 242 000. For current state of the art, see the URL

http://www.mersenne.org/status.htm. A highly efficient test for checking

numbers of the form 2^{n} − 1 is the Lucus-Lehmer test: define s(1) = 4 and ∀ n ∈ ℕ : s(n + 1) = s(n)^{2} − 2. Then for each odd number q, we have: 2^{q} − 1 is prime, if and only if 2^{q} − 1 divides s(q − 1).

Boolean Two_power_q_minus_1_is_Prime(integer q)
Begin
   integer  i;
   long integer  s;
   s := 4;
   for i from q-2 downto 1 do
      s := s*s-2 mod 2^q-1;
   done
   return  s = 0 ;  # comparision is either true or false
End

Because computation mod 2^{q} − 1 is extremly easy on today’s binary computers, the far mostly expensive operation of the test is the squaring of s. Nevertheless, at first one checks the given number 2^{q} − 1 for small prime divisors of the form 2jq + 1.

To get more examples, let us consider for any P, Q ∈ ℤ \ {0} the socalled generalized Lucas-Sequences, U_{n+1}(P, Q) = P · U_{n} − Q · U_{n−1} with U_{0} = 0 and U_{1} = 1, which are linear second-order recurrence relations. If D := P^{2} −4Q ≠ 0, their general solution is U_{n}(P, Q) = \frac{a^{n}−b^{n}}{a−b} with a = \frac{P +\sqrt{D}}{2} and b = \frac{P −\sqrt{D}}{2}, respectively a + b = P and ab = Q. Thus the Mersenne Numbers correspond to a = 2 and b = 1, respectively P = 3 and Q = 2. When looking for primes in such exponential increasing sequences one likes max{|a|, |b|} to be as small as possible to have a good chance to get many primes. Excluding the case D = 0, the smallest possibility is P = 1 and Q = −1, which defines the well known Fibonacci-sequence. Their numbers F_{n} has a = \frac{1+\sqrt{5}}{2} and b = \frac{1−\sqrt{5}}{2}. Except for n = 4, F_{n} can only be prime if n is a prime number, too. Furthermore each factor of F_{p} has the form 2jp ± 1 for a certain integer j. Thus for primes p ≤ 13 we have (2p−1)^{2} > F_{p} and so F_{p} must be also a prime number. It follows a table of the Fibonacci Primes with n < 2000000.

k n_{k} best fit error
1 3  4 0.8138
2 4  5 0.8037
3 5  7 0.7441
4 7  9 0.7716
5 11  12 0.8981
6 13  17 0.7861
7 17  22 0.7614
8 23  30 0.7630
9 29  41 0.7126
10 43  55 0.7826
11 47  74 0.6335
12 83  100 0.8287

k n_{k} best fit error
13 131  135 0.9687
14 137  183 0.7504
15 359  246 1.4564
16 431  333 1.2951
17 433  449 0.9637
18 449  607 0.7401
19 509  819 0.6215
20 569  1106 0.5146
21 571  1493 0.3825
22 2971  2016 1.4739
23 4723  2721 1.7355
24 5387  3674 1.4661

k  n_{k} best fit error
25 9311  4961 1.8769
26 9677  6698 1.4448
27 14431  9043 1.5959
28 25561  12209 2.0937
29 30757  16483 1.8660
30 35999  22254 1.6176
31 37511  30046 1.2484
32 50833  40566 1.2531
33 81839  54769 1.4943
34 104911* 73945 1.4188
35 130021* 99835 1.3024
36 148091* 134790 1.0987

k  n_{k} best fit error
37 201107* 181983 1.1051
38 397379* 245699 1.6173
39 433781* 331725 1.3077
40 590041* 447869 1.3174
41 593689* 604679 0.9818
42 604711* 816391 0.7407
43 931517* 1102229 0.8451
44 1049897* 1488145 0.7055
45 1285607* 2009180 0.6399
46 1636007* 2712641 0.6031
47 1803059* 3662401 0.4923
48 1968721* 4944694 0.3981
49 2904353*

Table 2 Exponents of known Fibonacci Primes and their approximation. The ∗ indicates an index of a probable-prime.

Similar like in the Mersenne case the best exponentiel approximation of n_{k} for 7 ≤ k ≤ 48 is given by n_{k} ≈ 2^{0.4331k+1.4490}. Here the correlation coefficient is ≈ 0.9929.

[figure: Fibonacci-Primes; p_{k}; k]

Figure 2 Exponents of known Fibonacci Primes with best fitting line.

Checking whether F_{p} is composite is done by searching for an integer g > 1 contradicting g^{F_{p}} ≡ g (mod F_{p}). If F_{p} is composite there are almost no witnesses g satisfying Fermat’s little Theorem. The smallest g for composite F_{p} are 112 if p=19, 48341 if p=31, 181332 if p=37 and 22982 if p=41.

Due to the fact that \sqrt{D} = \sqrt{5} in the “Fibonacci-case” is no integer, there is no simple base to represent these numbers efficiently. Therefore we look at the next smallest possibilities for m; this means m = 2. But apart from P = 3 and Q = 2 there is another option, namely P = 1 Q = −2 This leads to a = 2 and b = −1 and the numbers T_{n} = \frac{2^{n}−(−1)^{n}}{3} . Again, it is necessary that n must be 4 or a prime, for T_{n} to be prime. But because m is a power of 2 — existing digital computers use almost all binary representation of integers — this is highly welcome for practical computation and testing of such numbers for primability. E.g. the checking of Fermats little theorem for the considered T_{n} and a certain a ∈ ℕ — a^{T_{n}−1} ≡ 1 (mod T_{n}) for a ∈ ℕ — is greatly simplified because we can do the calculation modulo 2^{n} + 1 (n of cource is odd in our cases), which is extremly quick on binary computers. The results for n_{k} < 4031400, again with an geometric interpolation curve:

k n_{k} best fit error
1  3  1 2.4340
2  4  2 2.2999
3  5  2 2.0372
4  7  3 2.0212
5 11  5 2.2508
6 13  7 1.8850
7 17  10 1.7468
8 19  14 1.3835
9 23  19 1.1869
10 31  27 1.1336
11 43  39 1.1143
12 61  54 1.1202
13 79  77 1.0281
14 101  108 0.9314
15 127  153 0.8300
16 167  216 0.7734
17 191  305 0.6269
18 199  430 0.4628
19 313  607 0.5159
20 347  856 0.4053
21 701  1208 0.5802
22 1709  1705 1.0024

k  n_{k} best fit error
23  2617  2406 1.0878
24  3539  3395 1.0424
25  5807  4791 1.2121
26  10501  6760 1.5533
27  10691  9540 1.1207
28  11279  13462 0.8379
29  12391  18996 0.6523
30  14479  26806 0.5401
31  42737  37827 1.1298
32  83339* 53379 1.5613
33  95369* 75324 1.2661
34 117239* 106292 1.1030
35 127031* 149992 0.8469
36 138937* 211658 0.6564
37 141079* 298677 0.4723
38 267017* 421471 0.6335
39 269987* 594751 0.4539
40 374321* 839270 0.4460
41 986191* 1184319 0.8327
42 4031399* 1671228 2.4122
43 13347311* 2358320 5.6597
44 13372531* 3327895 4.0183

Table 3 Exponents of known T_{n} Primes and their approximation. The ∗ indicates an index of a probable-prime.

Here the best exponentiel approximation of n_{k} for 5 ≤ k ≤ 44 is given by n_{k} ≈ 2^{0.4969k−0.1952} yielding a correlation coefficient of ≈ 0.98976.

[figure: \frac{2^{p}+1}{3}-Primes; p_{k}; k]

Figure 3 Exponents of known \frac{2^{p}+1}{3}-Primes with best fitting line.

   At the URL http://www.primenumbers.net/prptop/ the current search

status can be seen. As in the ’Mersenne case’ all factors of T_{n} Primes have the form 2jn + 1 with certain integers j. Thus T_{p} must be prime if p < 11 and it is an odd prime. Further Henri Lifchitz deduced from Fermats little Theorem a very quick probable prime test for a given T_{p} which checks the necessary condition 25^{2^{p−1}} ≡ 25 (mod 2^{p} +1). Similar to the Lucas-Lehmer test for Mersenne-Primes, we get for T_{q} a very fast probabilistic test for odd q of the form:

Boolean Two_power_q_plus_1_divided_by_3_is_Prime(integer q)
Begin
   integer  i;
   long integer  s;
   s := 25;

for i from q-1 downto 1 do
s := s*s mod 2^q+1;
done
return  s = 25;  # comparision is either true or false
End

Again, because computation mod 2^{q} + 1 is extremly easy on today’s binary computers, the far mostly expensive is the squaring which can be done via Fast-Fourier-Transformations in time O(q ln q). Moreover, it is conjectured that this test is even sufficient, meaning any number q passing the test, must correspond to a prime number T_{q}.

A few remarks to the Lucas-sequence U_{n}(P, Q) if D ≤ 0. In the case D = 0 we get the solution U_{n} = n(\frac{P}{2} )^{n−1} which for P ≠ 2 is never a prime number. If D < 0 these numbers U_{n} = \begin{matrix}\sqrt{Q} & ^{n} \\ √ & \overline{−D} \end{matrix} sin(nϕ) rotate by the angle

ϕ = arctan \frac{\sqrt{−D}}{P} with increasing n. In this case the sequence U_{n} grows asymptotically like ≈ \sqrt{Q}^{n}. Thus, the cases P = 1, Q = 2 and P ∈ \{1, 2\}, Q = 3 are most interesting — in the case P ∈ {1, 3}, Q = 4 only finite many prime numbers occur because U_{n}(P, 4) is reducible.

Now we consider firstly the sequence with slowest growth, so P = 1, Q = 2 and D = −7 given U_{n}(1, 2) = \frac{(1+\sqrt{7}i)^{n}−(1−\sqrt{7}i)^{n}}{2^{n}\sqrt{7}i}. The prime numbers and their approximation for n_{k} < 121540 are:

k n_{k} best fit error
1 4+  3 1.2814
2 6+  4 1.6063
3 7  4 1.5661
4 8+  5 1.4957
5 9+  6 1.4062
6 10+  8 1.3057
7 11  9 1.2003
8 15+  11 1.3678
9 17  13 1.2955
10 19  16 1.2100
11 23  19 1.2240
12 25+  22 1.1119
13 26+  27 0.9663
14 29  32 0.9007
15 31  39 0.8046
16 47  46 1.0195
17 53  55 0.9607
18 65+  66 0.9847
19 67  79 0.8482
20 71  95 0.7511
21 73  113 0.6454
22 113  135 0.8349

k n_{k} best fit error
23 127  162 0.7841
24 199*  194 1.0268
25 257*  232 1.1082
26 349*  278 1.2576
27 421*  332 1.2678
28 433*  397 1.0897
29 449*  475 0.9443
30 691*  569 1.2144
31 761*  681 1.1177
32 823*  815 1.0102
33 991*  975 1.0165
34 1237*  1167 1.0603
35 1277*  1396 0.9148
36 1399*  1670 0.8375
37 1531*  1999 0.7659
38 1571*  2392 0.6568
39 3461*  2862 1.2092
40 3697*  3425 1.0794
41 4933*  4098 1.2036
42 6199*  4904 1.2640
43 7351*  5869 1.2526
44 9551*  7022 1.3601

k  n_{k} best fit error
45 9719*  8403 1.1566
46 11681* 10055 1.1617
47 12037* 12032 1.0004
48 14629* 14398 1.0160
49 14951* 17229 0.8678
50 19079* 20617 0.9254
51 20327* 24670 0.8239
52 22549* 29521 0.7638

k  n_{k} best fit error
53 30517* 35325 0.8639
54 51511* 42271 1.2186
55 52813* 50582 1.0441
56 60923* 60528 1.0065
57 73943* 72429 1.0209
58 79687* 86670 0.9194
59 91249*
60 117017*

Table 4 Exponents of known U_{n}(1, 2) Primes and their approximation. The + indicates an index which is not prime. The ∗ indicates an index of a probable-prime.

Because U_{2} =1 and U_{3} =U_{5} =U_{13} = −1 are the only units in the sequence, U_{n_{k}} is a prime number for the above non-prime indices n_{k} marked +.

[figure: U_{n}(1, 2)-Primes; p_{k}; k]

Figure 4 Exponents of known U_{n}(1, 2)-primes with best fitting line.

Here the best exponentiel approximation of n_{k} for 8 < k ≤ 59 is given by n_{k} ≈ 2^{0.2585k+1.3939} yielding a correlation coefficient of ≈ 0.9977.

For all P, Q ∈ ℤ \{0}, the general linear 2nd order difference equation, W_{n+1}(P, Q) = P ·W_{n} −Q·W_{n−1} with given W_{0}, W_{1} ∈ ℤ, has two independent solutions. In the case D ≠ 0 these may be given by U_{n}(P, Q) = \frac{a^{n}−b^{n}}{a−b} and by V_{n}(P, Q) = a^{n} + b^{n} with a = \frac{P +\sqrt{D}}{2} and b = \frac{P −\sqrt{D}}{2}. Especially holds V_{0} = 2 and V_{1} = U_{2} = P. Then any solution W_{n} can be writen as W_{n} = αU_{n} + βV_{n} with α = W_{1} − P \frac{W_{0}}{2} and β = \frac{W_{0}}{2} . Further, for all positive integers m and j holds U_{j} divides U_{mj} and if j is odd, then also V_{j} divides V_{mj}.

Next for those P and Q with gcd(P, Q) = 1 which 2nd order sequences U_{n}(P, Q) grows asymptotically slowest, a table follows indicating the number of primes and units in its corresponding sequence:

\begin{matrix}P & Q & D & # & r & n \\ 1 & 0 & 1 & 0 & * & ℕ \\ 1 & 1 & -3 & 0 & * & ℕ \\ 2 & 1 & 0 & ∞ & 1 \\ 1 & 2 & -7 & 33+7 & 1,2,3,5,13 \\ 1 & -1 & 5 & 21+1 & 1,2 \\ 2 & 3 & -8 & 17+6 & 1,3 \\ 1 & 3 & -11 & 18+6 & 1,2,5 \\ 1 & -2 & 9 & 21+3 & 1,2 \\ 3 & 2 & 1 & 14+4 & 1 \\ 3 & 4 & -7 & 5 & (-1) & 1 \\ 1 & 4 & -15 & 4 & (-3) & 1,2 \\ 4 & 5 & -4 & 15+3 & 1 \\ 3 & 5 & -11 & 9+1 & 1 \\ 2 & 5 & -16 & 24+5 & 1,3 \\ 1 & 5 & -19 & 15+2 & 1,2,7 \\ 1 & -3 & 13 & 10+3 & 1,2 \\ 2 & -1 & 8 & 17+6 & 1 \\ 1 & 6 & -23 & 13+3 & 1,2 \\ 3 & 1 & 5 & 1 & (1) & 1 \\ 5 & 7 & -3 & 11+5 & 1 \\ 4 & 7 & -12 & 12+4 & 1 \\ 3 & 7 & -19 & 17+3 & 1 \\ 2 & 7 & -24 & 12+2 & 1 \\ 1 & 7 & -27 & 11+4 & 1,2 \\ 5 & 8 & -7 & 3 & * & 1 \\ 3 & 8 & -23 & 14+3 & 1,3 \\ 1 & 8 & -31 & 13+0 & 1,2 \\ 2 & -3 & 16 & 11+4 & 1 \\ 4 & 3 & 4 & 6+3 & 1 \\ 5 & 6 & 1 & 11+6 & 1 \\ 5 & 9 & -11 & 2 & (-1) & 1 \end{matrix}  \begin{matrix}P & Q & D & # & r & n \\ 4 & 9 & -20 & 1 & (-2) & 1 \\ 2 & 9 & -32 & 3 & (-4) & 1 \\ 1 & 9 & -35 & 1 & (-5) & 1,2 \\ 3 & 10 & -31 & 11+1 & 1,3 \\ 1 & 10 & -39 & 11+4 & 1,2 \\ 3 & -1 & 13 & 12+0 & 1 \\ 6 & 11 & -8 & 10+2 & 1 \\ 5 & 11 & -19 & 10+1 & 1 \\ 4 & 11 & -28 & 11+3 & 1 \\ 3 & 11 & -35 & 6+1 & 1 \\ 2 & 11 & -40 & 14+2 & 1 \\ 1 & 11 & -43 & 5+3 & 1,2 \\ 5 & 12 & -23 & 14+2 & 1 \\ 1 & 12 & -47 & 13+3 & 1,2 \\ 7 & 13 & -3 & 14+1 & 1 \\ 6 & 13 & -16 & 8+1 & 1 \\ 5 & 13 & -27 & 11+2 & 1 \\ 4 & 13 & -36 & 10+1 & 1 \\ 3 & 13 & -43 & 10+1 & 1 \\ 2 & 13 & -48 & 10+2 & 1 \\ 1 & 13 & -51 & 7+2 & 1,2 \\ 4 & 1 & 12 & 0 & (2) & 1 \\ 5 & 14 & -31 & 6+3 & 1 \\ 3 & 14 & -47 & 8+3 & 1 \\ 1 & 14 & -55 & 8+1 & 1,2 \\ 7 & 15 & -11 & 12+3 & 1 \\ 4 & 15 & -44 & 10+2 & 1,3 \\ 2 & 15 & -56 & 8+0 & 1 \\ 1 & 15 & -59 & 7+5 & 1,2 \\ 7 & 12 & 1 & 12+3 & 1 \\ 5 & 4 & 9 & 1 & (1) & 1 \end{matrix}

\begin{matrix}P & Q & D & # & r & n \\ 7 & 16 & -15 & 1 & (-1) & 1 \\ 5 & 16 & -39 & 2 & (-3) & 1 \end{matrix}  \begin{matrix}P & Q & D & # & r & n \\ 3 & 16 & -55 & 2 & (-5) & 1 \\ 1 & 16 & -63 & 1 & (1) & 1,2 \end{matrix}

Table 5 # denotes the number of primes up to 1024 (plus these up to 4096) and the index n these U_{n} with U_{n} = ±1.

A sequence U_{n}(P, Q) is called reducible if ∃ r, P′, Q′ ∈ ℤ and k ∈ ℕ\{1} such that holds r^{n}U_{n}(P, Q) = U_{kn}(P′, Q′). Reducible U_{n}(P, Q) are for example U_{n}(1, 16) = \frac{1}{3·2^{n}} U_{2n}(3, 4) and U_{n}(4, 1) = 2^{−n}U_{2n}(2, −2) and U_{n}(5, 8) = (−1)^{n}U_{3n}(1, 2).

Next, we deduce reducibility of given P and Q in the case k = 2 explicitly. Notice with r = P ± 2\sqrt{Q} = (\sqrt{a} ± \sqrt{b})^{2} holds

r\frac{P ± \sqrt{D}}{2} = (\frac{r ± \sqrt{D}}{2})^{2}

Thus if \sqrt{Q} ∈ ℕ and r ≠ 0, which means D ≠ 0, we may transfrom

U_{n}(P, Q) = \frac{1}{r^{n}} U_{2n}(r, ±r\sqrt{Q}) = r^{n−1}U_{2n}(1, \frac{±\sqrt{Q}}{r}) with r = P ± 2\sqrt{Q}

and we get \frac{U_{n}(r,±r\sqrt{Q})}{r^{⌊ \frac{n}{2} ⌋}} divides U_{n}(P, Q). Only if this divisor equals ±1, which can be only finite many times, U_{n}(P, Q) can be a prime number.

Generally, if D < 0, a sufficient, but not necessary condition for the irreducibility of U_{n}(P, Q) is ^{k} \sqrt{Q} ∉ ℕ for all k ≥ 2. For k = 3 a necessary and sufficient condition is the existence of an integer root z := P′−P of the cubic equation z^{3}−12Qz−8P Q = 0. If D > 0 this occurs iff \frac{z}{2} = ^{3} \sqrt{Q}( ^{3} \sqrt{a}+^{3} \sqrt{b}) ∈ ℤ holds. Furthermore always Q′ = ^{3} \sqrt{Q}.

For which n is U_{n}(P, Q) = ±1? Per definition, always if n = 1. Because U_{2} = P if and only if P = ±1 then n = 2 is a solution. Next we have U_{3} = P^{2} − Q. Thus n = 3 iff P^{2} = Q ± 1. The case n = 4 iff P^{3} − 2P Q = ±1 and , because U_{2} divides U_{4} = ±1, then U_{2} = P must be also ±1. This means 1 − 2Q = ±1, leaving only the possibilities Q = 0 or Q = 1. Finally n = 5

is a root iff P^{4} − 3P^{2}Q + Q^{2} = ±1 which is equivalent to P^{2} = ^{\frac{3Q±\sqrt{5Q^{2}±4}}{2}}. So the Pell-equation 5Q^{2} ± 4 = x^{2} for x ∈ ℤ must hold. Its general solution has the form Q = F_{j} which F_{j} a Fibonacci number and x the companion Fibonacci number resulting in the condition P^{2} = F_{j±2}. An overview of all solutions n of U_{n}(P, Q) = ±1 follows:

n  necessary and sufficient condition
∗  P = 1, Q ∈ {0, 1}
1  1 = 1
2  P = 1
3  Q = P^{2} ± 1
5 P = 1, Q ∈ {0, 1, 2, 3} or P = 12, Q ∈ {55, 377}
7  P = 1, Q ∈ {0, 1, 5}
13  P = 1, Q ∈ {0, 1, 2}

Conjecture: Let m ∈ ℤ divide U_{p}(P, Q) with p a prime number. If p doesn’t divide D = P^{2} − 4Q then m = 2jp ± 1 with a certain j ∈ ℤ and if even \sqrt{D} ∈ ℕ then m = 2jp + 1. If p divides D then either p divides m, too or m = 2jp ± 1.

Proof: We have

U_{n}(P, Q) = \frac{2}{2^{n}} ∑_{i=0}^{⌊ \frac{n−1}{2} ⌋} \binom{n}{2i + 1}P^{n−2i−1}D^{i}

So if n is an odd prime number p, then all summands, except the very last, are divisible by p because \binom{p}{2i+1} is divisible by p and it holds

U_{p} ≡ D^{\frac{p−1}{2}} (mod p)

Clearly, if p divides D, then it divides U_{p}, too. In the case gcd(D, p) = 1 we know by Fermat’s little theorem that D^{p−1} ≡ 1 (mod p). Thus D^{\frac{p−1}{2}} ≡ ±1 (mod p) and if \sqrt{D} ∈ ℕ then \sqrt{D}^{p−1} ≡ 1 (mod p).

Have a look for n < 10000 at the primes in three selected sequences:

\begin{matrix}U_{n}(2, 5) = \\ k & n_{k} \\ 1 & 2 \\ 2 & 5 \\ 3 & 7 \\ 4 & 9 & + \\ 5 & 11 \\ 6 & 13 \\ 7 & 17 \\ 8 & 19 \\ 9 & 31 \\ 10 & 41 \\ 11 & 67 \\ 12 & 71 \\ 13 & 79 \\ 14 & 89 \\ 15 & 107 \end{matrix} \frac{(1+2i)^{n}−(1−2i)^{n}}{4i} ; k n_{k} ; 16 137 ; 17 149 ; 18 157 ; 19 197 ; 20 241 ; 21 311 ; 22 383 ; 23 479 ; 24 587 ; 25 1151 ; 26 1237 ; 27 1693 ; 28 3187 ; 29 3527 ; 30 5323

\begin{matrix}U_{n}(2, 3) = \\ k & n_{k} \\ 1 & 2 \\ 2 & 5 \\ 3 & 7 \\ 4 & 9 & + \\ 5 & 11 \\ 6 & 13 \\ 7 & 17 \\ 8 & 23 \\ 9 & 43 \\ 10 & 59 \\ 11 & 61 \\ 12 & 67 \\ 13 & 103 \end{matrix} \frac{(1+\sqrt{2}i)^{n}−(1−\sqrt{2}i)^{n}}{\sqrt{8}i} ; k n_{k} ; 14 389 ; 15 443 ; 16 449 ; 17 617 ; 18 683 ; 19 1061 ; 20 1471 ; 21 1709 ; 22 2297 ; 23 2473 ; 24 3181 ; 25 5801 ; 26 9883  \begin{matrix}U_{n}(1, 3) = \\ k & n_{k} \\ 1 & 3 \\ 2 & 4 & + \\ 3 & 7 \\ 4 & 10 & + \\ 5 & 13 \\ 6 & 17 \\ 7 & 19 \\ 8 & 25 & + \\ 9 & 47 \\ 10 & 61 \\ 11 & 79 \\ 12 & 127 \\ 13 & 131 \end{matrix} \frac{(1+\sqrt{11}i)^{n}−(1−\sqrt{11}i)^{n}}{2^{n}\sqrt{11}i} ; k n_{k} ; 14 179 ; 15 269 ; 16 283 ; 17 311 ; 18 881 ; 19 1279 ; 20 1489 ; 21 1747 ; 22 1879 ; 23 2131 ; 24 3361 ; 25 4451 ; 26 5387 ; 27 8563

12