|
@@ -1,8 +1,7 @@
|
1
|
1
|
# Finding lucas pseudoprimes
|
2
|
|
-from math import sqrt, gcd, floor, log
|
3
|
|
-from time import process_time
|
4
|
|
-from psimp import trd, hund_div
|
5
|
|
-from sieve import sieve
|
|
2
|
+from math import gcd, floor, log
|
|
3
|
+from psimp import trd
|
|
4
|
+
|
6
|
5
|
|
7
|
6
|
def isqrt(n):
|
8
|
7
|
""" Find the integer square root of n via newton's method, code via
|
|
@@ -17,6 +16,7 @@ def isqrt(n):
|
17
|
16
|
y = (x + n // x) // 2
|
18
|
17
|
return x
|
19
|
18
|
|
|
19
|
+
|
20
|
20
|
def hasIntSQRT(n):
|
21
|
21
|
"""Detect whether the square root of n is an integer,
|
22
|
22
|
i.e. whether the isqrt(n) is the true square root.
|
|
@@ -24,6 +24,7 @@ def hasIntSQRT(n):
|
24
|
24
|
isq = isqrt(n)
|
25
|
25
|
return isq * isq == n
|
26
|
26
|
|
|
27
|
+
|
27
|
28
|
def Dsequence():
|
28
|
29
|
"""Generate sequence 5, -7, 9, -11, 13, -15...
|
29
|
30
|
"""
|
|
@@ -35,6 +36,7 @@ def Dsequence():
|
35
|
36
|
yield -val
|
36
|
37
|
val = val + 2
|
37
|
38
|
|
|
39
|
+
|
38
|
40
|
def Legendre(a, p):
|
39
|
41
|
"""Function for calculating Legendre symbol.
|
40
|
42
|
|
|
@@ -47,11 +49,12 @@ def Legendre(a, p):
|
47
|
49
|
"""
|
48
|
50
|
if (p % 2 == 0):
|
49
|
51
|
raise ValueError("p must be odd, is {}".format(p))
|
50
|
|
- lv = pow(a, (p-1)//2, p)
|
|
52
|
+ lv = pow(a, (p-1)//2, p)
|
51
|
53
|
if lv == p - 1:
|
52
|
54
|
lv = -1
|
53
|
55
|
return lv
|
54
|
56
|
|
|
57
|
+
|
55
|
58
|
def Jacobi(a, n):
|
56
|
59
|
"""Function for calculating Jacobi symbol.
|
57
|
60
|
|
|
@@ -65,7 +68,7 @@ def Jacobi(a, n):
|
65
|
68
|
if (n % 2 == 0):
|
66
|
69
|
raise ValueError("n must be odd")
|
67
|
70
|
if a % n != a:
|
68
|
|
- return Jacobi(a%n, n)
|
|
71
|
+ return Jacobi(a % n, n)
|
69
|
72
|
if a != 0 and a % 2 == 0:
|
70
|
73
|
nm8 = n % 8
|
71
|
74
|
if (nm8 == 3 or nm8 == 5):
|
|
@@ -82,6 +85,7 @@ def Jacobi(a, n):
|
82
|
85
|
return -1 * Jacobi(n, a)
|
83
|
86
|
return Jacobi(n, a)
|
84
|
87
|
|
|
88
|
+
|
85
|
89
|
def FirstD(n):
|
86
|
90
|
"""Return first D in the sequence 5, -7, 9, -11, 13, -15... for which the
|
87
|
91
|
Jacobi symbol (D/n) is -1. Returns 0 if the number is a perfect square, or
|
|
@@ -107,6 +111,7 @@ def FirstD(n):
|
107
|
111
|
# Shouldn't fire
|
108
|
112
|
return 0
|
109
|
113
|
|
|
114
|
+
|
110
|
115
|
def Lucas(n, p, q):
|
111
|
116
|
"""Function for generating values of a lucas sequence with parameters p, q.
|
112
|
117
|
Via formula at <https://en.wikipedia.org/wiki/Lucas_sequence>
|
|
@@ -122,6 +127,7 @@ def Lucas(n, p, q):
|
122
|
127
|
Vc = Vn
|
123
|
128
|
return Uc, Vc
|
124
|
129
|
|
|
130
|
+
|
125
|
131
|
def LucasUn(n, p, q):
|
126
|
132
|
"""Return only the U numbers from a lucas sequence Un(P, Q). For example if
|
127
|
133
|
P = 1, and Q = -1 this will return the Fibonacci numbers; if P = 3 and Q =
|
|
@@ -129,6 +135,7 @@ def LucasUn(n, p, q):
|
129
|
135
|
"""
|
130
|
136
|
return Lucas(n, p, q)[0]
|
131
|
137
|
|
|
138
|
+
|
132
|
139
|
def LucasVn(n, p, q):
|
133
|
140
|
"""Return only V numbers from a local sequence Vn(P, Q).
|
134
|
141
|
"""
|
|
@@ -169,6 +176,7 @@ def LucasFast(k, p, q, n):
|
169
|
176
|
tp = tp // 2
|
170
|
177
|
return(uc, vc)
|
171
|
178
|
|
|
179
|
+
|
172
|
180
|
def LucasPrime(n):
|
173
|
181
|
"""Lucas probable prime test
|
174
|
182
|
(note: not the "strong" test.)
|
|
@@ -186,6 +194,7 @@ def LucasPrime(n):
|
186
|
194
|
lnum = LucasFast(n + 1, 1, q, n)
|
187
|
195
|
return lnum[0] == 0
|
188
|
196
|
|
|
197
|
+
|
189
|
198
|
def StrongLucasPrime(n):
|
190
|
199
|
"""Strong version of the Lucas probable prime test
|
191
|
200
|
"""
|
|
@@ -237,15 +246,3 @@ def StrongLucasPrime(n):
|
237
|
246
|
return True
|
238
|
247
|
ks = ks * 2
|
239
|
248
|
return False
|
240
|
|
-
|
241
|
|
-
|
242
|
|
-
|
243
|
|
-mp = 400000
|
244
|
|
-ap = sieve(mp)
|
245
|
|
-for i in range(mp):
|
246
|
|
- iip = i in ap
|
247
|
|
- lp = StrongLucasPrime(i)
|
248
|
|
- if iip != lp:
|
249
|
|
- print(iip, lp, i, hund_div(i), sep='\t')
|
250
|
|
-
|
251
|
|
-
|