Small improvements

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,25 +673,25 @@
   },
   {
    "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"
     }
    ],
    "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."
    ]
   },
   {

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__ = "[email protected]"
 
 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

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)

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.")

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

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)

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,15 +18,15 @@
     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
 
 
 # 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<p were already sieved out previously
             qq = 2 * q * (q + 1)
             p = 2 * q + 1
-            isprime[qq:: p] = [False] * ((B - qq) // p + 1)
+            isprime[qq:: p] = bytearray([False]) * ((B - qq) // p + 1)
 
     return tuple([2] + [2 * q + 1 for q in range(1, B + 1) if isprime[q]])
 
@@ -50,6 +51,10 @@ def prime_pi(x: float) -> 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
 

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "kryptools"
-version = "0.9.13"
+version = "0.9.14"
 authors = [
   { name="Gerald Teschl", email="[email protected]" },
 ]