12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697 |
- # 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)
-
- print(" ", end=" ")
- for k in range(1, 10):
- print("{}".format(k).rjust(2, " "), end=" ")
- print()
-
- for n in range(1, 21, 2):
- print("{}:".format(n).rjust(3, " "), end=" ")
- for k in range(1, 10):
- print("{}".format(Jacobi(k, n)).rjust(2, " "), end=" ")
- print()
-
-
-
-
-
|