Source Code for Module sympycore.arithmetic.number_theory

  1  """Provides algorithms from number theory.
 
  2  """ 
  3  from .numbers import FractionTuple, normalized_fraction, Complex, Float, div 
  4  
 
  5  __all__ = ['gcd', 'lcm', 'factorial',
 
  6             'integer_digits', 'real_digits',
 
  7             'multinomial_coefficients'] 
  8  
 
  9  __docformat__ = "restructuredtext en" 
 10  
 
 11 -def factorial(n, memo=[1, 1]): 
 12      """Return n factorial (for integers n >= 0 only).""" 
 13      if n < 0: 
 14          raise ValueError 
 15      k = len(memo) 
 16      if n < k: 
 17          return memo[n] 
 18      p = memo[-1] 
 19      while k <= n: 
 20          p *= k 
 21          k += 1 
 22          if k < 100: 
 23              memo.append(p) 
 24      return p 
 25  
 
 26 -def gcd(*args): 
 27      """Calculate the greatest common divisor (GCD) of the arguments.""" 
 28      L = len(args) 
 29      if L == 0: return 0 
 30      if L == 1: return args[0] 
 31      if L == 2: 
 32          a, b = args 
 33          while b: 
 34              a, b = b, a % b 
 35          return a 
 36      return gcd(gcd(args[0], args[1]), *args[2:]) 
 37  
 
 38 -def lcm(*args): 
 39      """Calculate the least common multiple (LCM) of the arguments.""" 
 40      L = len(args) 
 41      if L == 0: return 0 
 42      if L == 1: return args[0] 
 43      if L == 2: return div(args[0], gcd(*args))*args[1] 
 44      return lcm(lcm(args[0], args[1]), *args[2:]) 
 45  
 
 46  # TODO: this could use the faster implementation in mpmath
 
 47 -def integer_digits(n, base=10): 
 48      """Return a list of the digits of abs(n) in the given base.""" 
 49      assert base > 1 
 50      assert isinstance(n, (int, long)) 
 51      n = abs(n) 
 52      if not n: 
 53          return [0] 
 54      L = [] 
 55      while n: 
 56          n, digit = divmod(n, base) 
 57          L.append(int(digit)) 
 58      return L[::-1] 
 59  
 
 60  # TODO: this could (also?) be implemented as an endless generator
 
 61 -def real_digits(x, base=10, truncation=10): 
 62      """Return ``(L, d)`` where L is a list of digits of ``abs(x)`` in
 
 63      the given base and ``d`` is the (signed) distance from the leading
 
 64      digit to the radix point.
 
 65  
 
 66      For example, 1234.56 becomes ``([1, 2, 3, 4, 5, 6], 4)`` and 0.001
 
 67      becomes ``([1], -2)``. If, during the generation of fractional
 
 68      digits, the length reaches `truncation` digits, the iteration is
 
 69      stopped.""" 
 70      assert base > 1 
 71      assert isinstance(x, (int, long, FractionTuple)) 
 72      if x == 0: 
 73          return ([0], 1) 
 74      x = abs(x) 
 75      exponent = 0 
 76      while x < 1: 
 77          x *= base 
 78          exponent -= 1 
 79      integer, fraction = divmod(x, 1) 
 80      L = integer_digits(integer, base) 
 81      exponent += len(L) 
 82      if fraction: 
 83          p, q = fraction 
 84          for i in xrange(truncation - len(L)): 
 85              p = (p % q) * base 
 86              if not p: 
 87                  break 
 88              L.append(int(p//q)) 
 89      return L, exponent 
 90  
 
 91 -def binomial_coefficients(n): 
 92      """Return a dictionary containing pairs {(k1,k2) : C_kn} where
 
 93      C_kn are binomial coefficients and n=k1+k2.""" 
 94      d = {(0, n):1, (n, 0):1} 
 95      a = 1 
 96      for k in xrange(1, n//2+1): 
 97          a = (a * (n-k+1))//k 
 98          d[k, n-k] = d[n-k, k] = a 
 99      return d 
100  
 
101 -def binomial_coefficients_list(n): 
102      d = [1] * (n+1) 
103      a = 1 
104      for k in xrange(1, n//2+1): 
105          a = (a * (n-k+1))//k 
106          d[k] = d[n-k] = a 
107      return d 
108  
 
109 -def multinomial_coefficients(m, n, _tuple=tuple, _zip=zip): 
110      """Return a dictionary containing pairs ``{(k1,k2,..,km) : C_kn}``
 
111      where ``C_kn`` are multinomial coefficients such that
 
112      ``n=k1+k2+..+km``.
 
113  
 
114      For example:
 
115  
 
116      >>> print multinomial_coefficients(2,5)
 
117      {(3, 2): 10, (1, 4): 5, (2, 3): 10, (5, 0): 1, (0, 5): 1, (4, 1): 5}
 
118  
 
119      The algorithm is based on the following result:
 
120      
 
121         Consider a polynomial and it's ``m``-th exponent::
 
122         
 
123           P(x) = sum_{i=0}^m p_i x^k
 
124           P(x)^n = sum_{k=0}^{m n} a(n,k) x^k
 
125  
 
126         The coefficients ``a(n,k)`` can be computed using the
 
127         J.C.P. Miller Pure Recurrence [see D.E.Knuth, Seminumerical
 
128         Algorithms, The art of Computer Programming v.2, Addison
 
129         Wesley, Reading, 1981;]::
 
130         
 
131           a(n,k) = 1/(k p_0) sum_{i=1}^m p_i ((n+1)i-k) a(n,k-i),
 
132  
 
133         where ``a(n,0) = p_0^n``.
 
134      """ 
135  
 
136      if m==2: 
137          return binomial_coefficients(n) 
138      symbols = [(0,)*i + (1,) + (0,)*(m-i-1) for i in range(m)] 
139      s0 = symbols[0] 
140      p0 = [_tuple(aa-bb for aa,bb in _zip(s,s0)) for s in symbols] 
141      r = {_tuple(aa*n for aa in s0):1} 
142      r_get = r.get 
143      r_update = r.update 
144      l = [0] * (n*(m-1)+1) 
145      l[0] = r.items() 
146      for k in xrange(1, n*(m-1)+1): 
147          d = {} 
148          d_get = d.get 
149          for i in xrange(1, min(m,k+1)): 
150              nn = (n+1)*i-k 
151              if not nn: 
152                  continue 
153              t = p0[i] 
154              for t2, c2 in l[k-i]: 
155                  tt = _tuple([aa+bb for aa,bb in _zip(t2,t)]) 
156                  cc = nn * c2 
157                  b = d_get(tt) 
158                  if b is None: 
159                      d[tt] = cc 
160                  else: 
161                      cc = b + cc 
162                      if cc: 
163                          d[tt] = cc 
164                      else: 
165                          del d[tt] 
166          r1 = [(t, c//k) for (t, c) in d.iteritems()] 
167          l[k] = r1 
168          r_update(r1) 
169      return r 
170