From 6d63d22a03dff394f3c6e6d0f79b013619a65930 Mon Sep 17 00:00:00 2001 From: Gerald Teschl Date: Sun, 17 Nov 2024 17:22:06 +0100 Subject: [PATCH] Small improvements --- doc/kryptools.ipynb | 49 +++++++++++++++++++++-------------------- kryptools/__init__.py | 4 ++-- kryptools/dlp.py | 51 +++++++++++++++++++++++++++++++++---------- kryptools/dlp_bsgs.py | 29 +++++++++++++++++++++--- kryptools/factor.py | 2 +- kryptools/la.py | 12 +++++----- kryptools/primes.py | 15 ++++++++----- pyproject.toml | 2 +- tests/test_dlog.py | 5 +++-- 9 files changed, 113 insertions(+), 56 deletions(-) diff --git a/doc/kryptools.ipynb b/doc/kryptools.ipynb index 8130a06..c234f72 100644 --- a/doc/kryptools.ipynb +++ b/doc/kryptools.ipynb @@ -637,7 +637,7 @@ "id": "c3825781-b529-4fd3-8051-bc4d08850281", "metadata": {}, "source": [ - "The Prime-counting function $\\pi(x)$ is implemented by counting the number of primes returned by the sieve of Eratosthenes:" + "The Prime-counting function $\\pi(x)$ and the primorial $\\#x$ is implemented by counting the number of primes returned by the sieve of Eratosthenes:" ] }, { @@ -649,7 +649,7 @@ { "data": { "text/plain": [ - "8" + "(8, 210)" ] }, "execution_count": 21, @@ -658,9 +658,9 @@ } ], "source": [ - "from kryptools import prime_pi\n", + "from kryptools import prime_pi, primorial\n", "\n", - "prime_pi(20)" + "prime_pi(20), primorial(7)" ] }, { @@ -673,17 +673,17 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 126, "id": "2e105eeb-92a2-4b8a-86bc-3cb46c1d8b2e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "True" + "False" ] }, - "execution_count": 22, + "execution_count": 126, "metadata": {}, "output_type": "execute_result" } @@ -691,7 +691,7 @@ "source": [ "from kryptools import is_prime\n", "\n", - "is_prime(17)" + "is_prime(primorial(13)+1)" ] }, { @@ -887,7 +887,7 @@ { "data": { "text/plain": [ - "266" + "269" ] }, "execution_count": 29, @@ -1975,7 +1975,7 @@ "id": "72e95acf-ae6f-4d26-b22e-4baf15813016", "metadata": {}, "source": [ - "The coefficients can be from a field, e.g. `Fractions`. This can be done either by using `ring=Fraction` during creation or using the `map` method which applies a given function to all elements." + "The coefficients can be from a ring or field, e.g. `Fractions`. This can be done either by using `ring=Fraction` during creation or using the `map` method which applies a given function to all elements." ] }, { @@ -2145,21 +2145,12 @@ "id": "d9cfb165-6e94-413b-9a91-5c5aa3870aa3", "metadata": {}, "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[ 6, 0, 4 ]\n", - "[ 8, 1, 3 ]\n", - "[ 4, 3, 4 ]\n" - ] - }, { "data": { "text/plain": [ - "[ 1, 0, 0 ]\n", - "[ 0, 1, 0 ]\n", - "[ 0, 0, 1 ]" + "[ 6, 0, 4 ]\n", + "[ 8, 1, 3 ]\n", + "[ 4, 3, 4 ]" ] }, "execution_count": 68, @@ -2168,9 +2159,19 @@ } ], "source": [ + "from kryptools import eye\n", + "\n", "Mi = M.inv()\n", - "print(Mi)\n", - "M * Mi" + "assert M * Mi == eye(3, ring=gf)\n", + "Mi" + ] + }, + { + "cell_type": "markdown", + "id": "9efdd934-1fe0-4ca0-be46-50854d67ee81", + "metadata": {}, + "source": [ + "Note that `eye` will give you the identity matrix. There is also a function `zero` to create a zero matrix." ] }, { diff --git a/kryptools/__init__.py b/kryptools/__init__.py index 6838185..74f0191 100644 --- a/kryptools/__init__.py +++ b/kryptools/__init__.py @@ -5,11 +5,11 @@ __author__ = "Gerald Teschl" __copyright__ = "Copyright 2024, Gerald Teschl" __license__ = "MIT License" -__version__ = "0.9.13" +__version__ = "0.9.14" __email__ = "Gerald.Teschl@univie.ac.at" from .nt import egcd, cf, convergents, legendre_symbol, jacobi_symbol, sqrt_mod, euler_phi, carmichael_lambda, moebius_mu, is_carmichael_number, order, crt -from .primes import sieve_eratosthenes, prime_pi, is_prime, next_prime, previous_prime, next_safeprime, previous_safeprime, random_prime, random_strongprime, is_safeprime, random_safeprime, is_blumprime, random_blumprime, miller_rabin_test, lucas_test +from .primes import sieve_eratosthenes, prime_pi, primorial, is_prime, next_prime, previous_prime, next_safeprime, previous_safeprime, 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 diff --git a/kryptools/dlp.py b/kryptools/dlp.py index 14949a5..25ac778 100644 --- a/kryptools/dlp.py +++ b/kryptools/dlp.py @@ -4,63 +4,90 @@ """ -from math import gcd, log, sqrt +from math import gcd, log, log10, sqrt, floor from .nt import crt, order from .dlp_bsgs import dlog_bsgs from .dlp_qs import dlog_qs from .factor import factorint -def dlog_naive(a: int, b: int, n: int, m: int = None) -> int | None: +def dlog_naive(a: int, b: int, n: int, m: int = None, verbose: int = 0) -> int | None: """Compute the discrete log_a(b) in Z_n of an element a of order m by exhaustive search.""" a %= n b %= n if not m: m = n - 1 + if verbose: + print("Starting naive search: ", end="") + show = m // 100 + else: + show = 0 aa = 1 for i in range(m): + if show and not i % show and i: + print(".", end="") aa = (aa * a) % n if aa == b: + if verbose: + print("done") return i + 1 raise ValueError("DLP not solvable.") -def _dlog_switch(a: int, b: int, n: int, m: int) -> int: +def _dlog_switch(a: int, b: int, n: int, m: int, verbose: int = 0) -> int: """Compute the discrete log_a(b) in Z_n of an element a of order m choosing an appropriate method.""" if m < 1000: - return dlog_naive(a, b, n, m) + if verbose > 1: + print(f"Naive search: {a}^x = {b} (order={m})") + return dlog_naive(a, b, n, m, verbose=max(0, verbose - 1)) # compare the theoreticaly expected running times ob bsgs and ic; the constant 6 is determined experimentally if log(m) - 6 < 2 * sqrt(log(n) * log(log(n))): - return dlog_bsgs(a, b, n, m) - return dlog_qs(a, b, n, m) + if verbose > 1: + print(f"BS GS: {a}^x == {b} (order={m})") + return dlog_bsgs(a, b, n, m, verbose=max(0, verbose - 1)) + if verbose > 1: + print(f"Quadratic sieve: {a}^x == {b} (order={m})") + return dlog_qs(a, b, n, m, verbose=max(0, verbose - 1)) -def _dlog_ph(a: int, b: int, n: int, q: int, k: int) -> int: +def _dlog_ph(a: int, b: int, n: int, q: int, k: int, verbose: int = 0) -> int: """Compute the discrete log_a(b) in Z_n of an element a of order q^k using Pohlig-Hellman reduction.""" + if verbose > 1: + print(f"PH reduction: {a}^x = {b} (order={q}^{k})") if k == 1 or q**k < 10000: - return _dlog_switch(a, b, n, q**k) + return _dlog_switch(a, b, n, q**k, verbose = verbose) aj = pow(a, q ** (k - 1), n) a1 = aj bj = pow(b, q ** (k - 1), n) - xj = _dlog_switch(a1, bj, n, q) + xj = _dlog_switch(a1, bj, n, q, verbose = verbose) for j in range(2, k + 1): aj = pow(a, q ** (k - j), n) bj = pow(b, q ** (k - j), n) - yj = _dlog_switch(a1, bj * pow(aj, -xj, n) % n, n, q) # pylint: disable=E1130 + yj = _dlog_switch(a1, bj * pow(aj, -xj, n) % n, n, q, verbose = verbose) # pylint: disable=E1130 xj = xj + q ** (j - 1) * yj % q**j return xj -def dlog(a: int, b: int, n: int, m: int | None = None) -> int: +def dlog(a: int, b: int, n: int, m: int | None = None, verbose: int = 0) -> int: """Compute the discrete log_a(b) in Z_n of an element `a` of order `m` using Pohlig-Hellman reduction.""" a %= n b %= n assert gcd(a, n) == 1, "a and n must be coprime." assert gcd(b, n) == 1, "b and n must be coprime." if m: + if verbose: + print("Factoring the order.") mf = factorint(m) else: + if verbose: + print("Determining the order: ", end="") m, mf = order(a, n, True) + if verbose: + print(m) + if verbose: + print(f"Solving the DLP in Z_{n} ({floor(log10(n)) + 1} digits)") + print(f"Solving the DLP: {a}^x = {b}") + print(f"Order: {m}") if pow(b, m, n) != 1: raise ValueError("DLP not solvable.") # We first use Pohlig-Hellman to split m into powers of prime factors @@ -69,7 +96,7 @@ def dlog(a: int, b: int, n: int, m: int | None = None) -> int: for pj, kj in mf.items(): aj = pow(a, m // pj**kj, n) bj = pow(b, m // pj**kj, n) - l = _dlog_ph(aj, bj, n, pj, kj) + l = _dlog_ph(aj, bj, n, pj, kj, verbose = verbose) mm += [pj**kj] ll += [l] return crt(ll, mm) diff --git a/kryptools/dlp_bsgs.py b/kryptools/dlp_bsgs.py index 8e3cd39..97c887b 100644 --- a/kryptools/dlp_bsgs.py +++ b/kryptools/dlp_bsgs.py @@ -5,25 +5,48 @@ from math import isqrt -def dlog_bsgs(a: int, b: int, n: int, m: int = None) -> int: +def dlog_bsgs(a: int, b: int, n: int, m: int = None, verbose: int = 0) -> int: """Compute the discrete log_a(b) in Z_n of an element a of order m using Shanks' baby-step-giant-step algorithm.""" a %= n b %= n if not m: m = n - 1 mm = 1 + isqrt(m - 1) + # initialize baby_steps table + if verbose: + print("Computing baby steps", end="") + show = mm // 100 + if show: + print(": ", end="") + else: + show = 0 baby_steps = {} baby_step = 1 for k in range(mm): + if show and not k % show and k: + print(".", end="") baby_steps[baby_step] = k + if b == baby_step: + if verbose: + print(".") + return k baby_step = baby_step * a % n # now take the giant steps + if verbose: + print(".") + print("Computing giant steps", end="") + if show: + print(": ", end="") giant_stride = pow(a, -mm, n) giant_step = b - for l in range(mm): + for l in range(1, mm): + giant_step = (giant_step * giant_stride) % n + if show and not k % show and k: + print(".", end="") if giant_step in baby_steps: + if verbose: + print(".") return l * mm + baby_steps[giant_step] - giant_step = giant_step * giant_stride % n raise ValueError("DLP not solvable.") diff --git a/kryptools/factor.py b/kryptools/factor.py index 1c4965a..042acc2 100644 --- a/kryptools/factor.py +++ b/kryptools/factor.py @@ -42,7 +42,7 @@ def _factor_fermat(n: int, maxsteps: int = 10, verbose: int = 0) -> list: print("") -def factorint(n: int, verbose: int = 0, trial_bnd: int = 2500) -> dict: +def factorint(n: int, trial_bnd: int = 2500, verbose: int = 0) -> dict: "Factor an ineger `n` into prime factors." prime_factors = {} if not isinstance(n, int): diff --git a/kryptools/la.py b/kryptools/la.py index d1508c9..f97d308 100644 --- a/kryptools/la.py +++ b/kryptools/la.py @@ -315,19 +315,19 @@ def delta(i, j): return self.__class__([[delta(i, j) for j in range(n)] for i in range(m)]) -def zeros(m: int, n: int = None, zero=0) -> "Matrix": +def zeros(m: int, n: int = None, zero=0, ring=None) -> "Matrix": "Returns a zero matrix of the given dimension." if not n: n = m - return Matrix([[zero for j in range(n)] for i in range(m)]) + return Matrix([[zero for j in range(n)] for i in range(m)], ring=ring) -def eye(m: int, n: int = None) -> "Matrix": +def eye(m: int, n: int = None, zero=0, one=1, ring=None) -> "Matrix": "Returns an identity matrix of the given dimension." def delta(i, j): if i == j: - return 1 - return 0 + return one + return zero if not n: n = m - return Matrix([[delta(i, j) for j in range(n)] for i in range(m)]) + return Matrix([[delta(i, j) for j in range(n)] for i in range(m)], ring=ring) diff --git a/kryptools/primes.py b/kryptools/primes.py index 42f8257..8ecb06d 100644 --- a/kryptools/primes.py +++ b/kryptools/primes.py @@ -1,7 +1,8 @@ """ Tools for prime numbers: sieve_eratosthenes(B) a tuple of all primes up to including B - prime_pi(n) number of primes below or equal n + prime_pi(x) number of primes below or equal x + primorial(x) product of all primes below or equal x is_prime(n) test if n is probably prime next_prime(n) find the next prime larger or equal n previous_prime(n) find the previous prime smaller or equal n @@ -17,7 +18,7 @@ miller_rabin_test(n, b) Miller-Rabin primality test with base b lucas_test(n) strong Lucas test """ -from math import isqrt, gcd, floor, ceil +from math import isqrt, gcd, floor, ceil, prod #from random import getrandbits as randbits from secrets import randbits from .nt import jacobi_symbol @@ -25,7 +26,7 @@ # Erathostenes -def sieve_eratosthenes(B: int) -> tuple: +def sieve_eratosthenes(B: float) -> tuple: """Returns a tuple of all primes up to (including) `B`.""" if B < 2: return () @@ -35,13 +36,13 @@ def sieve_eratosthenes(B: int) -> tuple: B1 = (isqrt(B) - 1)//2 B = (B - 1)//2 # to begin with, all odd numbers are potentially prime - isprime = [True] * (B + 1) + isprime = bytearray([True]) * (B + 1) # bytearray makes it slightly faster # sieve out the primes p=2*q+1 starting at q = 1 for q in range(1, B1 + 1): if isprime[q]: # sieve out all multiples; numbers p*r with r

int: """Prime-counting function pi(x). Returns the number of primes below or equal `x`.""" return len(sieve_eratosthenes(x)) +def primorial(x: float) -> int: + """Primorial. Returns the product of all primes below or equal `x`.""" + return prod(sieve_eratosthenes(x)) + # Primality testing diff --git a/pyproject.toml b/pyproject.toml index 81c6510..f28dc2b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "kryptools" -version = "0.9.13" +version = "0.9.14" authors = [ { name="Gerald Teschl", email="gerald.teschl@univie.ac.at" }, ] diff --git a/tests/test_dlog.py b/tests/test_dlog.py index 0cff052..44dc13e 100644 --- a/tests/test_dlog.py +++ b/tests/test_dlog.py @@ -1,6 +1,6 @@ import pytest from random import randint, seed -from kryptools import dlog +from kryptools import dlog, primorial seed(0) @@ -12,7 +12,8 @@ def test_dlogt(): for data in ( [557639, 278819, 2], [24570203447, 12285101723, 2], - [28031135240181527, 14015567620090763, 2]): + [28031135240181527, 14015567620090763, 2], + [primorial(379)+1, primorial(379)//3, 383]): p, m, a = data # m is the order of a in Z_p x = randint(2, m - 1) b = pow(a, x, p)