diff --git a/imp/crypto/cyphers.py b/imp/crypto/cyphers.py index b4b1297..b6896eb 100644 --- a/imp/crypto/cyphers.py +++ b/imp/crypto/cyphers.py @@ -1,32 +1,137 @@ +''' +Implementation of various substitution cyphers, some of which are +already in the standard library but lack features/options I require. + +Terminology: +1. "Simple": simple substitution cyphers operates on single letters +2. "Polygraphic": polygraphic substitution cyphers operate on groups + of letters larger than length 1. +3. "Monoalphabetic cypher": a monoalphabetic cypher applies the same fixed + substitutions to the entire message (ie ROT13) +4. "Polyalphabetic cypher": a polyalphabetic cypher applies different + substitutions to the entire message (ie vigenere) +''' + from imp.constants import ALPHA_LOWER, ALPHA_UPPER -ERRMSG_ROT_CHARSET_DUPL = 'Bad charsets for ROT: found duplicate character \'{0}\'' -ERRMSG_ROT_CHAR_UNKNOWN = 'No charset for ROT has character \'{0}\' (exception thrown since `ignore_noncharset=False`)' +''' +Constant Declarations +''' +# Error Message Format Strings +ERRMSG_CHARSET_DUPL = 'Bad charsets for {0}: found duplicate character \'{1}\'' +ERRMSG_NOT_IN_CHARSET = 'No charset for {0} has character \'{1}\'' +# Charsets +CHARSET_SUB_DEFAULT = [ALPHA_LOWER, ALPHA_UPPER] ''' -Substitution Cyphers +========================================== +Simple Monoalphabetic Substitution Cyphers +========================================== ''' -def ROT(plaintext: str, +def ROT(plaintext: str | list, n: int, - charsets: list[str] = [ALPHA_LOWER, ALPHA_UPPER], - ignore_noncharset: bool = True) -> str: - cyphertext = '' + charsets: list[str] | list[list] = CHARSET_SUB_DEFAULT, + ignore_noncharset: bool = True, + as_string: bool = True) -> str: + cyphertext = [] for c in plaintext: - index = None + i = None + active_charset = None for charset in charsets: try: i = charset.index(c) # if we reached here then no ValueError occured - if index is not None: - raise ValueError(ERRMSG_ROT_CHARSET_DUPL.format(c)) - break + if active_charset is not None: + raise ValueError(ERRMSG_CHARSET_DUPL.format('ROT', c)) + active_charset = charset except ValueError: pass - if index is not None: - cyphertext += charset[(i + n) % len(charset)] + + if i is not None: + cyphertext.append(active_charset[(i + n) % len(active_charset)]) elif not ignore_noncharset: - raise ValueError(ERRMSG_ROT_CHAR_UNKNOWN.format(c)) + raise ValueError(ERRMSG_NOT_IN_CHARSET.format('ROT', c)) + else: + # if ignore_noncharset then just append c + cyphertext.append(c) + if as_string: + return ''.join(cyphertext) return cyphertext # Common ROT aliases -def ROT13(plaintext: str) -> str: return ROT(plaintext, 13) -def caesar(plaintext: str) -> str: return ROT(plaintext, 1) +def ROT13(plaintext: str, charsets: list[str] = CHARSET_SUB_DEFAULT) -> str: + return ROT(plaintext, 13, charsets=charsets) +def caesar(plaintext: str, charsets: list[str] = CHARSET_SUB_DEFAULT) -> str: + return ROT(plaintext, 1, charsets=charsets) + +# The Atbash Cypher maps a character to its reverse (ie a -> z, b -> y, etc) +def atbash(plaintext: str | list, + charsets: list[str] | list[list] = CHARSET_SUB_DEFAULT, + ignore_noncharset: bool = True, + as_string: bool = True) -> str: + cyphertext = [] + for c in plaintext: + i = None + active_charset = None + for charset in charsets: + try: + i = charset.index(c) + # if we reached here then no ValueError occured + if active_charset is not None: + raise ValueError(ERRMSG_CHARSET_DUPL.format('Atbash', c)) + active_charset = charset + except ValueError: pass + + if i is not None: + L = len(active_charset) + cyphertext.append(active_charset[(L - i - 1) % L]) + elif not ignore_noncharset: + raise ValueError(ERRMSG_NOT_IN_CHARSET.format('Atbash', c)) + else: + # if ignore_noncharset then just append c + cyphertext.append(c) + if as_string: + return ''.join(cyphertext) + return cyphertext + +''' +========================================= +Simple Polyalphabetic Substituion Cyphers +========================================= +''' +def vigenere(plaintext: str | list, + key: str | list, + charsets: list[str] | list[list] = CHARSET_SUB_DEFAULT, + key_charset: str | list = ALPHA_UPPER, + decrypt: bool = False, + ignore_noncharset: bool = True, + as_string: bool = True) -> str: + cyphertext = [] + # map the key characters to their rotations + rots = [key_charset.index(k) for k in key] + pos = 0 + for c in plaintext: + i = None + active_charset = None + for charset in charsets: + try: + i = charset.index(c) + # if we reached here then no ValueError occured + if active_charset is not None: + raise ValueError(ERRMSG_CHARSET_DUPL.format('Vigenere', c)) + active_charset = charset + except ValueError: pass + + if i is not None: + rot = rots[pos % len(rots)] + if decrypt: rot = -rot + cyphertext.append(active_charset[(i + rot) % len(active_charset)]) + pos += 1 + elif not ignore_noncharset: + raise ValueError(ERRMSG_NOT_IN_CHARSET.format('Vigenere', c)) + else: + # if ignore_noncharset then just append c + cyphertext.append(c) + if as_string: + return ''.join(cyphertext) + return cyphertext + diff --git a/imp/extern/__init__.py b/imp/extern/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/imp/extern/primefac.py b/imp/extern/primefac.py new file mode 100644 index 0000000..a6320b1 --- /dev/null +++ b/imp/extern/primefac.py @@ -0,0 +1,987 @@ +''' +Credit: Lucas A. Brown (@lucasaugustus) for his primefac module (v2.0.12): + https://github.com/lucasaugustus/primefac +primefac is under an MIT license. Moreover I've attached Lucas' credits statement: + "Significant parts of this code are derived or outright copied from other + people's work. In particular, the SIQS code was derived mostly verbatim + from https://github.com/skollmann/PyFactorise by Stephan Kollmann, while + the functions to manipulate points on elliptic curves were copied from a + reply to the blog post at http://programmingpraxis.com/2010/04/23/. The + rest, I believe, is my own work, but I may have forgotten something." + - Lucas A. Brown (primefac.py) +''' + +from math import sqrt, isqrt, log2, ceil, gcd, inf, prod as iterprod +from itertools import takewhile, compress, count +from multiprocessing import Process, Queue as mpQueue +from random import randrange#, seed; seed(0) +from time import time +from datetime import datetime, timedelta + +def primegen(limit=inf): + """ + Generates primes < limit almost lazily by a segmented sieve of Eratosthenes. + Memory usage depends on the sequence of prime gaps. Unconditionally, we can + say that memory usage is O(p^0.7625), where p is the most-recently-yielded + prime. On the Riemann Hypothesis and Cramer's conjecture, we can bring this + down to O(p^0.75 * log(p)) and O(sqrt(p) * log(p)^2), respectively. + Input: limit -- a number (default = inf) + Output: sequence of integers + Examples: + >>> list(islice(primegen(), 20)) + [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71] + >>> list(primegen(73)) + [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71] + """ + # We don't sieve 2, so we ought to be able to get sigificant savings by halving the length of the sieve. + # But the tiny extra computation involved in that seems to exceed the savings. + yield from takewhile(lambda x: x < limit, (2,3,5,7,11,13,17,19,23,29,31,37,41,43,47)) + pl, pg = [3,5,7], primegen() + for p in pl: next(pg) + n = next(pg); nn = n*n + while True: + n = next(pg) + ll, nn = nn, n*n + sl = (nn - ll) + sieve = bytearray([True]) * sl + for p in pl: + k = (-ll) % p + sieve[k::p] = bytearray([False]) * ((sl-k)//p + 1) + if nn > limit: break # TODO bring this condition up to the while statement + yield from compress(range(ll,ll+sl,2), sieve[::2]) + pl.append(n) + yield from takewhile(lambda x: x < limit, compress(range(ll,ll+sl,2), sieve[::2])) + +def ilog(x, b): # TODO: investigate optimization starting from x.bin_length() * 2 // b + """ + Greatest integer l such that b**l <= x + Input: x, b -- integers + Output: An integer + Examples: + >>> ilog(263789, 10) + 5 + >>> ilog(1023, 2) + 9 + """ + l = 0 + while x >= b: + x //= b + l += 1 + return l + # TODO possible optimization route: x.bit_length() == ilog(x, 2) + 1; we can therefore use x.bit_length() * 2 // b as a + # 1st approximation to ilog(x, b), then compute pow(b, x.bit_length() * 2 // b), then compare that to x and adjust. + +def modinv(a, m): + """ + Returns the inverse of a modulo m, normalized to lie between 0 and m-1. If + a is not coprime to m, return None. + Input: + a -- an integer coprime to m + n -- a positive integer + Output: None or an integer x between 0 and m-1 such that (a * x) % m == 1 + Examples: + >>> [modinv(1,1), modinv(2,5), modinv(5,8), modinv(37,100), modinv(12,30)] + [0, 3, 5, 73, None] + """ + return pow(a, -1, m) + +def introot(n, r=2): # TODO Newton iteration? + """ + For returns the rth root of n, rounded to the nearest integer in the + direction of zero. Returns None if r is even and n is negative. + Input: + n -- an integer + r -- a natural number or None + Output: An integer + Examples: + >>> [introot(-729, 3), introot(-728, 3), introot(1023, 2), introot(1024, 2)] + [-9, -8, 31, 32] + """ + if n < 0: return None if r%2 == 0 else -introot(-n, r) + if n < 2: return n + if r == 1: return n + if r == 2: return isqrt(n) + #if r % 2 == 0: return introot(isqrt(n), r//2) # TODO Check validity of this line. + lower = upper = 1 << (n.bit_length() // r) + while lower ** r > n: lower >>= 2 + while upper ** r <= n: upper <<= 2 + while lower != upper - 1: + mid = (lower + upper) // 2 + m = mid**r + if m == n: return mid + elif m < n: lower = mid + elif m > n: upper = mid + return lower + +def ispower(n, r=0): + """ + If r == 0: + If n is a perfect power, return a tuple containing largest integer (in + terms of magnitude) that, when squared/cubed/etc, yields n as the first + component and the relevant power as the second component. + If n is not a perfect power, return None. + If r > 0: + We check whether n is a perfect rth power; we return its rth root if it + is and None if it isn't. + Input: + n -- an integer + r -- an integer + Output: An integer, a 2-tuple of integers, or None + Examples: + >>> [ispower(n) for n in [64, 25, -729, 1729]] + [(8, 2), (5, 2), (-9, 3), None] + >>> [ispower(64, r) for r in range(7)] + [(8, 2), 64, 8, 4, None, None, 2] + """ + #if r == 0: return any(ispower(n, r) for r in primegen(n.bit_length()+1)) + #return n == introot(n, r) ** r + if r == 0: + if n in (0, 1, -1): return (n, 1) + for r in primegen(n.bit_length()+1): + x = ispower(n, r) + if x is not None: return (x, r) + return None + # TODO tricks for special cases + if (r == 2) and (n & 2): return None + if (r == 3) and (n & 7) in (2,4,6): return None + x = introot(n, r) + return None if x is None else (x if x**r == n else None) + +def jacobi(a, n): + """ + The Jacobi symbol (a|n). + Input: + a -- any integer + n -- odd integer + Output: -1, 0, or 1 + Examples: + >>> [jacobi(a, 15) for a in [-10, -7, -4, -2, -1, 0, 1, 2, 4, 7, 10]] + [0, 1, -1, -1, -1, 0, 1, 1, 1, -1, 0] + >>> [jacobi(a, 13) for a in [-10, -9, -4, -2, -1, 0, 1, 2, 4, 9, 10]] + [1, 1, 1, -1, 1, 0, 1, -1, 1, 1, 1] + >>> [jacobi(a, 11) for a in [-10, -9, -4, -2, -1, 0, 1, 2, 4, 9, 10]] + [1, -1, -1, 1, -1, 0, 1, -1, 1, 1, -1] + """ + if (n%2 == 0) or (n < 0): return None # n must be a positive odd number TODO delete this check? + if (a == 0) or (a == 1): return a + a, t = a%n, 1 + while a != 0: + while not a & 1: + a //= 2 + if n & 7 in (3, 5): t *= -1 + a, n = n, a + if (a & 3 == 3) and (n & 3) == 3: t *= -1 + a %= n + return t if n == 1 else 0 + +def isprime(n, tb=(3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59)): # TODO optimize the basis, possibly by varying it with n + """ + BPSW primality test variant: we use the strong Lucas PRP test and preface + the computation with trial division for speed. No composites are known to + pass the test, though it is suspected that infinitely many will do so. + There are definitely no such errors below 2^64. + This function is mainly a streamlined version of bpsw(). + Input: + n -- integer. Number to be examined. + tb -- iterable of primes. Basis for trial division. + Output: True if probably prime; False if definitely composite. + Examples: + >>> [n for n in range(91) if isprime(1000*n+1)] + [3, 4, 7, 9, 13, 16, 19, 21, 24, 28, 51, 54, 55, 61, 69, 70, 76, 81, 88, 90] + >>> [isprime(rpn("38 ! 1 +")), isprime(rpn("38 ! 1 -"))] + [False, True] + """ + # 1. Do some trial division with tb as the basis. + if n % 2 == 0 or n < 3: return n == 2 + for p in tb: + if n % p == 0: return n == p + + # 2. If sprp(n,2) fails, return False. If it succeeds, continue. + t, s = (n - 1) // 2, 1 + while t % 2 == 0: t //= 2; s += 1 + #assert 1 + 2**s * t == n + x = pow(2, t, n) + if x != 1 and x != n - 1: + for j in range(1, s): + x = pow(x, 2, n) + if x == 1: return False + elif x == n - 1: break + else: return False + + # 3. Select parameters for slprp. + for D in count(5, 4): + j = jacobi(D, n) + if j == 0: return D == n + if j == -1: break + D = -2 - D + j = jacobi(D, n) + if j == 0: return -D == n + if j == -1: break + if D == -13 and ispower(n,2): return False # If n is square, then this loop amounts to very slow trial division. + + # Now run slprp(n, 1, (1 - D) // 4) and return the result. + b = (1 - D) // 4 + if 1 < gcd(n, b) < n: return False + s, t = 1, (n + 1) // 2 + while t % 2 == 0: s += 1; t //= 2 + v, w, q, Q = 2, 1, 1, 1 + for k in bin(t)[2:]: + q = (q*Q) % n + if k == '1': Q, v, w = (q*b) % n, (w*v - q) % n, (w*w - 2*q*b) % n + else: Q, w, v = q , (w*v - q) % n, (v*v - 2*q ) % n + # assert ( (2*w-v) * modinv(D,n) ) % n, v == lucasmod(t, 1, b, n) + if v == 0 or ( (2*w-v) * modinv(D,n) ) % n == 0: return True + q = pow(b, t, n) + for _ in range(1, s): + v = (v*v - 2*q) % n + if v == 0: return True + q = (q*q) % n + return False + +def pollardrho_brent(n, verbose=False): + """ + Factors integers using Brent's variation of Pollard's rho algorithm. + If n is prime, we immediately return n; if not, we keep chugging until a + nontrivial factor is found. This function calls the randomizer; two + successive calls may therefore return two different results. + Input: n -- number to factor + Output: A factor of n. + Examples: + >>> n = rpn("20 ! 1 +"); f = pollardrho_brent(n); n % f + 0 + """ + if isprime(n): return n + g = n + while g == n: + y, c, m, g, r, q = randrange(1, n), randrange(1, n), randrange(1, n), 1, 1, 1 + while g==1: + x, k = y, 0 + for i in range(r): y = (y**2 + c) % n + while k < r and g == 1: + ys = y + for i in range(min(m, r-k)): + y = (y**2 + c) % n + q = q * abs(x-y) % n + g, k = gcd(q, n), k+m + r *= 2 + if g==n: + while True: + ys = (ys**2+c)%n + g = gcd(x-ys, n) + if g > 1: break + return g + +def pollard_pm1(n, B1=100, B2=1000, verbose=False): # TODO: What are the best default bounds and way to increment them? + """ + Integer factoring function. Uses Pollard's p-1 algorithm. Note that this + is only efficient if the number to be factored has a prime factor p such + that p-1's largest prime factor is "small". In this implementation, that + tends to mean less than 10,000,000 or so. + Input: + n -- number to factor + B1 -- Natural number. Bound for phase 1. Default == 100. + B2 -- Natural number > B1. Bound for phase 2. Default == 1000. + Output: A factor of n. + Examples: + >>> pollard_pm1(rpn("28 ! 1 - 239 //")) + 1224040923709997 + """ + if isprime(n): return n + m = ispower(n) + if m: return m[0] + while True: + pg = primegen() + q = 2 # TODO: what about other initial values of q? + p = next(pg) + while p <= B1: q, p = pow(q, p**ilog(B1, p), n), next(pg) + g = gcd(q-1, n) + if 1 < g < n: return g + while p <= B2: q, p = pow(q, p, n), next(pg) + g = gcd(q-1, n) + if 1 < g < n: return g + # These bounds failed. Increase and try again. + B1 *= 10 + B2 *= 10 + +def mlucas(v, a, n): + # Helper for williams_pp1(). Multiplies along a Lucas sequence mod n. + v1, v2 = v, (v**2 - 2) % n + for bit in bin(a)[3:]: v1, v2 = ((v1**2 - 2) % n, (v1*v2 - v) % n) if bit == "0" else ((v1*v2 - v) % n, (v2**2 - 2) % n) + return v1 +def williams_pp1(n, verbose=False): # TODO: experiment with different values of v0, and implement the two-phase version + """ + Integer factoring function. Uses Williams' p+1 algorithm, single-stage + variant. Note that this is only efficient when the number to be factored + has a prime factor p such that p+1's largest prime factor is "small". + Input: n -- integer to factor + Output: Integer. A nontrivial factor of n. + Example: + >>> williams_pp1(315951348188966255352482641444979927) + 12403590655726899403 + """ + if isprime(n): return n + m = ispower(n) + if m: return m[0] + for v in count(3): + for p in primegen(): + e = ilog(isqrt(n), p) + if e == 0: break + for _ in range(e): v = mlucas(v, p, n) + g = gcd(v - 2, n) + if 1 < g < n: return g + if g == n: break + +def ecadd(p1, p2, p0, n): + # Helper for ecm(). Adds two points on a Montgomery curve mod n. + x1,z1 = p1; x2,z2 = p2; x0,z0 = p0 + t1, t2 = (x1-z1)*(x2+z2), (x1+z1)*(x2-z2) + return (z0*pow(t1+t2,2,n) % n, x0*pow(t1-t2,2,n) % n) +def ecdub(p, A, n): + # Helper for ecm(). Doubles a point on a Montgomery curve mod n. + x, z = p; An, Ad = A + t1, t2 = pow(x+z,2,n), pow(x-z,2,n) + t = t1 - t2 + return (t1*t2*4*Ad % n, (4*Ad*t2 + t*An)*t % n) +def ecmul(m, p, A, n): + # Helper for ecm(). Multiplies a point on a Montgomery curve mod n. + if m == 0: return (0, 0) + elif m == 1: return p + else: + q = ecdub(p, A, n) + if m == 2: return q + b = 1 + while b < m: b *= 2 + b //= 4 + r = p + while b: + if m&b: q, r = ecdub(q, A, n), ecadd(q, r, p, n) + else: q, r = ecadd(r, q, p, n), ecdub(r, A, n) + b //= 2 + return r +def secm(n, B1, B2, seed): + """ + Seeded ECM. Helper function for ecm(). Returns a possibly-trivial divisor + of n given two bounds and a seed. Uses the two-phase algorithm on + Montgomery curves. See https://wp.me/prTJ7-zI and https://wp.me/prTJ7-A7 + for more details. Most of the code for this function's "helpers" were + shamelessly copied from the first of those links. + Input: + n -- Integer to factor + B1 -- Integer. Number of iterations for the first phase. + B2 -- Integer. Number of iterations for the second phase. + seed -- Integer. Selects the specific curve we'll be working on. + Output: Integer. A possibly-trivial factor of n. + Examples: + >>> secm(rpn('24 ! 1 -'), 100, 1000, 22) + 991459181683 + """ + u, v = (seed**2 - 5) % n, 4*seed % n + p = pow(u, 3, n) + Q, C = (pow(v-u,3,n)*(3*u+v) % n, 4*p*v % n), (p, pow(v,3,n)) + pg = primegen() + p = next(pg) + while p <= B1: Q, p = ecmul(p**ilog(B1, p), Q, C, n), next(pg) + g = gcd(Q[1], n) + if 1 < g < n: return g + while p <= B2: + # There is a trick that can speed up the second stage. Instead of multiplying each prime by Q, we iterate over i from + # B1 + 1 to B2, adding 2Q at each step; when i is prime, the current Q can be accumulated into the running solution. + # Again, we defer the calculation of the GCD until the end of the iteration. TODO: Implement and compare performance. + Q = ecmul(p, Q, C, n) + g *= Q[1] + g %= n + p = next(pg) + return gcd(g, n) +def ecmparams(n): # TODO: Better parameters. + counter = 0 + for i in count(): + for _ in range(2*i+1): + yield (2**i, 10 * 2**i, randrange(6,n), counter) + counter += 1 + for j in range(i+1): + yield (2**j, 10 * 2**j, randrange(6,n), counter) + counter += 1 +def ecm(n, paramseq=ecmparams, nprocs=1, verbose=False): + """ + "Modern" integer factoring via elliptic curves. Uses Montgomery curves, the + two-phase algorithm, and (optionally) multiple processes. The hard work is + done by secm(); this function just does the managerial work of pulling a + sequence of parameters out of a generator and feeding them into secm(). + Input: + n -- number to factor + paramseq -- sequence of parameters to feed into secm(). It must be an + infinite generator of 4-tuples (a,b,c,d), where a is the + number of iterations for the first phase, b is the number of + iterations for the second phase, c is a seed to select the + curve to work on, and d is an auxiliary used to count the + parameters generated so far. We need a < b and 6 <= c < n. + nprocs -- number of processes to use. Default == 1. Setting this > 1 + is discouraged on "small" inputs because managing multiple + processes incurs significant overhead. + Output: + A factor of n. + Note that if the parameter sequence calls the randomizer (which is + currently the default behavior), then two successive calls may therefore + return two different results. + Examples: + >>> n = 625793187653 * 991459181683 # = 620448401733239439359999 = 24! - 1 + >>> f = ecm(n) + >>> (n//f) * f + 620448401733239439359999 + """ + g = n % 6 + if g % 2 == 0: return 2 + if g % 3 == 0: return 3 + if isprime(n): return n + m = ispower(n) + if m: return m[0] + if nprocs == 1: + for (B1,B2,seed,i) in paramseq(n): + f = secm(n, B1, B2, seed) + if 1 < f < n: return f + assert nprocs != 1 + def factory(params, output): output.put(secm(*params)) + ps, facs, procs = paramseq(n), mpQueue(), [] + procs = [Process(target=factory, args=((n,)+next(ps)[:3], facs)) for _ in range(nprocs)] + for p in procs: p.start() + while True: + g = facs.get() + if 1 < g < n: + for p in procs: p.terminate() + return g + for p in range(nprocs): # TODO: Try doing this with a process pool. + if not procs[p].is_alive(): + del procs[p] + break + procs.append(Process(target=factory, args=((n,)+next(ps)[:3], facs))) + procs[-1].start() + +def sqrtmod_prime(a, p): + """ + Solves x**2 == a (mod p) for x. We assume that p is a prime and a is a + quadratic residue modulo p. If either of these assumptions is false, the + return value is meaningless. + The Cipolla-Lehmer section is my own. The rest appears to be derived from + https://codegolf.stackexchange.com/a/9088. + Input: + a -- natural number + p -- prime number + Output: whole number less than p + Examples: + >>> sqrtmod_prime(4, 5) + 3 + >>> sqrtmod_prime(13, 23) + 6 + >>> sqrtmod_prime(997, 7304723089) + 761044645 + """ + a %= p + if p%4 == 3: return pow(a, (p+1) >> 2, p) + elif p%8 == 5: + v = pow(a << 1, (p-5) >> 3, p) + return (a*v*(((a*v*v<<1)%p)-1))%p + elif p%8 == 1: + # CranPom ex 2.31, pg 112 / 126. Pretty sure this amounts to Cipolla-Lehmer. + if a == 0: return 0 # Necessary to avoid an infinite loop in the legendre section + h = 2 + # while legendre(h*h - 4*a, p) != -1: + while pow(h*h - 4*a, (p-1) >> 1, p) != p - 1: h += 1 # TODO compare speed vs random selection + #return ( lucasmod((p+1)//2, h, a, p)[1] * modinv(2, p) ) % p + k, v, w, q, Q = (p+1)//2, 2, h % p, 1, 1 + for kj in bin(k)[2:]: + q = (q*Q) % p + if kj == '1': Q, v, w = (q*a) % p, (w*v - h*q) % p, (w*w - 2*q*a) % p + else: Q, w, v = q , (w*v - h*q) % p, (v*v - 2*q ) % p + return (v*k) % p + else: return a # p == 2 + +def siqs(n, verbose=False): + """ + Uses the Self-Initializing Quadratic Sieve to extract a factor of n. + We make some calls to randrange, so the output may change from call to call. + This is derived from https://github.com/skollmann/PyFactorise. + Input: + n -- number to factor + verbose -- if True, print progress reports. Default == False. + Output: + If n is prime, we return n. + If n is composte, we return a nontrivial factor of n. + If n is too small, we raise a ValueError. + Examples: + >>> siqs(factorial(24) - 1) in (625793187653, 991459181683) + True + """ + + if (not isinstance(n, int)) or n < 0: raise ValueError("Number must be a positive integer.") + + if n < 2**64: + if verbose: print("Number is too small for SIQS. Using Pollard Rho instead.") + return pollardrho_brent(n) + + if verbose: print("Factorizing %d (%d digits)..." % (n, len(str(n)))) + + if isprime(n): return n + + if verbose: print("Number is composite.") + if verbose: print("Checking whether it is a perfect power...") + perfect_power = ispower(n) + if perfect_power: + if verbose: print(n, "=", "%d^%d." % (perfect_power[0], perfect_power[1])) + return perfect_power[0] + + if verbose: print("Not a perfect power.") + if verbose: print("Using SIQS on %d (%d digits)..." % (n, len(str(n)))) + + if verbose: starttime, sievetime, latime = time(), 0, 0 + + # Choose parameters nf (sieve of factor base) and m (for sieving in [-m,m]. + # Using similar parameters as msieve-1.52 + dig = len(str(n)) + if dig <= 34: nf, m = 200, 65536 + elif dig <= 36: nf, m = 300, 65536 + elif dig <= 38: nf, m = 400, 65536 + elif dig <= 40: nf, m = 500, 65536 + elif dig <= 42: nf, m = 600, 65536 + elif dig <= 44: nf, m = 700, 65536 + elif dig <= 48: nf, m = 1000, 65536 + elif dig <= 52: nf, m = 1200, 65536 + elif dig <= 56: nf, m = 2000, 65536 * 3 + elif dig <= 60: nf, m = 4000, 65536 * 3 + elif dig <= 66: nf, m = 6000, 65536 * 3 + elif dig <= 74: nf, m = 10000, 65536 * 3 + elif dig <= 80: nf, m = 30000, 65536 * 3 + elif dig <= 88: nf, m = 50000, 65536 * 3 + elif dig <= 94: nf, m = 60000, 65536 * 9 + else: nf, m = 100000, 65536 * 9 + + class FactorBasePrime: + """A factor base prime for the SIQS""" + def __init__(self, p, tmem, lp): + self.p, self.soln1, self.soln2, self.tmem, self.lp, self.ainv = p, None, None, tmem, lp, None + + # Compute and return nf factor base primes suitable for a SIQS on n. + factor_base = [] + for p in primegen(): + if pow(n, (p-1) >> 1, p) == 1: # if n is a quadratic residue mod p + t = sqrtmod_prime(n, p) + lp = round(log2(p)) # This gets rounded for the sake of speed. + factor_base.append(FactorBasePrime(p, t, lp)) + if len(factor_base) >= nf: break + + if verbose: print("Factor base size: %d primes.\nLargest prime in base: %d." % (nf, factor_base[-1].p)) + + npolys, i_poly, prev_cnt, relcount, smooth_relations, required_relations_ratio = 0, 0, 0, 0, [], 1.05 + + p_min_i, p_max_i = None, None + for (i,fb) in enumerate(factor_base): + if p_min_i is None and fb.p >= 400: # 400 is a tunable parameter. + p_min_i = i + if p_max_i is None and fb.p > 4000: # 4000 is a tunable parameter. + p_max_i = i - 1 + break + + # The following may happen if the factor base is small, make sure that we have enough primes. + if p_max_i is None: p_max_i = len(factor_base) - 1 + if p_min_i is None or p_max_i - p_min_i < 20: p_min_i = min(p_min_i, 5) # TODO This line is problematic for some small n. + + target = sqrt(2 * float(n)) / m + target1 = target / sqrt((factor_base[p_min_i].p + factor_base[p_max_i].p) / 2) + + while True: + + if verbose: temptime = time() + + if verbose: print("*** Phase 1: Finding smooth relations ***") + required_relations = round(nf * required_relations_ratio) + if verbose: print("Target: %d relations" % required_relations) + enough_relations = False + while not enough_relations: + + if i_poly == 0: # Compute the first of a set of polynomials + # find q such that the product of factor_base[q_i] is about sqrt(2n)/m; try a few sets to find a good one + best_q, best_a, best_ratio = None, None, None + for _ in range(30): + a, q = 1, [] + + while a < target1: + p_i = 0 + while p_i == 0 or p_i in q: p_i = randrange(p_min_i, p_max_i+1) + a *= factor_base[p_i].p + q.append(p_i) + + ratio = a / target + + # ratio too small seems to be not good + if (best_ratio is None or (ratio >= 0.9 and ratio < best_ratio) or best_ratio < 0.9 and ratio > best_ratio): + best_q, best_a, best_ratio = q, a, ratio + a, q = best_a, best_q + + s, B = len(q), [] + for l in range(s): + fb_l = factor_base[q[l]] + q_l = fb_l.p + assert a % q_l == 0 + gamma = (fb_l.tmem * modinv(a // q_l, q_l)) % q_l + if gamma > q_l // 2: gamma = q_l - gamma + B.append(a // q_l * gamma) + + b = sum(B) % a + b_orig = b + if (2 * b > a): b = a - b + + assert 0 < b and 2 * b <= a and (b * b - n) % a == 0 + + g1, g2, g3, ga, gb = b * b - n, 2 * a * b, a * a, a, b_orig + ha, hb = a, b + for fb in factor_base: + if a % fb.p != 0: + fb.ainv = modinv(a, fb.p) + fb.soln1 = (fb.ainv * (fb.tmem - b)) % fb.p + fb.soln2 = (fb.ainv * (-fb.tmem - b)) % fb.p + + else: # Compute the (i+1)-th polynomial, given that g is the i-th polynomial. + #v = lowest_set_bit(i) + 1 + #z = -1 if ceil(i / (2 ** v)) % 2 == 1 else 1 + z = -1 if ceil(i_poly / (1 + (i_poly ^ (i_poly-1)))) % 2 == 1 else 1 + #b = (g.b + 2 * z * B[v - 1]) % g.a + b = (gb + 2 * z * B[(i_poly & (-i_poly)).bit_length() - 1]) % ga + a = ga + b_orig = b + if (2 * b > a): b = a - b + assert (b * b - n) % a == 0 + + g1, g2, g3, ga, gb = b * b - n, 2 * a * b, a * a, a, b_orig + ha, hb = a, b + for fb in factor_base: + if a % fb.p != 0: + fb.soln1 = (fb.ainv * ( fb.tmem - b)) % fb.p + fb.soln2 = (fb.ainv * (-fb.tmem - b)) % fb.p + + + i_poly += 1 + npolys += 1 + if i_poly >= 2 ** (len(B) - 1): i_poly = 0 + # BEGIN SIEVING. Most of our time is spent between here and the "END SIEVING" comment. + sieve_array = [0] * (2 * m + 1) + for fb in factor_base: + if fb.soln1 is None: continue + p, lp = fb.p, fb.lp + a_start_1 = fb.soln1 - ((m + fb.soln1) // p) * p + if p > 20: + for a in range(a_start_1 + m, 2 * m + 1, p): sieve_array[a] += lp + a_start_2 = fb.soln2 - ((m + fb.soln2) // p) * p + for a in range(a_start_2 + m, 2 * m + 1, p): sieve_array[a] += lp + # Perform the trial division step of the SIQS. + limit = round(log2(m * sqrt(float(n)))) - 25 # 25 is a tunable parameter. The rounding is for speed. + for (i,sa) in enumerate(sieve_array): + if sa >= limit: + x = i - m + gx = (g3 * x + g2) * x + g1 + # Determine whether gx can be fully factorized into primes from the factor base. + # If so, store the indices of the factors from the factor base as divisors_idx. If not, set that to None. + a, divisors_idx = gx, [] + for (fbi,fb) in enumerate(factor_base): + if a % fb.p == 0: + exp = 0 + while a % fb.p == 0: + a //= fb.p + exp += 1 + divisors_idx.append((fbi, exp)) + if a == 1: + u = ha * x + hb + v = gx + assert (u * u) % n == v % n + smooth_relations.append((u, v, divisors_idx)) + break + relcount = len(smooth_relations) + if relcount >= required_relations: + enough_relations = True + break + # END SIEVING. Most of our time is spent between here and the "BEGIN SIEVING" comment. + + if verbose and relcount > 0 and (relcount >= required_relations or i_poly % 8 == 0 or relcount > prev_cnt): + frac = relcount / required_relations + t = time() - starttime # Time spent so far + ett = t / frac # Estimated total time + ettc = ett - t # Estimated time to completion + eta = datetime.isoformat(datetime.now() + timedelta(seconds=ettc), sep=' ', timespec='seconds') + print('\b'*256 + "Found %d rels (%.1f%%) on %d polys. ETA %ds (%s). " % (relcount, 100*frac, npolys, ettc, eta), end='', flush=True) + prev_cnt = relcount + + if verbose: + sievetime += time() - temptime + print("\nSieving time: %f seconds." % sievetime) + temptime = time() + + print("*** Phase 2: Linear Algebra ***") + print("Building matrix for linear algebra step...") + M = [0] * nf + mask = 1 + for sr in smooth_relations: + for (j,exp) in sr[2]: + if exp % 2: M[j] += mask + mask <<= 1 + + if verbose: print("Finding perfect squares using matrix and factors from perfect squares...") + # Perform the linear algebra step of the SIQS: do fast Gaussian elimination to determine pairs of perfect squares mod n. + # Use the optimisations described in [Çetin K. Koç and Sarath N. Arachchige. 'A Fast Algorithm for Gaussian Elimination + # over GF(2) and its Implementation on the GAPP.' Journal of Parallel and Distributed Computing 13.1 (1991): 118-122]. + if verbose: gaussstart = time() + row_is_marked = bytearray([False]) * relcount + pivots = [-1] * nf + for j in range(nf): + M_j = M[j] + i = (M_j & (-M_j)).bit_length() - 1 #i = -1 if M[j] == 0 else lowest_set_bit(M[j]) + # i is now the row of the first nonzero entry in column j, or -1 if no such row exists. + if i > -1: + pivots[j] = i + row_is_marked[i] = True + for k in range(nf): + if (M[k] >> i) & 1 and k != j: # test M[i][k] == 1 + M[k] ^= M_j # add column j to column k mod 2 + if verbose and ((nf - j) % 100 == 0 or nf == j + 1): + frac = (j+1) / nf + t = time() - gaussstart # Time spent so far + ett = t / frac # Estimated total time + ettc = ett - t # Estimated time to completion + eta = datetime.isoformat(datetime.now() + timedelta(seconds=ettc), sep=' ', timespec='seconds') + print('\b'*256 + "Gaussian elimination: %d/%d (%.1f%%). ETA %ds (%s). " % (j+1, nf, 100*frac, ettc, eta), end='', flush=True) + if verbose: print("\nGaussian elimination time: %f seconds" % (time() - gaussstart)) + attempts = 0 + for i in range(relcount): + if not row_is_marked[i]: + square_indices = [i] + for j in range(nf): + if (M[j] >> i) & 1: # test M[i][j] == 1 + square_indices.append(pivots[j]) + # Given the solution encoded by square_indices, try to find a factor of n, and return it. + attempts += 1 + if verbose: print('\b'*42 + "Attempt #%d" % attempts, end='', flush=True) + # Given on of the solutions returned by siqs_solve_matrix and the corresponding smooth relations, + # calculate the pair (sqrt1, sqrt2) such that sqrt1^2 = sqrt2^2 (mod n). + sqrt1, sqrt2 = 1, 1 + for idx in square_indices: + sqrt1 *= smooth_relations[idx][0] + sqrt2 *= smooth_relations[idx][1] + sqrt2 = isqrt(sqrt2) + assert (sqrt1 * sqrt1) % n == (sqrt2 * sqrt2) % n + factor = gcd(sqrt1 - sqrt2, n) + if 1 != factor != n: + if verbose: + print(" succeeded.") + latime += time() - temptime + print("Linear algebra time: %f seconds." % latime) + totaltime = time() - starttime + print("Total time: %f seconds." % (totaltime)) + return factor + + if verbose: + print('\b'*42 + "All %d attempts failed." % attempts) + latime += time() - temptime + print("Linear algebra time: %f seconds." % latime) + + print("Failed to find a solution. Finding more relations...") + required_relations_ratio += 0.05 + +def multifactor(n, methods=(pollardrho_brent, pollard_pm1, williams_pp1, ecm, siqs), verbose=False): + """ + Integer factoring function. Uses several methods in parallel. Waits for a + function to return, kills the rest, and reports. Note that two successive + calls may return different results depending on which method finishes first + and whether any methods call the randomizer. + Input: + n -- number to factor + methods -- list of functions to run. + Output: A factor of n. + Examples: + >>> methods = [pollardrho_brent, pollard_pm1, williams_pp1, ecm, siqs] + >>> n = rpn("24 ! 1 -"); f = multifactor(n, methods); n%f + 0 + """ + # Note that the multiprocing incurs relatively significant overhead. Only call this if n is proving difficult to factor. + def factory(method, n, verbose, output): output.put((method(n, verbose=verbose), str(method).split()[1])) + factors = mpQueue() + procs = [Process(target=factory, args=(m, n, verbose, factors)) for m in methods] + for p in procs: p.start() + (f, g) = factors.get() + for p in procs: p.terminate() + names = {"pollardrho_brent":"prb", "pollard_pm1":"p-1", "williams_pp1":"p+1"} + return (f, names.get(g, g)) + +def primefac(n, trial=1000, rho=42000, verbose=False, methods=(pollardrho_brent,)): + """ + Generates the prime factors of the input. Factors that appear x times are + yielded x times. + Input: + n -- the number to be factored + trial -- Trial division is performed up to this limit and no further. + We trial divide by 2 and 3 whether the user wants to or not. + Default == 1000. + rho -- How long to have Pollard's rho chug before declaring a cofactor + difficult. Default == 42,000 iterations. Floating-point inf is + an acceptable value. + methods -- Use these methods on difficult cofactors. If the tuple has + more than 1 element, we have multifactor handle it. Calling + multifactor has a high overhead, so when the tuple has a + single element, we call that function directly. The default + is (pollardrho_brent,). Each function f in methods must + accept a single number n as its argument. If n is prime, + f(n) must return n. If n is composite, f(n) must return a + number strictly between 1 and n that evenly divides n. + Giving up is not allowed. + Output: Prime factors of n + Factoring: + We use a three-stage factoring algorithm. + 1. Trial divide with all primes less than or equal to the specified + limit. We trial divide by 2 and 3 regardless of the limit. + 2. Run Pollard's rho algorithm on whatever remains. This algorithm + may split a number into two composite cofactors. Such cofactors + remain here until they survive the specified number of rounds of + the rho algorithm. + 3. Subject each remaining cofactor to a set of up to five factoring + methods in parallel: + Pollard's rho algorithm with Brent's improvement, + Pollard's p-1 method, + Williams' p+1 method, + the elliptic curve method, + and the self-initializing quadratic sieve. + Examples: + >>> list(primefac(1729)) + [7, 13, 19] + >>> list(sorted(primefac(rpn("24 ! 1 -")))) + [625793187653, 991459181683] + """ + # Obtains a complete factorization of n, yielding the prime factors as they are obtained. + # If the user explicitly specifies a splitting method, use that method. Otherwise, + # 1. Pull out small factors with trial division. + # 2. Do a few rounds of Pollard's Rho algorithm. + # 3. Launch multifactor on the remainder. Note that multifactor's multiprocessing incurs relatively significant overhead. + + if verbose: print("\nFactoring %d (%d digits):" % (n, len(str(n)))) + if n < 0: + if verbose: print("Trial division: -1") + yield -1; n *= -1 + if n < 2: return + if isprime(n): + if verbose: print("Number is prime.") + yield n + return + if verbose: print("Number is composite.") + + # Trial division + factors, nroot = [], isqrt(n) + for p in primegen(max(4,trial)+1): # Note that we check for 2 and 3 whether the user wants to or not. + if verbose: print('\b'*80 + "Trial division:", p, end='', flush=False) + if n%p == 0: + while n%p == 0: + if verbose: print('\b'*80 + "Trial division:", p) + yield p + n //= p + nroot = isqrt(n) + if p > nroot: + if n != 1: + if verbose: + laststr = "Trial division: " + str(p) + print('\b'*80 + str(n) + " " * len(laststr)) + yield n + return + if isprime(n): + if verbose: + laststr = "Trial division: " + str(p) + print('\b'*80 + str(n) + " " * len(laststr)) + yield n + return + + if verbose: print("\nTrial division finished.") + + # TODO: Fermat's method? + + # Pollard's rho + factors, difficult = [n], [] + while len(factors) != 0: + rhocount = 0 + n = factors.pop() + if verbose: print("Beginning Pollard Rho on %d." % n) + try: + g = n + while g == n: + x, c, g = randrange(1, n), randrange(1, n), 1 + y = x + while g==1: + if rhocount >= rho: raise Exception + rhocount += 1 + x = (x**2 + c) % n + y = (y**2 + c) % n + y = (y**2 + c) % n + g = gcd(x-y, n) + # We now have a nontrivial factor g of n. If we took too long to get here, we're actually at the except statement. + f = n // g + if verbose: print("Pollard Rho split %d into %d and %d." % (n, g, f)) + if isprime(g): + if verbose: print(g, "is prime.") + yield g + else: + if verbose: print(g, "is composite.") + factors.append(g) + if isprime(f): + if verbose: print(f, "is prime.") + yield f + else: + if verbose: print(f, "is composite.") + factors.append(f) + except Exception: + if verbose: print("Declaring", n, "difficult.") + difficult.append(n) # Factoring n took too long. We'll have multifactor chug on it. + + # TODO: P-1 by itself? + # TODO: ECM by itself? + + names = {"pollardrho_brent":"prb", "pollard_pm1":"p-1", "williams_pp1":"p+1", "ecm":"ecm", "siqs":"siqs"} + factors = difficult + if len(methods) == 1: + method = methods[0] + name = names[str(method).split()[1]] + while len(factors) != 0: + n = factors.pop() + if verbose: print("Beginning %s on %d." % (name, n)) + f = method(n, verbose=verbose) + g = n // f + if verbose: print("%s split %d into %d and %d." % (name, n, f, g)) + if isprime(f): + if verbose: print(f, "is prime.") + yield f + else: + if verbose: print(f, "is composite.") + factors.append(f) + if isprime(g): + if verbose: print(g, "is prime.") + yield g + else: + if verbose: print(g, "is composite.") + factors.append(g) + else: + while len(factors) != 0: + n = factors.pop() + if verbose: + #methodstring = ", ".join(names[str(m).split()[1]] for m in methods) + if len(methods) == 2: methodstring = " and ".join(names[str(m).split()[1]] for m in methods) + else: + assert len(methods) > 2 + methodstring = ", ".join(names[str(m).split()[1]] for m in methods[:-1]) + methodstring += ", and " + names[str(methods[-1]).split()[1]] + if verbose: print("Beginning %s on %d." % (methodstring, n)) + f, name = multifactor(n, methods=methods, verbose=verbose) + g = n // f + if verbose: print("%s split %d into %d and %d." % (name, n, f, g)) + if isprime(f): + if verbose: print(f, "is prime.") + yield f + else: + if verbose: print(f, "is composite.") + factors.append(f) + if isprime(g): + if verbose: print(g, "is prime.") + yield g + else: + if verbose: print(g, "is composite.") + factors.append(g) diff --git a/imp/math/numbers/__init__.py b/imp/math/numbers/__init__.py new file mode 100644 index 0000000..22a0d21 --- /dev/null +++ b/imp/math/numbers/__init__.py @@ -0,0 +1,82 @@ +''' +Terminology: + Although "divisor" and "factor" mean the same thing. + When Celeste discusses "divisors of n" it is implied to + mean "proper divisors of n + n itself", and "factors" are + the "prime proper divisors of n". +''' + +from imp.extern.primefac import primefac + +def factors(n: int) -> int: + pfactors: list[tuple[int, int]] = [] + # generate primes and progressively store them in pfactors + pfgen = primefac(n) + watching = next(pfgen) + mult = 1 + # ASSUMPTION: prime generation is (non-strict) monotone increasing + while True: + p = next(pfgen, None) + if p == watching: + mult += 1 + else: + pfactors.append((watching, mult)) + watching = p # reset + mult = 1 # reset + if p is None: + break + return pfactors + +def factors2divisors(pfactors: list[tuple[int, int]], + sorted: bool = True) -> list[int]: + ''' + Generates all divisors < n of an integer n given its prime factorisation. + Input: prime factorisation of n (excluding 1 and n, and duplicates) + in the typical form: list[(prime, multiplicity)] + ''' + divisors = [1] + for (prime, multiplicity) in pfactors: + extension = [] + for i in range(1, multiplicity+1): + term = prime**i + extension.extend(list([divisor*term for divisor in divisors])) + divisors.extend(extension) + if sorted: divisors.sort() + return divisors + +def factors2aliquots(pfactors: list[tuple[int, int]]) -> list[int]: + return factors2divisors(pfactors)[:-1] + +# "aliquots(n)" is an alias for "divisors(n)" +def aliquots(n: int) -> int: + ''' + Returns all aliquot parts (proper divisors) of + an integer n, that is all divisors 0 < d <= n. + ''' + return factors2aliquots(factors(n)) + +def divisors(n: int) -> int: + ''' + Returns all divisors 0 < d < n of an integer n. + ''' + return factors2divisors(factors(n)) + +def aliquot_sum(n: int) -> int: + return sum(aliquots(n)) + + +def littleomega(n: int) -> int: + ''' + The Little Omega function counts the number of + distinct prime factors of an integer n. + Ref: https://en.wikipedia.org/wiki/Prime_omega_function + ''' + return len(factors(n)) + +def bigomega(n: int) -> int: + ''' + The Big Omega function counts the total number of + prime factors (including multiplicity) of an integer n. + Ref: https://en.wikipedia.org/wiki/Prime_omega_function + ''' + return sum(factor[1] for factor in factors(n)) diff --git a/imp/math/numbers/functions.py b/imp/math/numbers/functions.py new file mode 100644 index 0000000..8addd51 --- /dev/null +++ b/imp/math/numbers/functions.py @@ -0,0 +1,3 @@ +def factorial(n: int) -> int: + if n == 0: return 1 + return n * factorial(n-1) diff --git a/imp/math/numbers/kinds.py b/imp/math/numbers/kinds.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/imp/math/numbers/kinds.py @@ -0,0 +1 @@ + diff --git a/imp/math/primes.py b/imp/math/primes.py new file mode 100644 index 0000000..b130a0c --- /dev/null +++ b/imp/math/primes.py @@ -0,0 +1,47 @@ +from math import gcd, inf + +from imp.math.numbers import bigomega, factors +from imp.extern.primefac import ( + isprime, + primegen as Primes, +) + +def coprime(n: int, m: int) -> bool: + return gcd(n, m) == 1 + +def almostprime(n: int, k: int) -> bool: + ''' + A natural n is "k-almost prime" if it has exactly + k prime factors (including multiplicity). + ''' + return (bigomega(n) == k) + +def semiprime(n: int) -> bool: + ''' + A semiprime number is one that is 2-almost prime. + Ref: https://en.wikipedia.org/wiki/Semiprime + ''' + return almostprime(n, 2) + +def eulertotient(x: int | list) -> int: + ''' + Evaluates Euler's Totient function. + Input: `x: int` is prime factorised by Lucas A. Brown's primefac.py + else `x: list` is assumed to the prime factorisation of `x: int` + ''' + pfactors = x if isinstance(x, list) else factors(n) + return prod((p-1)*(p**(e-1)) for (p, e) in pfactors) +# def eulertotient(n: int) -> int: +# ''' +# Uses trial division to compute +# Euler's Totient (Phi) Function. +# ''' +# 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 diff --git a/imp/math/util.py b/imp/math/util.py index 2012e71..d0d058b 100644 --- a/imp/math/util.py +++ b/imp/math/util.py @@ -1,3 +1,6 @@ +from collections.abc import Iterable +from itertools import chain, combinations + def clamp(n: int, min: int, max: int) -> int: if n < min: return min @@ -20,3 +23,7 @@ def xor_bytes(A: bytes, B: bytes) -> bytes: def xor_str(A: str, B: str) -> str: return ''.join([chr(ord(a) ^ ord(b)) for (a, b) in zip(A, B)]) + +def powerset(iterable: Iterable) -> Iterable: + s = list(iterable) + return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))