Some extensions, better docs

Author: teschlg

doc/kryptools.ipynb

@@ -6,7 +6,7 @@
 from .ec import EC_Weierstrass
 from .factor import factorint
 from .la import Matrix
-from .lat import hermite_nf, gram_schmidt, lll
+from .lat import gram_det, hadamard_ratio, hermite_nf, gram_schmidt, babai_round_cvp, babai_plane_cvp, lagrange_lr, lll
 from .nt import cf, convergents, jacobi_symbol, sqrt_mod, euler_phi, order, carmichael_lambda
 from .poly import Poly
 from .primes import sieve_eratosthenes, isprime

kryptools/__init__.py

@@ -3,9 +3,10 @@
 """
 
 from math import gcd
-from random import randint
+from random import randint, seed
 from .primes import sieve_eratosthenes
 from .dlp_qs import determine_factorbound, determine_trialdivison_bounds, is_smooth
+seed(0)
 
 def dlog_ic(a: int, b: int, n: int, m: int, pollard: bool = True, verbose: int = 0) -> int:
     """Compute the discrete log_a(b) in Z_p of an element a of prime order m using Index Calculus."""

kryptools/dlp_ic.py

@@ -2,10 +2,11 @@
 Discrete log solvers: Quadratic sieve
 """
 
+from math import exp, log, sqrt, gcd, isqrt
+from random import randint, seed
 from .nt import sqrt_mod
 from .primes import sieve_eratosthenes
-from math import exp, log, sqrt, gcd
-
+seed(0)
 
 def determine_factorbound(n: int) -> (int, int):
     """Determines the optimal factor bound and the expected number of trys until a for a given n."""

kryptools/dlp_qs.py

