123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135 |
- # Finding lucas pseudoprimes
- from math import sqrt, gcd
-
- def isqrt(n):
- """ Find the integer square root of n via newton's method, code via
- stackoverflow:
-
- (https://stackoverflow.com/questions/15390807/integer-square-root-in-python)
- """
- x = n
- y = (x + 1) // 2
- while y < x:
- x = y
- y = (x + n // x) // 2
- return x
-
- def hasIntSQRT(n):
- """Detect whether the square root of n is an integer,
- i.e. whether the isqrt(n) is the true square root.
- """
- isq = isqrt(n)
- return isq * isq == n
-
- def Dsequence():
- """Generate sequence 5, -7, 9, -11, 13, -15...
- """
- val = 5
- while True:
- if val % 4 == 1:
- yield val
- else:
- yield -val
- val = val + 2
-
- def Legendre(a, p):
- """Function for calculating Legendre symbol.
-
- Note that this is only supposed to be defined if p is an odd prime, but we
- don't know in advance if that is true - will therefore only throw and error
- if p is even.
-
- Uses original power definition from
- <https://en.wikipedia.org/wiki/Legendre_symbol>
- """
- if (p % 2 == 0):
- raise ValueError("p must be odd, is {}".format(p))
- lv = pow(a, (p-1)//2, p)
- if lv == p - 1:
- lv = -1
- return lv
-
- def Jacobi(a, n):
- """Function for calculating Jacobi symbol.
-
- Note that this is only defined for positive odd integers n.
-
- Uses algorithm from
- <https://en.wikipedia.org/wiki/Jacobi_symbol#Calculating_the_Jacobi_symbol>
- """
- if n < 1:
- raise ValueError("n must be positive")
- if (n % 2 == 0):
- raise ValueError("n must be odd")
- if a % n != a:
- return Jacobi(a%n, n)
- if a != 0 and a % 2 == 0:
- nm8 = n % 8
- if (nm8 == 3 or nm8 == 5):
- return -1 * Jacobi(a//2, n)
- else:
- return Jacobi(a//2, n)
- if a == 1:
- return 1
- if gcd(a, n) != 1:
- return 0
- if a == 0 and n == 1:
- return 1
- if n % 4 == 3 and a % 4 == 3:
- return -1 * Jacobi(n, a)
- return Jacobi(n, a)
-
- def FirstD(n):
- """Return first D in the sequence 5, -7, 9, -11, 13, -15... for which the
- Jacobi symbol (D/n) is -1
- """
- if hasIntSQRT(n):
- raise ValueError("n must not be a square")
- 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:
- return D
- # Shouldn't fire
- return 0
-
- def Lucas(n, p, q):
- """Function for generating values of a lucas sequence with parameters p, q.
- Via formula at <https://en.wikipedia.org/wiki/Lucas_sequence>
- """
- if n < 0:
- raise ValueError("n must be a non-negative integer")
- if n == 0:
- return 0, 2
- Unm1, Vnm1 = Lucas(n - 1, p, q)
- Un = (p * Unm1 + Vnm1) // 2
- Vn = ((p * p - 4 * q) * Unm1 + p * Vnm1) // 2
- return Un, Vn
-
- def LucasUn(n, p, q):
- """Return only the U numbers from a lucas sequence Un(P, Q). For example if
- P = 1, and Q = -1 this will return the Fibonacci numbers; if P = 3 and Q =
- 2 then it returns the Mersenne numbers etc.
- """
- return Lucas(n, p, q)[0]
-
- def LucasVn(n, p, q):
- """Return only V numbers from a local sequence Vn(P, Q).
- """
- return Lucas(n, p, q)[1]
-
-
- def LucasFast(k, p, q, n):
- if k < 0:
- 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)
- dk = floor(log(k, 2))
- tp = 2**dk
- for i in range(dk):
-
-
-
-
|