implemented significant improvements to the primality subsection

imp/math/factor/__init__.py

@@ -0,0 +1,2 @@
+from .pftrialdivision import trial_division
+from .factordb import DBResult, FType, FCertainty

imp/math/factor/factordb.py

@@ -0,0 +1,161 @@
+'''
+Simple interface for https://factordb.com inspired by
+https://github.com/ihebski/factordb
+
+TODO:
+1. Implement primality certificate generation, read this:
+https://reference.wolfram.com/language/PrimalityProving/ref/ProvablePrimeQ.html
+'''
+
+import requests
+from enum import Enum, StrEnum
+
+_FDB_URI = 'https://factordb.com'
+# Generated by https://www.asciiart.eu/text-to-ascii-art
+# using "ANSI Shadow" font and "Box drawings double" border with 1 H. Padding
+_BANNER = '''
+╔═══════════════════════════════════════════════════╗
+║ ███████╗ ██████╗████████╗██████╗ ██████╗ ██████╗  ║
+║ ██╔════╝██╔════╝╚══██╔══╝██╔══██╗██╔══██╗██╔══██╗ ║
+║ █████╗  ██║        ██║   ██████╔╝██║  ██║██████╔╝ ║
+║ ██╔══╝  ██║        ██║   ██╔══██╗██║  ██║██╔══██╗ ║
+║ ██║     ╚██████╗   ██║   ██║  ██║██████╔╝██████╔╝ ║
+║ ╚═╝      ╚═════╝   ╚═╝   ╚═╝  ╚═╝╚═════╝ ╚═════╝  ║
+╚═══════════════════════════════════════════════════╝
+╚═══ cli by imbored ═══╝
+'''.strip()
+
+# Enumeration of number types based on their factorisation
+class FType(Enum):
+    Unit      = 1
+    Composite = 2
+    Prime     = 3
+    Unknown   = 4
+
+class FCertainty(Enum):
+    Certain = 1
+    Partial = 2
+    Unknown = 3
+
+# FactorDB result status codes
+class DBStatus(StrEnum):
+    C    = 'C'
+    CF   = 'CF'
+    FF   = 'FF'
+    P    = 'P'
+    PRP  = 'PRP'
+    U    = 'U'
+    Unit = 'Unit' # just for 1
+    N    = 'N'
+    Add  = '*'
+
+    def is_unknown(self) -> bool:
+        return (self in [DBStatus.U, DBStatus.N, DBStatus.Add])
+
+    def classify(self) -> tuple[FType, FCertainty]:
+        return _STATUS_MAP[self]
+    
+    def msg_verbose(self) -> str:
+        return _STATUS_MSG_VERBOSE[self]
+
+# Map of DB Status codes to their factorisation type and certainty 
+_STATUS_MAP = {
+    DBStatus.Unit: (FType.Unit,      FCertainty.Certain),
+    DBStatus.C:    (FType.Composite, FCertainty.Unknown),
+    DBStatus.CF:   (FType.Composite, FCertainty.Partial),
+    DBStatus.FF:   (FType.Composite, FCertainty.Certain),
+    DBStatus.P:    (FType.Prime,     FCertainty.Certain),
+    DBStatus.PRP:  (FType.Prime,     FCertainty.Partial),
+    DBStatus.U:    (FType.Unknown,   FCertainty.Unknown),
+    DBStatus.N:    (FType.Unknown,   FCertainty.Unknown),
+    DBStatus.Add:  (FType.Unknown,   FCertainty.Unknown),
+}  
+
+# Reference: https://factordb.com/status.html
+# NOTE: my factor messages differ slightly from the reference
+_STATUS_MSG_VERBOSE = {
+    DBStatus.Unit: 'Unit, trivial factorisation',
+    DBStatus.C:    'Composite, no factors known',
+    DBStatus.CF:   'Composite, *partially* factors',
+    DBStatus.FF:   'Composite, fully factored',
+    DBStatus.P:    'Prime',
+    DBStatus.PRP:  'Probable prime',
+    DBStatus.U:    'Unknown (*but in database)',
+    DBStatus.N:    'Not in database (-not added due to your settings)',
+    DBStatus.Add:  'Not in database (+added during request)',
+}
+
+# Struct for storing database results with named properties
+class DBResult:
+    def __init__(self, 
+                 status: DBStatus, 
+                 factors: tuple[tuple[int, int]]) -> None:
+        self.status = status
+        self.factors = factors
+        self.ftype, self.certainty = self.status.classify()
+
+def _make_cookie(fdbuser: str | None) -> dict[str, str]:
+    return {} if fdbuser is None else {'fdbuser': fdbuser}
+
+def _get_key(by_id: bool):
+    return 'id' if by_id else 'query'
+
+def _api_query(n: int, 
+               fdbuser: str | None, 
+               by_id: bool = False) -> requests.models.Response:
+    key = _get_key(by_id)
+    uri = f'{_FDB_URI}/api?{key}={n}'
+    return requests.get(uri, cookies=_make_cookie(fdbuser))
+
+def _report_factor(n: int,
+                   factor: int,
+                   fdbuser: str | None,
+                   by_id: bool = False) -> requests.models.Response:
+    key = _get_key(by_id)
+    uri = f'{_FDB_URI}/reportfactor.php?{key}={n}&factor={factor}'
+    return requests.get(uri, cookies=_make_cookie(fdbuser))
+
+'''
+Attempts a query to FactorDB, returns a DBResult object 
+on success, or None if the request failed due to the
+get request raising a RequestException.
+'''
+def query(n: int, 
+          token: str | None = None, 
+          by_id: bool = False,
+          cli: bool = False) -> DBResult | None:
+    if cli:
+        print(_BANNER)
+    try:
+        resp = _api_query(n, token, by_id=by_id)
+    except requests.exceptions.RequestException:
+        return None
+    content = resp.json()
+    result = DBResult(
+        DBStatus(content['status']),
+        tuple((int(F[0]), F[1]) for F in content['factors'])
+    )
+    if cli:
+        print(f'Status: {result.status.msg_verbose()}')
+        print(result.factors)
+
+    # ensure the unit has the trivial factorisation (for consistency)
+    if result.status == DBStatus.Unit:
+        result.factors = ((1, 1),)
+    return result
+
+'''
+Reports a known factor to FactorDB, also tests it is
+actually a factor to avoid wasting FactorDBs resources. 
+'''
+def report(n: int, 
+           factor: int, 
+           by_id: int,
+           token: str | None = None) -> None:
+    try:
+        resp = _report_factor(n, factor, token)
+    except requests.exceptions.RequestException:
+        return None
+    content = resp.json()
+    print(content)
+

