diff --git a/kryptools/__init__.py b/kryptools/__init__.py index 15eeb94..a99d7e9 100644 --- a/kryptools/__init__.py +++ b/kryptools/__init__.py @@ -10,6 +10,7 @@ from .nt import egcd, cf, convergents, legendre_symbol, jacobi_symbol, sqrt_mod, euler_phi, carmichael_lambda, moebius_mu, order, crt from .primes import sieve_eratosthenes, prime_pi, is_prime, next_prime, previous_prime, random_prime, random_strongprime, is_safeprime, random_safeprime, is_blumprime, random_blumprime, miller_rabin_test, lucas_test +from .intfuncs import iroot, ilog, perfect_square, perfect_power, prime_power from .factor import factorint, divisors from .dlp import dlog from .ec import EC_Weierstrass diff --git a/kryptools/factor.py b/kryptools/factor.py index 0146f44..bd1dfc7 100644 --- a/kryptools/factor.py +++ b/kryptools/factor.py @@ -172,16 +172,21 @@ def add_factors(m: int, mm: tuple) -> None: def divisors(n:int, proper:bool = False) -> int: """Returns an unsorted list of all divisors of an integer `n`.""" + if not isinstance(n, int): + raise ValueError("Number must be an integer!") + n = abs(n) + if proper: + divisors = [] + else: + divisors = [1] + if n <= 1: + return divisors facctordict= factorint(n) primes = [p for p in facctordict] nprimes = len(primes) multiplicites = [facctordict[p] for p in primes] exponents = [0] * nprimes - if proper: - divisors = [] - else: - divisors = [1] i = 0 while True: exponents[i] += 1 diff --git a/kryptools/intfuncs.py b/kryptools/intfuncs.py new file mode 100644 index 0000000..875aa4e --- /dev/null +++ b/kryptools/intfuncs.py @@ -0,0 +1,97 @@ +""" +Integer Functions: + iroot(k, n) integer k'th root of n + ilog(b, n) integer log of n with respect ot base b + perfect_square(n) test if a number is a perfect square + perfect_power(n) test if a number is a perfect power + prime_power(n) test if a number is a prime power +""" +from math import isqrt +from .primes import sieve_eratosthenes, is_prime + + +def iroot(k: int, n: int) -> int: + """For a given positive integers `n` and `k`, finds the largest integer `r` such that `r**k <= n`.""" + if not isinstance(k, int) or k < 1: + raise ValueError(f"k={k} must be a positive integer.") + if not isinstance(n, int) or n < 0: + raise ValueError(f"n={n} must be a nonegative integer.") + if n <= 1 or k == 1: + return n + if k == 2: + return isqrt(n) + rr = n + kk = k + while kk > 1: + kk //= 2 + rr = isqrt(rr) + r = rr + 1 + k1 = k-1 + while rr < r: + r = rr + rr = (k1 * rr + n // rr ** k1) // k # Newton iteration + return r + +def ilog(b: int, n:int) -> int: + """For a given positive integer `n`, finds the largest integer `l` such that `b**l <= n`.""" + # https://stackoverflow.com/questions/39190815/how-to-make-perfect-power-algorithm-more-efficient/39191163#39191163 + if not isinstance(b, int) or b < 2: + raise ValueError(f"base b={b} must a positive integer larger then one.") + if not isinstance(n, int) or n <= 0: + raise ValueError(f"n={n} must be a positive integer.") + if n == 1: + return 0 + if b == 2: + return n.bit_length() - 1 + lo, blo, hi, bhi = 0, 1, 1, b + while bhi < n: + lo, blo, hi, bhi = hi, bhi, hi + hi, bhi * bhi + while 1 < (hi - lo): + mid = (lo + hi) // 2 + bmid = blo * pow(b, mid - lo) + if n < bmid: + hi, bhi = mid, bmid + elif bmid < n: + lo, blo = mid, bmid + else: return mid + if bhi == n: return hi + return lo + + +def perfect_square(n: int) -> int|None: + """Returns the root if `n` is a perfect square and None else.""" + s = isqrt(n) + if s*s == n: + return s + +def perfect_power(n: int) -> tuple|None: + """Returns integers `(m, p)` with `m**p == n` if `n` is a perfect power and None else.""" + if n == 0: + return 0, 1 + if n == 1: + return 1, 1 + for p in sieve_eratosthenes(n.bit_length() - 1): + m = iroot(p, n) + if pow(m, p) == n: + return m, p + +def prime_power(n: int) -> tuple|None: + """Returns integers `(p, k)` with `p**k == n` if `n` is a prime power and None else.""" + if is_prime(n): + return (n, 1) + r = n + k = 1 + while r.bit_length() > 20: + k += 1 + r = iroot(k, n) + if r**k == n and is_prime(r): + return (r, k) + for p in sieve_eratosthenes(r - 1): + k = 0 + while n % p == 0: + k += 1 + n //= p + if k > 1: + if n == 1: + return (p, k) + return None diff --git a/kryptools/primes.py b/kryptools/primes.py index 48736eb..d13ba39 100644 --- a/kryptools/primes.py +++ b/kryptools/primes.py @@ -13,9 +13,6 @@ random_blumprime(l) find a pseudorandom Blum prime with bit length at least l miller_rabin_test(n, b) Miller-Rabin primality test with base b lucas_test(n) strong Lucas test - perfect_square(n) test if a number is a perfect square - perfect_power(n) test if a number is a perfect power - iroot(k, n) integer k'th root of n """ from math import isqrt, gcd, floor, ceil from random import randint @@ -110,40 +107,6 @@ def _lucas_sequence(n, D, k): Qk %= n return U % n, V % n, Qk - -def perfect_square(n: int) -> int|None: - """Returns the root if `n` is a perfect square and None else.""" - s = isqrt(n) - if s*s == n: - return s - -def iroot(k: int, n: int) -> int: - """Compute the integer k'th root of n.""" - if not isinstance(k, int) or n < 1: - raise ValueError(f"{k} is not a positive integer") - if k == 1: - return n - if k == 2: - return isqrt(n) - r = isqrt(n) - while r**k >= n: - rr = r - r = isqrt(r) - r =rr + 1 - k1 = k-1 - while rr < r: - r = rr - rr = (k1 * rr + n // rr ** k1) // k # Newton iteration - return r - -def perfect_power(n: int) -> tuple|None: - """Returns integers `(m, p)` with `m**p == n` if `n` is a perfect power and None else.""" - for p in sieve_eratosthenes(n.bit_length() - 1): - m = iroot(p, n) - if pow(m, p) == n: - return (m, p) - - # https://en.wikipedia.org/wiki/Lucas_pseudoprime#Implementing_a_Lucas_probable_prime_test def lucas_test(n: int) -> bool: @@ -152,7 +115,8 @@ def lucas_test(n: int) -> bool: return False if n % 2 == 0: return n == 2 - if perfect_square(n) is not None: # the search for D will not succeed in this case + s = isqrt(n) + if s * s == n: # the search for D will not succeed in this case return False # write n = k 2^s - 1 @@ -189,6 +153,8 @@ def lucas_test(n: int) -> bool: def is_prime(n: int, trialdivision: bool = True) -> bool: """Test if an integer `n` if probable prime.""" + if not isinstance(n, int) or n < 2: + return False # Test small primes if trialdivision: for p in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, diff --git a/pyproject.toml b/pyproject.toml index fa0ff26..af1969b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,6 +6,7 @@ authors = [ ] description = "Implemenation of same basic algorithms used in cryptography." readme = "README.md" +license = { text = "MIT" } requires-python = ">=3.10" classifiers = [ "Programming Language :: Python :: 3", @@ -15,5 +16,6 @@ classifiers = [ [project.urls] Homepage = "https://github.com/teschlg/kryptools" +Source = "https://github.com/teschlg/kryptools" Issues = "https://github.com/teschlg/kryptools/issues" Docs = "https://github.com/teschlg/kryptools/tree/main/doc" \ No newline at end of file diff --git a/tests/test_intfuncs.py b/tests/test_intfuncs.py new file mode 100644 index 0000000..97e68ed --- /dev/null +++ b/tests/test_intfuncs.py @@ -0,0 +1,39 @@ +import pytest +from random import randint, seed +from kryptools import sieve_eratosthenes +from kryptools import iroot, ilog, perfect_square, perfect_power, prime_power +seed(0) + + +def test_iroot(): + assert iroot(3, 0) == 0 + assert iroot(3, 1) == 1 + assert iroot(1, 3) == 3 + for _ in range(100): + n = randint(0, 10**10) + for k in range(2, 5): + r = iroot(k, n) + assert r ** k <= n and (r+1) ** k > n + +def test_ilog(): + for n in range(1, 10000): + for b in range(2,5): + l = ilog(b, n) + assert b**l <= n and b**(l+1) > n + for _ in range(1000): + n = randint(1, 10**10) + for b in range(2,5): + l = ilog(b, n) + assert b**l <= n and b**(l+1) > n + +def test_perfect_power(): + for n in (3* 7, 4 * 3, 4 * 3 * 5 * 7): + assert perfect_power(n) == None + for k in sieve_eratosthenes(11): + assert perfect_power(n**k) == (n, k) + +def test_prime_power(): + for p in sieve_eratosthenes(30): + for k in range(1,6): + assert prime_power(p**k) == (p, k) + \ No newline at end of file diff --git a/tests/test_primes.py b/tests/test_primes.py index e40162d..d40aba7 100644 --- a/tests/test_primes.py +++ b/tests/test_primes.py @@ -1,12 +1,39 @@ import pytest -from kryptools import sieve_eratosthenes, prime_pi, next_prime, previous_prime - +from math import isqrt +from random import randint, seed +from kryptools import sieve_eratosthenes, prime_pi, next_prime, previous_prime, is_prime, is_safeprime, miller_rabin_test, lucas_test +seed(0) #https://oeis.org/A000040 OEIS_A000010 = [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, - 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, - 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, - 269, 271 ] + 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, + 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, + 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, + 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, + 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, + 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, + 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, + 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, + 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, + 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, + 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223 ] + +def sieve(max: int) -> list: + is_prime = [True] * (max + 1) # to begin with, all numbers are potentially prime + is_prime[0 :: 2] = [False] * (max // 2 + 1) + is_prime[1] = False + is_prime[2] = True + + # sieve out the primes starting at 3 in steps of 2 (ignoring even numbers) + for p in range(3, isqrt(max - 1) + 2, 2): + if is_prime[p]: # sieve out all multiples; note that numbers p*q with q

max_sieve: + prime_test = sieve(OEIS_A014233[2]) + +def test_miller_rabin(): + for i, n in enumerate(OEIS_A014233[0:3]): + bases = OEIS_A000010[:i+1] + assert miller_rabin_test(n, bases) == True + if i < 1: + for m in range(3, n, 2): + assert miller_rabin_test(m, bases) == prime_test[m] + else: + for _ in range(500 * i): + r = randint(2, n) + assert miller_rabin_test(m, bases) == prime_test[m] + +OEIS_A217255 = [ 5459, 5777, 10877, 16109, 18971, 22499, 24569, 25199, 40309, 58519, + 75077, 97439, 100127, 113573, 115639, 130139, 155819, 158399, 161027, 162133, 176399, + 176471, 189419, 192509, 197801, 224369, 230691, 231703, 243629, 253259, 268349, + 288919, 313499, 324899 ] + +if OEIS_A217255[30] > max_sieve: + prime_test = sieve(OEIS_A217255[30]) + +def test_lucas(): + for m in range(2, OEIS_A217255[30]+1): + if not (lucas_test(m) == prime_test[m]): + assert m in OEIS_A217255