intfunc.py 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530
  1. """
  2. The routines here were removed from numbers.py, power.py,
  3. digits.py and factor_.py so they could be imported into core
  4. without raising circular import errors.
  5. Although the name 'intfunc' was chosen to represent functions that
  6. work with integers, it can also be thought of as containing
  7. internal/core functions that are needed by the classes of the core.
  8. """
  9. import math
  10. import sys
  11. from functools import lru_cache
  12. from .sympify import sympify
  13. from .singleton import S
  14. from sympy.external.gmpy import (gcd as number_gcd, lcm as number_lcm, sqrt,
  15. iroot, bit_scan1, gcdext)
  16. from sympy.utilities.misc import as_int, filldedent
  17. def num_digits(n, base=10):
  18. """Return the number of digits needed to express n in give base.
  19. Examples
  20. ========
  21. >>> from sympy.core.intfunc import num_digits
  22. >>> num_digits(10)
  23. 2
  24. >>> num_digits(10, 2) # 1010 -> 4 digits
  25. 4
  26. >>> num_digits(-100, 16) # -64 -> 2 digits
  27. 2
  28. Parameters
  29. ==========
  30. n: integer
  31. The number whose digits are counted.
  32. b: integer
  33. The base in which digits are computed.
  34. See Also
  35. ========
  36. sympy.ntheory.digits.digits, sympy.ntheory.digits.count_digits
  37. """
  38. if base < 0:
  39. raise ValueError('base must be int greater than 1')
  40. if not n:
  41. return 1
  42. e, t = integer_log(abs(n), base)
  43. return 1 + e
  44. def integer_log(n, b):
  45. r"""
  46. Returns ``(e, bool)`` where e is the largest nonnegative integer
  47. such that :math:`|n| \geq |b^e|` and ``bool`` is True if $n = b^e$.
  48. Examples
  49. ========
  50. >>> from sympy import integer_log
  51. >>> integer_log(125, 5)
  52. (3, True)
  53. >>> integer_log(17, 9)
  54. (1, False)
  55. If the base is positive and the number negative the
  56. return value will always be the same except for 2:
  57. >>> integer_log(-4, 2)
  58. (2, False)
  59. >>> integer_log(-16, 4)
  60. (0, False)
  61. When the base is negative, the returned value
  62. will only be True if the parity of the exponent is
  63. correct for the sign of the base:
  64. >>> integer_log(4, -2)
  65. (2, True)
  66. >>> integer_log(8, -2)
  67. (3, False)
  68. >>> integer_log(-8, -2)
  69. (3, True)
  70. >>> integer_log(-4, -2)
  71. (2, False)
  72. See Also
  73. ========
  74. integer_nthroot
  75. sympy.ntheory.primetest.is_square
  76. sympy.ntheory.factor_.multiplicity
  77. sympy.ntheory.factor_.perfect_power
  78. """
  79. n = as_int(n)
  80. b = as_int(b)
  81. if b < 0:
  82. e, t = integer_log(abs(n), -b)
  83. # (-2)**3 == -8
  84. # (-2)**2 = 4
  85. t = t and e % 2 == (n < 0)
  86. return e, t
  87. if b <= 1:
  88. raise ValueError('base must be 2 or more')
  89. if n < 0:
  90. if b != 2:
  91. return 0, False
  92. e, t = integer_log(-n, b)
  93. return e, False
  94. if n == 0:
  95. raise ValueError('n cannot be 0')
  96. if n < b:
  97. return 0, n == 1
  98. if b == 2:
  99. e = n.bit_length() - 1
  100. return e, trailing(n) == e
  101. t = trailing(b)
  102. if 2**t == b:
  103. e = int(n.bit_length() - 1)//t
  104. n_ = 1 << (t*e)
  105. return e, n_ == n
  106. d = math.floor(math.log10(n) / math.log10(b))
  107. n_ = b ** d
  108. while n_ <= n: # this will iterate 0, 1 or 2 times
  109. d += 1
  110. n_ *= b
  111. return d - (n_ > n), (n_ == n or n_//b == n)
  112. def trailing(n):
  113. """Count the number of trailing zero digits in the binary
  114. representation of n, i.e. determine the largest power of 2
  115. that divides n.
  116. Examples
  117. ========
  118. >>> from sympy import trailing
  119. >>> trailing(128)
  120. 7
  121. >>> trailing(63)
  122. 0
  123. See Also
  124. ========
  125. sympy.ntheory.factor_.multiplicity
  126. """
  127. if not n:
  128. return 0
  129. return bit_scan1(int(n))
  130. @lru_cache(1024)
  131. def igcd(*args):
  132. """Computes nonnegative integer greatest common divisor.
  133. Explanation
  134. ===========
  135. The algorithm is based on the well known Euclid's algorithm [1]_. To
  136. improve speed, ``igcd()`` has its own caching mechanism.
  137. If you do not need the cache mechanism, using ``sympy.external.gmpy.gcd``.
  138. Examples
  139. ========
  140. >>> from sympy import igcd
  141. >>> igcd(2, 4)
  142. 2
  143. >>> igcd(5, 10, 15)
  144. 5
  145. References
  146. ==========
  147. .. [1] https://en.wikipedia.org/wiki/Euclidean_algorithm
  148. """
  149. if len(args) < 2:
  150. raise TypeError("igcd() takes at least 2 arguments (%s given)" % len(args))
  151. return int(number_gcd(*map(as_int, args)))
  152. igcd2 = math.gcd
  153. def igcd_lehmer(a, b):
  154. r"""Computes greatest common divisor of two integers.
  155. Explanation
  156. ===========
  157. Euclid's algorithm for the computation of the greatest
  158. common divisor ``gcd(a, b)`` of two (positive) integers
  159. $a$ and $b$ is based on the division identity
  160. $$ a = q \times b + r$$,
  161. where the quotient $q$ and the remainder $r$ are integers
  162. and $0 \le r < b$. Then each common divisor of $a$ and $b$
  163. divides $r$, and it follows that ``gcd(a, b) == gcd(b, r)``.
  164. The algorithm works by constructing the sequence
  165. r0, r1, r2, ..., where r0 = a, r1 = b, and each rn
  166. is the remainder from the division of the two preceding
  167. elements.
  168. In Python, ``q = a // b`` and ``r = a % b`` are obtained by the
  169. floor division and the remainder operations, respectively.
  170. These are the most expensive arithmetic operations, especially
  171. for large a and b.
  172. Lehmer's algorithm [1]_ is based on the observation that the quotients
  173. ``qn = r(n-1) // rn`` are in general small integers even
  174. when a and b are very large. Hence the quotients can be
  175. usually determined from a relatively small number of most
  176. significant bits.
  177. The efficiency of the algorithm is further enhanced by not
  178. computing each long remainder in Euclid's sequence. The remainders
  179. are linear combinations of a and b with integer coefficients
  180. derived from the quotients. The coefficients can be computed
  181. as far as the quotients can be determined from the chosen
  182. most significant parts of a and b. Only then a new pair of
  183. consecutive remainders is computed and the algorithm starts
  184. anew with this pair.
  185. References
  186. ==========
  187. .. [1] https://en.wikipedia.org/wiki/Lehmer%27s_GCD_algorithm
  188. """
  189. a, b = abs(as_int(a)), abs(as_int(b))
  190. if a < b:
  191. a, b = b, a
  192. # The algorithm works by using one or two digit division
  193. # whenever possible. The outer loop will replace the
  194. # pair (a, b) with a pair of shorter consecutive elements
  195. # of the Euclidean gcd sequence until a and b
  196. # fit into two Python (long) int digits.
  197. nbits = 2 * sys.int_info.bits_per_digit
  198. while a.bit_length() > nbits and b != 0:
  199. # Quotients are mostly small integers that can
  200. # be determined from most significant bits.
  201. n = a.bit_length() - nbits
  202. x, y = int(a >> n), int(b >> n) # most significant bits
  203. # Elements of the Euclidean gcd sequence are linear
  204. # combinations of a and b with integer coefficients.
  205. # Compute the coefficients of consecutive pairs
  206. # a' = A*a + B*b, b' = C*a + D*b
  207. # using small integer arithmetic as far as possible.
  208. A, B, C, D = 1, 0, 0, 1 # initial values
  209. while True:
  210. # The coefficients alternate in sign while looping.
  211. # The inner loop combines two steps to keep track
  212. # of the signs.
  213. # At this point we have
  214. # A > 0, B <= 0, C <= 0, D > 0,
  215. # x' = x + B <= x < x" = x + A,
  216. # y' = y + C <= y < y" = y + D,
  217. # and
  218. # x'*N <= a' < x"*N, y'*N <= b' < y"*N,
  219. # where N = 2**n.
  220. # Now, if y' > 0, and x"//y' and x'//y" agree,
  221. # then their common value is equal to q = a'//b'.
  222. # In addition,
  223. # x'%y" = x' - q*y" < x" - q*y' = x"%y',
  224. # and
  225. # (x'%y")*N < a'%b' < (x"%y')*N.
  226. # On the other hand, we also have x//y == q,
  227. # and therefore
  228. # x'%y" = x + B - q*(y + D) = x%y + B',
  229. # x"%y' = x + A - q*(y + C) = x%y + A',
  230. # where
  231. # B' = B - q*D < 0, A' = A - q*C > 0.
  232. if y + C <= 0:
  233. break
  234. q = (x + A) // (y + C)
  235. # Now x'//y" <= q, and equality holds if
  236. # x' - q*y" = (x - q*y) + (B - q*D) >= 0.
  237. # This is a minor optimization to avoid division.
  238. x_qy, B_qD = x - q * y, B - q * D
  239. if x_qy + B_qD < 0:
  240. break
  241. # Next step in the Euclidean sequence.
  242. x, y = y, x_qy
  243. A, B, C, D = C, D, A - q * C, B_qD
  244. # At this point the signs of the coefficients
  245. # change and their roles are interchanged.
  246. # A <= 0, B > 0, C > 0, D < 0,
  247. # x' = x + A <= x < x" = x + B,
  248. # y' = y + D < y < y" = y + C.
  249. if y + D <= 0:
  250. break
  251. q = (x + B) // (y + D)
  252. x_qy, A_qC = x - q * y, A - q * C
  253. if x_qy + A_qC < 0:
  254. break
  255. x, y = y, x_qy
  256. A, B, C, D = C, D, A_qC, B - q * D
  257. # Now the conditions on top of the loop
  258. # are again satisfied.
  259. # A > 0, B < 0, C < 0, D > 0.
  260. if B == 0:
  261. # This can only happen when y == 0 in the beginning
  262. # and the inner loop does nothing.
  263. # Long division is forced.
  264. a, b = b, a % b
  265. continue
  266. # Compute new long arguments using the coefficients.
  267. a, b = A * a + B * b, C * a + D * b
  268. # Small divisors. Finish with the standard algorithm.
  269. while b:
  270. a, b = b, a % b
  271. return a
  272. def ilcm(*args):
  273. """Computes integer least common multiple.
  274. Examples
  275. ========
  276. >>> from sympy import ilcm
  277. >>> ilcm(5, 10)
  278. 10
  279. >>> ilcm(7, 3)
  280. 21
  281. >>> ilcm(5, 10, 15)
  282. 30
  283. """
  284. if len(args) < 2:
  285. raise TypeError("ilcm() takes at least 2 arguments (%s given)" % len(args))
  286. return int(number_lcm(*map(as_int, args)))
  287. def igcdex(a, b):
  288. """Returns x, y, g such that g = x*a + y*b = gcd(a, b).
  289. Examples
  290. ========
  291. >>> from sympy.core.intfunc import igcdex
  292. >>> igcdex(2, 3)
  293. (-1, 1, 1)
  294. >>> igcdex(10, 12)
  295. (-1, 1, 2)
  296. >>> x, y, g = igcdex(100, 2004)
  297. >>> x, y, g
  298. (-20, 1, 4)
  299. >>> x*100 + y*2004
  300. 4
  301. """
  302. g, x, y = gcdext(int(a), int(b))
  303. return x, y, g
  304. def mod_inverse(a, m):
  305. r"""
  306. Return the number $c$ such that, $a \times c = 1 \pmod{m}$
  307. where $c$ has the same sign as $m$. If no such value exists,
  308. a ValueError is raised.
  309. Examples
  310. ========
  311. >>> from sympy import mod_inverse, S
  312. Suppose we wish to find multiplicative inverse $x$ of
  313. 3 modulo 11. This is the same as finding $x$ such
  314. that $3x = 1 \pmod{11}$. One value of x that satisfies
  315. this congruence is 4. Because $3 \times 4 = 12$ and $12 = 1 \pmod{11}$.
  316. This is the value returned by ``mod_inverse``:
  317. >>> mod_inverse(3, 11)
  318. 4
  319. >>> mod_inverse(-3, 11)
  320. 7
  321. When there is a common factor between the numerators of
  322. `a` and `m` the inverse does not exist:
  323. >>> mod_inverse(2, 4)
  324. Traceback (most recent call last):
  325. ...
  326. ValueError: inverse of 2 mod 4 does not exist
  327. >>> mod_inverse(S(2)/7, S(5)/2)
  328. 7/2
  329. References
  330. ==========
  331. .. [1] https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
  332. .. [2] https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
  333. """
  334. c = None
  335. try:
  336. a, m = as_int(a), as_int(m)
  337. if m != 1 and m != -1:
  338. x, _, g = igcdex(a, m)
  339. if g == 1:
  340. c = x % m
  341. except ValueError:
  342. a, m = sympify(a), sympify(m)
  343. if not (a.is_number and m.is_number):
  344. raise TypeError(
  345. filldedent(
  346. """
  347. Expected numbers for arguments; symbolic `mod_inverse`
  348. is not implemented
  349. but symbolic expressions can be handled with the
  350. similar function,
  351. sympy.polys.polytools.invert"""
  352. )
  353. )
  354. big = m > 1
  355. if big not in (S.true, S.false):
  356. raise ValueError("m > 1 did not evaluate; try to simplify %s" % m)
  357. elif big:
  358. c = 1 / a
  359. if c is None:
  360. raise ValueError("inverse of %s (mod %s) does not exist" % (a, m))
  361. return c
  362. def isqrt(n):
  363. r""" Return the largest integer less than or equal to `\sqrt{n}`.
  364. Parameters
  365. ==========
  366. n : non-negative integer
  367. Returns
  368. =======
  369. int : `\left\lfloor\sqrt{n}\right\rfloor`
  370. Raises
  371. ======
  372. ValueError
  373. If n is negative.
  374. TypeError
  375. If n is of a type that cannot be compared to ``int``.
  376. Therefore, a TypeError is raised for ``str``, but not for ``float``.
  377. Examples
  378. ========
  379. >>> from sympy.core.intfunc import isqrt
  380. >>> isqrt(0)
  381. 0
  382. >>> isqrt(9)
  383. 3
  384. >>> isqrt(10)
  385. 3
  386. >>> isqrt("30")
  387. Traceback (most recent call last):
  388. ...
  389. TypeError: '<' not supported between instances of 'str' and 'int'
  390. >>> from sympy.core.numbers import Rational
  391. >>> isqrt(Rational(-1, 2))
  392. Traceback (most recent call last):
  393. ...
  394. ValueError: n must be nonnegative
  395. """
  396. if n < 0:
  397. raise ValueError("n must be nonnegative")
  398. return int(sqrt(int(n)))
  399. def integer_nthroot(y, n):
  400. """
  401. Return a tuple containing x = floor(y**(1/n))
  402. and a boolean indicating whether the result is exact (that is,
  403. whether x**n == y).
  404. Examples
  405. ========
  406. >>> from sympy import integer_nthroot
  407. >>> integer_nthroot(16, 2)
  408. (4, True)
  409. >>> integer_nthroot(26, 2)
  410. (5, False)
  411. To simply determine if a number is a perfect square, the is_square
  412. function should be used:
  413. >>> from sympy.ntheory.primetest import is_square
  414. >>> is_square(26)
  415. False
  416. See Also
  417. ========
  418. sympy.ntheory.primetest.is_square
  419. integer_log
  420. """
  421. x, b = iroot(as_int(y), as_int(n))
  422. return int(x), b