imp/math/factor/pftrialdivision.py

@@ -0,0 +1,46 @@
+'''
+The trial division algorithm is essentially the idea that
+all factors of an integer n are less than or equal to isqrt(n),
+where isqrt is floor(sqrt(n)).
+
+Moreover, if p divides n, then all other factors of n must
+be factors of n//p. Hence they must be <= isqrt(n//p).
+'''
+from math import isqrt # integer square root
+
+# Returns the "multiplicity" of a prime factor
+def pf_multiplicity(n: int, p: int) -> int:
+    mult = 0
+    while n % p == 0: 
+        n //= p
+        mult += 1
+    return n, mult
+
+'''
+Trial division prime factorisation algorithm.
+Returns a list of tuples (p, m) where p is
+a prime factor and m is its multiplicity.
+'''
+def trial_division(n: int) -> list[tuple[int, int]]:
+    factors = []
+
+    # determine multiplicity of the only even prime (2)
+    n, mult_2 = pf_multiplicity(n, 2)
+    if mult_2: factors.append((2, mult_2))
+    # determine odd factors and their multiplicities
+    p = 3
+    mult = 0
+    limit = isqrt(n)
+    while p <= limit:
+        n, mult = pf_multiplicity(n, p)
+        if mult:
+            factors.append((p, mult))
+            limit = isqrt(n) # recalculate limit
+            mult = 0 # reset
+        else:
+            p += 2
+    # if n is still greater than 1, then n is a prime factor
+    if n > 1:
+        factors.append((n, 1))
+    return factors
+            

imp/math/primes.py