@@ -0,0 +1,90 @@
+"""
+Integer factorization: Dixon's method
+"""
+
+from math import isqrt, gcd, sqrt, log, exp, ceil
+from .primes import sieve_eratosthenes
+from .nt import legendre_symbol, sqrt_mod
+from .factor_qs import bytexor, byteset, bytetest
+
+
+def factor_dixon(n: int) -> list:
+    """Find factors of n using the method of Dixon."""
+    # first determine the bound B for the factorbase: Choosing B=p^(1/u) Canfield-Erdös-Pomerance gives us
+    # the expected running time |B|^2 u^u = u^(u+2) p^(2/u)/log(n). There is no explicit expression for the optimum, hence
+    # we use Newton
+    u = 2 * sqrt(log(n) / log(log(n)))  # asymptotic value
+    for _ in range(3):
+        u = (2 * log(n) + u * u * (2 + log(u))) / (
+            2 + 3 * u + 2 * u * log(u)
+        )  # Newton iteration
+    B = int(exp(log(n) / u))
+    
+    # B = int(exp(0.5 * sqrt( log(n) * log(log(n)) )*( 1 + 1/log(log(n)) )))
+    factorbase = []
+    for p in sieve_eratosthenes(B):  # compute the factorbase
+        ls = legendre_symbol(n, p)
+        if ls == 0:  # we already found a factor;-)
+            if n == p:
+                return n
+            return p
+        if ls == 1:  # we only take primes such that n is quadratic residue
+            factorbase.append(p)
+    lf = len(factorbase) + 1  # length of the factorbase (including -1)
+    lfb = ceil(lf / 8)  # the number of bytes we need to store a relation
+    m = isqrt(n - 1) + 1
+    def process_relation(j: int, relation: bytes):
+        nonlocal relation_no, values, relations
+        relation_no += 1
+        rhs = bytearray(b"\x00") * lfb # construct the k'th row of the identity matrix
+        byteset(rhs, relation_no)
+        relation += rhs  # extend the relation with this row
+        values[relation_no] = j  # store the value which lead to the relation
+        # do the Gauss elimination
+        index = lf # this will be the index of the first nonzero entry
+        # print(f'{j:3}', ' '.join(f'{b:08b}' for b in reversed(relation)))
+        for i in range(lf):
+            if bytetest(relation, i) and relations[i] != None:  # make this entry zero if we can (Gauss elimination)
+                bytexor(relation, relations[i])
+            if bytetest(relation, i) and index == lf:  # is this the index of the first nonzero entry?
+                index = i
+        # print(f'{j:3}', ' '.join(f'{b:08b}' for b in reversed(relation)))
+        if index == lf:  # the new relation is linearly dependent: we have found a linear combination of the 0 vector
+            # now we need to determine the factors
+            u = 1  # product over all values m - j such that f(m-j) is B-smooth
+            v = 1  # sqrt of the product over all f(m-j) which are B-smooth
+            w = 1  # save nonsquare terms for the next round
+            for i in range(lf):
+                if bytetest(relation, lfb * 8 + i):  # select the relations which sum to the 0 vector
+                    ui = m + values[i]
+                    u = (u * ui) % n
+                    vi = ui ** 2 - n
+                    d = gcd(w, vi)
+                    v = (v * d) % n  # sqrt of the part which is already square
+                    w = (w // d) * (vi // d)  # this part is not square yet
+            v = v * isqrt(w) % n
+            res = gcd(u - v, n)
+            if 1 < res < n:
+                return res
+            relation_no -= 1  # this one did not work, try again
+        else:
+            relations[index] = relation
+        return None
+
+    relation_no = -1
+    relations = [None] * lf  # here we will store the relations in case we found a B-smooth number
+    values = [None] * lf  # here we will store the values leading to the B-smooth numbers
+    for j in range(1,m): # loop over j=0,1,-1,2,-2,...
+        if j & 1 == 0:
+            j >>= 1
+        else:
+            j >>= 1
+            j *= -1
+        relation = is_smooth((j + m) ** 2 - n, factorbase, lfb)  # test if f(j) is B-smooth
+        if relation is None:
+            continue
+        res = process_relation(j, relation)
+        if res:
+            return(res)
+    return n
+

kryptools/factor_dix.py

@@ -0,0 +1,30 @@
+"""
+Integer factorization: Fermat's method
+"""
+
+from math import isqrt
+
+
+def factor_fermat(n: int) -> list:
+    """Find factors of n using the method of Fermat."""
+    factors = []
+    # Fermat only works if n has two factors which are either both even or both odd
+    while n % 2 == 0:
+        factors.append(2)
+        n //= 2
+    a = isqrt(n - 1) + 1
+    step =2
+    if n % 3 == 2:  # if n % 3 = 2, then a must be a multiple of 3
+        a += 2 - ((a - 1) % 3)
+        step = 3
+    elif (n % 4 == 1) ^ (a & 1):  # if n % 4 = 1,3 then a must be odd, even, respectively
+        a += 1
+    while a <= (n + 9) // 6:
+        b = isqrt(a * a - n)
+        if b * b == a * a - n:
+            factors.append(a - b)
+            factors.append(a + b)
+            return factors
+        a += step
+    factors.append(n)
+    return factors

kryptools/factor_fmt.py

@@ -6,12 +6,6 @@
 from fractions import Fraction
 from .la import Matrix, zeros, eye
 
-def babai_round_cvp(x: Matrix, U: Matrix) -> Matrix:
-    "Babai's rounding algorithm for solving the CVP."
-    s = U.inv() * x
-    k = s.map(round)
-    return U * k
-
 def hermite_nf(M: Matrix) -> Matrix:
     "Compute the Hermite normal form of a matrix M."
     n, m = M.cols, M.rows
@@ -68,13 +62,21 @@ def gram_schmidt(U: Matrix) -> (Matrix, Matrix):
     return Us, M
 
 def gram_det(U: Matrix) -> float:
+    "Compute the Gram determinant of a matrix."
     Us = gram_schmidt(U)[0]
     return prod([Us[:, i].norm() for i in range(U.rows)])
 
 def hadamard_ratio(M: Matrix) -> float:
+    "Compute the Hadamard ratio of a matrix."
     m = M.rows
     return (gram_det(M) / prod([M[:, i].norm() for i in range(m)])) ** (1 / m)
 
+def babai_round_cvp(x: Matrix, U: Matrix) -> Matrix:
+    "Babai's rounding algorithm for solving the CVP."
+    s = U.inv() * x
+    k = s.applyfunc(round)
+    return U * k
+
 def babai_plane_cvp(x: Matrix, U: Matrix) -> Matrix:
     "Babai's closest plane algorithm for solving the CVP."
     Us = gram_schmidt(U)[0]
@@ -83,6 +85,18 @@ def babai_plane_cvp(x: Matrix, U: Matrix) -> Matrix:
         y = y - round(y.dot(Us[:, k]) / norm2(Us[:, k])) * U[:, k]
     return (x - y).applyfunc(round)
 
+def lagrange_lr(V: Matrix) -> Matrix:
+    "Lagrange lattice reduction."
+    assert (V.rows, V.cols) == (2, 2)
+    v1, v2 = V[:, 0], V[:, 1]
+    if norm2(v1) > norm2(v2):
+        v1, v2 = v2, v1
+    v3 = v2 - round(v1.dot(v2) / norm2(v1)) * v1
+    while norm2(v3) < norm2(v1):
+        v2, v1 = v1, v3
+        v3 = v2 - round(v1.dot(v2) / norm2(v1)) * v1
+    return Matrix([list(v1), list(v3)]).transpose()
+
 def lll(V: Matrix, delta: float = 0.75, sort: bool = True) -> Matrix:
     "lll algorithm for lattice reduction"
 

kryptools/lat.py

@@ -1,6 +1,6 @@
 [project]
 name = "kryptools"
-version = "0.1"
+version = "0.2"
 authors = [
   { name="Gerald Teschl", email="gerald.teschl@univie.ac.at" },
 ]