working fast lucas algo

Note: non-fast lucas stack overflows, will re-implement

lucas.py

@@ -1,5 +1,5 @@
 # Finding lucas pseudoprimes
-from math import sqrt, gcd
+from math import sqrt, gcd, floor, log
 
 def isqrt(n):
     """ Find the integer square root of n via newton's method, code via
@@ -88,7 +88,7 @@ def FirstD(n):
     if n < 1 or n % 2 == 0:
         raise ValueError("n must be a positive odd number")
     for D in Dsequence():
-        if Jacobi(D/n) == -1:
+        if Jacobi(D, n) == -1:
             return D
     # Shouldn't fire
     return 0
@@ -124,11 +124,44 @@ def LucasFast(k, p, q, n):
         raise ValueError("k must be a non-negative integer")
     if n < 1:
         raise ValueError("n must be a postive integer")
-    uc, vc = Lucas(1, p, q)
+    if n % 2 == 0:
+        raise ValueError("n must be odd")
+    uc, vc = Lucas(0, p, q)
     dk = floor(log(k, 2))
     tp = 2**dk
-    for i in range(dk):
-
+    krunn = k
+    ks = 0
+    for i in range(dk + 1):
+        un = (uc * vc)
+        vn = (pow(vc, 2, n) - 2 * pow(q, ks, n))
+        uc = un % n
+        vc = vn % n
+        ks = ks * 2
+        if krunn >= tp:
+            krunn = krunn - tp
+            ks = ks + 1
+            un = p * uc + vc
+            if un % 2 != 0:
+                un = un + n
+            vn = (p * p - 4 * q) * uc + p * vc
+            if vn % 2 != 0:
+                vn = vn + n
+            uc = (un // 2) % n
+            vc = (vn // 2) % n
+        tp = tp // 2
+    return(uc, vc)
+
+
+for i in range(1, 20000, 2):
+    if hasIntSQRT(i):
+        continue
+    d = FirstD(i)
+    q = (1 - d) // 4
+    fu, fv = LucasFast(i + 1, 1, q, i)
+    su, sv = Lucas(i + 1, 1, q)
+    su = su % i
+    sv = sv % i
+    print("{}: {}, {}; {}".format(i, fu == su, fv == sv, Lucas(i + 1, 1, q)))