@@ -1,67 +0,0 @@
-from math import gcd
-
-'''
-Euler's Totient (Phi) Function
-'''
-def totient(n: int) -> int:
-    phi = int(n > 1 and n)
-    for p in range(2, int(n ** .5) + 1):
-        if not n % p:
-            phi -= phi // p
-            while not n % p:
-                n //= p
-    #if n is > 1 it means it is prime
-    if n > 1: phi -= phi // n 
-    return phi
-
-'''
-Tests the primality of an integer using its totient.
-NOTE: If totient(n) has already been calculated 
-      then pass it as the optional phi parameter. 
-'''
-def is_prime(n: int, phi: int = None) -> bool:
-    return n - 1 == (phi if phi is not None else totient(n))
-
-'''
-Prime number generator function.
-Returns the tuple (p, phi(p)) where p is prime
-    and phi is Euler's totient function.
-'''
-def prime_gen(yield_phi: bool = False) -> int | tuple[int, int]:
-    n = 1
-    while True:
-        n += 1
-        phi = totient(n)
-        if is_prime(n, phi=phi):
-            if yield_phi:
-                yield (n, phi)
-            else:
-                yield n
-
-'''
-Returns the prime factorisation of a number.
-Returns a list of tuples (p, m) where p is
-a prime factor and m is its multiplicity.
-NOTE: uses a trial division algorithm
-'''
-def prime_factors(n: int) -> list[tuple[int, int]]:
-    phi = totient(n)
-    if is_prime(n, phi=phi):
-        return [(n, 1)]
-    factors = []
-    for p in prime_gen(yield_phi=False):
-        if p >= n:
-            break
-        # check if divisor
-        multiplicity = 0
-        while n % p == 0:
-            n //= p
-            multiplicity += 1
-        if multiplicity:
-            factors.append((p, multiplicity))
-            if is_prime(n):
-                break
-    if n != 1:
-        factors.append((n, 1))
-    return factors
-            

imp/math/primes/__init__.py

@@ -0,0 +1,48 @@
+from math import inf, isqrt # integer square root
+from itertools import takewhile, compress
+
+SMALL_PRIMES = (2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59)
+
+'''
+Euler's Totient (Phi) Function 
+Implemented in O(nloglog(n)) using the Sieve of Eratosthenes.
+'''
+def eulertotient(n: int) -> int:
+    phi = int(n > 1 and n)
+    for p in range(2, isqrt(n) + 1):
+        if not n % p:
+            phi -= phi // p
+            while not n % p:
+                n //= p
+    #if n is > 1 it means it is prime
+    if n > 1: phi -= phi // n 
+    return phi
+
+'''
+Tests the primality of an integer using its totient.
+NOTE: If totient(n) has already been calculated 
+      then pass it as the optional phi parameter. 
+'''
+def is_prime(n: int, phi: int = None) -> bool:
+    return n - 1 == (phi if phi is not None else eulertotient(n))
+
+# Taken from Lucas A. Brown's primefac.py (some variables renamed)
+def primegen(limit=inf) -> int:
+    ltlim = lambda x: x < limit
+    yield from takewhile(ltlim, SMALL_PRIMES)
+    pl, prime = [3,5,7], primegen()
+    for p in pl: next(prime)
+    n = next(prime); nn = n*n
+    while True:
+        n = next(prime); ll, nn = nn, n*n
+        delta = nn - ll
+        sieve = bytearray([True]) * delta
+        for p in pl:
+            k = (-ll) % p
+            sieve[k::p] = bytearray([False]) * ((delta-k)//p + 1)
+        if nn > limit: break
+        yield from compress(range(ll,ll+delta,2), sieve[::2])
+        pl.append(n)
+    yield from takewhile(ltlim, compress(range(ll,ll+delta,2), sieve[::2]))
+
+

imp/math/util.py

@@ -13,3 +13,10 @@ def clamp_min(n: int, max: int) -> int:
 
 def digits(n: int) -> int:
     return len(str(n))
+
+# NOTE: assumes A and B are equal length
+def xor_bytes(A: bytes, B: bytes) -> bytes:
+    return b''.join([(a ^ b).to_bytes() for (a, b) in zip(A, B)])
+def xor_str(A: str, B: str) -> str:
+    return ''.join([chr(ord(a) ^ ord(b)) for (a, b) in zip(A, B)])
+