Small improvements and further tests

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

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

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

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,

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"

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)
+

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<p were already sieved out previously
+			pp = p * p
+			p2 = p + p
+			is_prime[pp :: p2] = [False] * ((max - pp) // p2 + 1)
+	return is_prime
+
+prime_test = sieve(OEIS_A000010[-1])
+assert [ p for p in range(OEIS_A000010[-1]+1) if prime_test[p] ] == OEIS_A000010
 
 def test_sieve_eratosthenes():
 	for i, p in enumerate(sieve_eratosthenes(OEIS_A000010[-1])):
@@ -25,3 +52,51 @@ def test_sieve_eratosthenes():
 def test_prime_pi():
 	for n in range(1, len(OEIS_A000720)+1):
 		assert prime_pi(n) == OEIS_A000720[n - 1]
+
+max_sieve = 25326001
+prime_test = sieve(max_sieve)
+
+def test_is_prime():
+	for n in range(max_sieve + 1):
+		assert is_prime(n) == prime_test[n]
+
+OEIS_A005385 = [ 5, 7, 11, 23, 47, 59, 83, 107, 167, 179, 227, 263, 347, 359, 383, 467,
+	479, 503, 563, 587, 719, 839, 863, 887, 983, 1019, 1187, 1283, 1307, 1319, 1367, 1439,
+	1487, 1523, 1619, 1823, 1907, 2027, 2039, 2063, 2099, 2207, 2447, 2459, 2579, 2819,
+	2879, 2903, 2963 ]
+
+def test_is_safeprime():
+	for n in range(OEIS_A000720[-1]+2):
+		assert is_safeprime(n) == (n in OEIS_A005385)
+
+OEIS_A014233 = [ 2047, 1373653, 25326001, 3215031751, 2152302898747, 3474749660383,
+	341550071728321, 341550071728321, 3825123056546413051, 3825123056546413051,
+	3825123056546413051, 318665857834031151167461, 3317044064679887385961981 ]
+
+if OEIS_A014233[2] > 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