Source code for thompson.number_theory

r"""At the bottom level of the package hierarchy are various helper functions which implement standard number-theoretic algorithms.
 
.. testsetup::
	
	from random                 import randint
	from thompson.number_theory import *
"""

import fractions

from collections import namedtuple
from functools   import reduce
from operator    import mul
from numbers     import Number

__all__ = ['lcm', 'gcd', 'extended_gcd', 'solve_linear_diophantine', 'solve_linear_congruence',
	'divisors', 'prod']

BaseSolutionSet = namedtuple('BaseSolutionSet', 'base increment')
[docs]class SolutionSet(BaseSolutionSet): r"""Solution sets to the orbit sharing test :math:`u\phi^k = v`. These are either - empty (no such :math:`k`) exists; - a singleton; or - a linear sequence a + b*n. Create solution sets as follows: .. doctest:: :options: -ELLIPSIS >>> print(SolutionSet(5, 2)) {..., 3, 5, 7, 9, 11, 13, ...} >>> print(SolutionSet.singleton(4)) {4} >>> print(SolutionSet.empty_set()) {} >>> print(SolutionSet.the_integers()) {..., -1, 0, 1, 2, 3, 4, ...} """ __slots__ = () def __new__(cls, base, increment): if increment is not None: increment = abs(increment) self = super(SolutionSet, cls).__new__(cls, base, increment) return self @classmethod def singleton(cls, value): return cls(value, None) @classmethod def empty_set(cls): return cls(None, None) @classmethod def the_integers(cls): return cls(0, 1)
[docs] def is_sequence(self): """Returns True if this set contains more than one distinct element; otherwise returns False. >>> SolutionSet.empty_set().is_sequence() False >>> SolutionSet.singleton(2).is_sequence() False >>> SolutionSet(base=4, increment=3).is_sequence() True >>> SolutionSet.the_integers().is_sequence() True """ return self.base is not None and self.increment is not None
[docs] def is_singleton(self): """Returns True if this contains precisely one element, otherwise False. >>> SolutionSet.empty_set().is_singleton() False >>> SolutionSet.singleton(2).is_singleton() True >>> SolutionSet(base=4, increment=3).is_singleton() False >>> SolutionSet.the_integers().is_singleton() False """ return self.base is not None and self.increment is None
[docs] def is_empty(self): """Returns True if this set is empty, otherwise False. >>> SolutionSet.empty_set().is_empty() True >>> SolutionSet.singleton(2).is_empty() False >>> SolutionSet(base=4, increment=3).is_empty() False >>> SolutionSet.the_integers().is_empty() False """ return self.base is None
[docs] def is_the_integers(self): """Returns True if this set contains every integer; otherwise returns False. >>> SolutionSet.empty_set().is_the_integers() False >>> SolutionSet.singleton(2).is_the_integers() False >>> SolutionSet(base=4, increment=3).is_the_integers() False >>> SolutionSet.the_integers().is_the_integers() True """ return self.is_sequence() and abs(self.increment) == 1
[docs] def __contains__(self, other): #other in self """Returns true if this set contains an *other* number. >>> 1024 in SolutionSet.empty_set() False >>> 1024 in SolutionSet.singleton(128) False >>> 1024 in SolutionSet(0, 256) True >>> 1024 in SolutionSet.the_integers() True """ if not isinstance(other, Number): return NotImplemented if self.is_empty(): return False if self.is_singleton(): return other == self.base return self.base % self.increment == other % self.increment
[docs] def __and__(self, other): #self = self & other """The ``&`` operator (usually used for bitwise and) stands for intersection of sets. .. doctest:: :options: -ELLIPSIS >>> phi = SolutionSet.empty_set() >>> Z = SolutionSet.the_integers() >>> singleton = SolutionSet.singleton >>> print(phi & phi) {} >>> print(phi & singleton(1)) {} >>> print(phi & SolutionSet(2, 3)) {} >>> print(singleton(1) & singleton(1)) {1} >>> print(singleton(1) & singleton(2)) {} >>> print(singleton(8) & SolutionSet(4, 2)) {8} >>> print(SolutionSet(1, 3) & SolutionSet(2, 3)) {} >>> print(SolutionSet(1, 3) & SolutionSet(1, 2)) {..., -5, 1, 7, 13, 19, 25, ...} >>> print(SolutionSet(1, 18) & SolutionSet(5, 24)) {} >>> print(SolutionSet(1, 18) & SolutionSet(13, 24)) {..., -35, 37, 109, 181, 253, 325, ...} >>> print(SolutionSet(1, 3) & Z) {..., -2, 1, 4, 7, 10, 13, ...} >>> print(Z & Z) {..., -1, 0, 1, 2, 3, 4, ...} """ if not isinstance(other, SolutionSet): return NotImplemented if self.is_empty() or other.is_empty(): return SolutionSet.empty_set() if self.is_the_integers(): return other if other.is_the_integers(): return self if self.is_singleton(): if self.base in other: return self return SolutionSet.empty_set() if other.is_singleton(): if other.base in self: return other return SolutionSet.empty_set() #Solve s.base + s.inc * x = o.base + o.inc * y for (x,y) a = self.increment b = -other.increment c = other.base - self.base #Equation above is equiv to ax + by = c --- a linear Diophantine equation. result = solve_linear_diophantine(a, b, c) if result is None: return SolutionSet.empty_set() base, inc, lcm = result new_base = self.base + base[0] * self.increment new_base %= lcm assert new_base in self and new_base in other return SolutionSet(new_base, lcm)
def __str__(self): if self.is_empty(): return '{}' if self.is_singleton(): return '{{{0}}}'.format(self.base) values = (self.base + i * self.increment for i in range(-1, 5)) values = (str(num) for num in values) return "{{..., {0}, ...}}".format(", ".join(values))
gcd = fractions.gcd gcd.__doc__ += """ >>> gcd(12, 8) 4 >>> gcd(0, 50) 50 >>> gcd(7, 101) 1 """
[docs]def extended_gcd(a, b): """From `this exposition of the extended gcd algorithm <http://anh.cs.luc.edu/331/notes/xgcd.pdf>`_. Computes :math:`d = \gcd(a, b)` and returns a triple :math:`(d, x, y)` where :math:`d = ax + by`. >>> extended_gcd(14, 33) (1, -7, 3) >>> extended_gcd(12, 32) (4, 3, -1) """ prevx, x = 1, 0; prevy, y = 0, 1 while b: q = a//b x, prevx = prevx - q*x, x y, prevy = prevy - q*y, y a, b = b, a % b return a, prevx, prevy
[docs]def lcm(*args): """Computes the least common multiple of the given *args*. If a single argument is provided, the least common multiple of its elements is computed. >>> n = randint(1, 1000000) >>> lcm(n) == n True >>> lcm(2, 13) 26 >>> lcm(*range(1, 10)) #1, 2, ..., 9 2520 >>> lcm(range(1, 10)) #Don't have to unpack the arguments 2520 """ #See also http://stackoverflow.com/a/147539 if len(args) == 0: raise ValueError('No arguments provided.') elif len(args) == 1: try: iterator = iter(args[0]) except TypeError: return args[0] else: args = list(iterator) return reduce(_lcm2, args)
def _lcm2(a, b): return a*b // gcd(a, b)
[docs]def solve_linear_diophantine(a, b, c): r"""Solves the equation :math:`ax + by = c` for integers :math:`x, y \in\mathbb{Z}`. :rtype: ``None`` if no solution exists; otherwise a triple *(base, inc, lcm)* where - *base* is a single solution :math:`(x_0, y_0)`; - *inc* is a pair :math:`(\delta x, \delta y)` such that if :math:`(x, y)` is a solution, so is :math:`(x, y) \pm (\delta x, \delta y)`; and - *lcm* is the value of :math:`\operatorname{lcm}(a, b)`. """ d, x, y = extended_gcd(a, b) if c % d != 0: #Solution exists iff d divides c return None scale = c // d x *= scale y *= scale assert a*x + b*y == c base = (x, y) inc = (b // d, a // d) lcm = a * b // d return base, inc, lcm
[docs]def solve_linear_congruence(a, b, n): """Solves the congruence :math:`ax \equiv b \pmod n`. :rtype: a :class:`~thompson.orbits.SolutionSet` representing all such solutions :math:`x`. .. doctest:: :options: -ELLIPSIS >>> x = solve_linear_congruence(6, 12, 18) >>> print(x) {..., -1, 2, 5, 8, 11, 14, ...} >>> y = solve_linear_congruence(5, 7, 11) >>> print(y) {..., -3, 8, 19, 30, 41, 52, ...} >>> print(x & y) {..., -25, 8, 41, 74, 107, 140, ...} """ d = gcd(a, n) #1. A solution exists iff d divides b. if b % d != 0: return SolutionSet.empty_set() #2. Divide through by d to obtain an equivalent equation with smaller modulus. a //= d b //= d n //= d #3. At this stage, gcd(a, n) = 1. Appeal to Euclid to compute x = a^-1 mod n. d, x, y = extended_gcd(a, n) assert a*x + n*y == d == 1 #4. We have one concrete solution. Add multiples of n to get the rest. base = (x*b) % n return SolutionSet(base, n)
[docs]def divisors(n, include_one=True): """An iterator that yields the positive divisors :math:`d \mid n`. :param bool include_one: set to False to exlude :math:`d = 1`. .. doctest:: >>> list(divisors(12)) [1, 2, 3, 4, 6, 12] >>> list(divisors(125, include_one=False)) [5, 25, 125] """ if include_one: yield 1 for i in range(2, n): if n % i == 0: yield i yield n
[docs]def prod(iterable): """Handy function for computing the product of an iterable collection of numbers. From `Stack Overflow <http://stackoverflow.com/a/7948307>`_. >>> prod(range(1, 5)) 24 """ return reduce(mul, iterable, 1)