Przeglądaj źródła

Non-recursive Lucas; faster FirstD; time comparisons

Petra Lamborn 5 lat temu
rodzic
commit
62753715e8
1 zmienionych plików z 36 dodań i 12 usunięć
  1. 36
    12
      lucas.py

+ 36
- 12
lucas.py Wyświetl plik

1
 # Finding lucas pseudoprimes
1
 # Finding lucas pseudoprimes
2
 from math import sqrt, gcd, floor, log
2
 from math import sqrt, gcd, floor, log
3
+from time import process_time
3
 
4
 
4
 def isqrt(n):
5
 def isqrt(n):
5
     """ Find the integer square root of n via newton's method, code via
6
     """ Find the integer square root of n via newton's method, code via
83
     """Return first D in the sequence 5, -7, 9, -11, 13, -15... for which the
84
     """Return first D in the sequence 5, -7, 9, -11, 13, -15... for which the
84
     Jacobi symbol (D/n) is -1
85
     Jacobi symbol (D/n) is -1
85
     """
86
     """
86
-    if hasIntSQRT(n):
87
-        raise ValueError("n must not be a square")
88
     if n < 1 or n % 2 == 0:
87
     if n < 1 or n % 2 == 0:
89
         raise ValueError("n must be a positive odd number")
88
         raise ValueError("n must be a positive odd number")
89
+    # Flag to see if it's been checked whether or not n is
90
+    # square
91
+    haschecked = False
90
     for D in Dsequence():
92
     for D in Dsequence():
91
         if Jacobi(D, n) == -1:
93
         if Jacobi(D, n) == -1:
92
             return D
94
             return D
95
+        # This is supposed to be faster
96
+        if D > 30 and not haschecked:
97
+            if hasIntSQRT(n):
98
+                raise ValueError("n must not be a square")
93
     # Shouldn't fire
99
     # Shouldn't fire
94
     return 0
100
     return 0
95
 
101
 
99
     """
105
     """
100
     if n < 0:
106
     if n < 0:
101
         raise ValueError("n must be a non-negative integer")
107
         raise ValueError("n must be a non-negative integer")
102
-    if n == 0:
103
-        return 0, 2
104
-    Unm1, Vnm1 = Lucas(n - 1, p, q)
105
-    Un = (p * Unm1 + Vnm1) // 2
106
-    Vn = ((p * p - 4 * q) * Unm1 + p * Vnm1) // 2
107
-    return Un, Vn
108
+    Uc = 0
109
+    Vc = 2
110
+    for i in range(n):
111
+        Un = (p * Uc + Vc) // 2
112
+        Vn = ((p * p - 4 * q) * Uc + p * Vc) // 2
113
+        Uc = Un
114
+        Vc = Vn
115
+    return Uc, Vc
108
 
116
 
109
 def LucasUn(n, p, q):
117
 def LucasUn(n, p, q):
110
     """Return only the U numbers from a lucas sequence Un(P, Q). For example if
118
     """Return only the U numbers from a lucas sequence Un(P, Q). For example if
152
     return(uc, vc)
160
     return(uc, vc)
153
 
161
 
154
 
162
 
155
-for i in range(1, 20000, 2):
163
+for i in range(1, 20, 2):
156
     if hasIntSQRT(i):
164
     if hasIntSQRT(i):
157
         continue
165
         continue
158
     d = FirstD(i)
166
     d = FirstD(i)
161
     su, sv = Lucas(i + 1, 1, q)
169
     su, sv = Lucas(i + 1, 1, q)
162
     su = su % i
170
     su = su % i
163
     sv = sv % i
171
     sv = sv % i
164
-    print("{}: {}, {}; {}".format(i, fu == su, fv == sv, Lucas(i + 1, 1, q)))
165
-
166
-
172
+    if (fu != su) or (fv != sv):
173
+        print("{}: {}, {}; {} {}".format(i, fu, fv, su, sv))
174
+print("Done!")
175
+
176
+i = 19999999999999999999999999999999999999999999999999999999999
177
+t1 = process_time()
178
+hasIntSQRT(i)
179
+te = process_time() - t1
180
+print(te)
181
+t1 = process_time()
182
+d = FirstD(i)
183
+te = process_time() - t1
184
+print(te)
185
+q = (1 - d) // 4
186
+t1 = process_time()
187
+LucasFast(i + 1, 1, q, i)
188
+te = process_time() - t1
189
+print(te)
190
+t1 = process_time()
167
 
191
 
168
 
192