diff --git a/imp/math/factor/__init__.py b/imp/math/factor/__init__.py new file mode 100644 index 0000000..e4a5de0 --- /dev/null +++ b/imp/math/factor/__init__.py @@ -0,0 +1,2 @@ +from .pftrialdivision import trial_division +from .factordb import DBResult, FType, FCertainty diff --git a/imp/math/factor/factordb.py b/imp/math/factor/factordb.py new file mode 100644 index 0000000..b24382e --- /dev/null +++ b/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) + diff --git a/imp/math/factor/pftrialdivision.py b/imp/math/factor/pftrialdivision.py new file mode 100644 index 0000000..8e1b3d0 --- /dev/null +++ b/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 + diff --git a/imp/math/primes.py b/imp/math/primes.py deleted file mode 100644 index 513f61e..0000000 --- a/imp/math/primes.py +++ /dev/null @@ -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 - diff --git a/imp/math/primes/__init__.py b/imp/math/primes/__init__.py new file mode 100644 index 0000000..ed6d3ca --- /dev/null +++ b/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])) + + diff --git a/imp/math/util.py b/imp/math/util.py index bcb9e3a..2012e71 100644 --- a/imp/math/util.py +++ b/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)]) +