Coverage for /usr/lib/python3/dist-packages/sympy/ntheory/factor_.py: 23%
907 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2Integer factorization
3"""
5from collections import defaultdict
6from functools import reduce
7import random
8import math
10from sympy.core import sympify
11from sympy.core.containers import Dict
12from sympy.core.evalf import bitcount
13from sympy.core.expr import Expr
14from sympy.core.function import Function
15from sympy.core.logic import fuzzy_and
16from sympy.core.mul import Mul
17from sympy.core.numbers import igcd, ilcm, Rational, Integer
18from sympy.core.power import integer_nthroot, Pow, integer_log
19from sympy.core.singleton import S
20from sympy.external.gmpy import SYMPY_INTS
21from .primetest import isprime
22from .generate import sieve, primerange, nextprime
23from .digits import digits
24from sympy.utilities.iterables import flatten
25from sympy.utilities.misc import as_int, filldedent
26from .ecm import _ecm_one_factor
28# Note: This list should be updated whenever new Mersenne primes are found.
29# Refer: https://www.mersenne.org/
30MERSENNE_PRIME_EXPONENTS = (2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203,
31 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209, 44497, 86243, 110503, 132049,
32 216091, 756839, 859433, 1257787, 1398269, 2976221, 3021377, 6972593, 13466917, 20996011, 24036583,
33 25964951, 30402457, 32582657, 37156667, 42643801, 43112609, 57885161, 74207281, 77232917, 82589933)
35# compute more when needed for i in Mersenne prime exponents
36PERFECT = [6] # 2**(i-1)*(2**i-1)
37MERSENNES = [3] # 2**i - 1
40def _ismersenneprime(n):
41 global MERSENNES
42 j = len(MERSENNES)
43 while n > MERSENNES[-1] and j < len(MERSENNE_PRIME_EXPONENTS):
44 # conservatively grow the list
45 MERSENNES.append(2**MERSENNE_PRIME_EXPONENTS[j] - 1)
46 j += 1
47 return n in MERSENNES
50def _isperfect(n):
51 global PERFECT
52 if n % 2 == 0:
53 j = len(PERFECT)
54 while n > PERFECT[-1] and j < len(MERSENNE_PRIME_EXPONENTS):
55 # conservatively grow the list
56 t = 2**(MERSENNE_PRIME_EXPONENTS[j] - 1)
57 PERFECT.append(t*(2*t - 1))
58 j += 1
59 return n in PERFECT
62small_trailing = [0] * 256
63for j in range(1,8):
64 small_trailing[1<<j::1<<(j+1)] = [j] * (1<<(7-j))
67def smoothness(n):
68 """
69 Return the B-smooth and B-power smooth values of n.
71 The smoothness of n is the largest prime factor of n; the power-
72 smoothness is the largest divisor raised to its multiplicity.
74 Examples
75 ========
77 >>> from sympy.ntheory.factor_ import smoothness
78 >>> smoothness(2**7*3**2)
79 (3, 128)
80 >>> smoothness(2**4*13)
81 (13, 16)
82 >>> smoothness(2)
83 (2, 2)
85 See Also
86 ========
88 factorint, smoothness_p
89 """
91 if n == 1:
92 return (1, 1) # not prime, but otherwise this causes headaches
93 facs = factorint(n)
94 return max(facs), max(m**facs[m] for m in facs)
97def smoothness_p(n, m=-1, power=0, visual=None):
98 """
99 Return a list of [m, (p, (M, sm(p + m), psm(p + m)))...]
100 where:
102 1. p**M is the base-p divisor of n
103 2. sm(p + m) is the smoothness of p + m (m = -1 by default)
104 3. psm(p + m) is the power smoothness of p + m
106 The list is sorted according to smoothness (default) or by power smoothness
107 if power=1.
109 The smoothness of the numbers to the left (m = -1) or right (m = 1) of a
110 factor govern the results that are obtained from the p +/- 1 type factoring
111 methods.
113 >>> from sympy.ntheory.factor_ import smoothness_p, factorint
114 >>> smoothness_p(10431, m=1)
115 (1, [(3, (2, 2, 4)), (19, (1, 5, 5)), (61, (1, 31, 31))])
116 >>> smoothness_p(10431)
117 (-1, [(3, (2, 2, 2)), (19, (1, 3, 9)), (61, (1, 5, 5))])
118 >>> smoothness_p(10431, power=1)
119 (-1, [(3, (2, 2, 2)), (61, (1, 5, 5)), (19, (1, 3, 9))])
121 If visual=True then an annotated string will be returned:
123 >>> print(smoothness_p(21477639576571, visual=1))
124 p**i=4410317**1 has p-1 B=1787, B-pow=1787
125 p**i=4869863**1 has p-1 B=2434931, B-pow=2434931
127 This string can also be generated directly from a factorization dictionary
128 and vice versa:
130 >>> factorint(17*9)
131 {3: 2, 17: 1}
132 >>> smoothness_p(_)
133 'p**i=3**2 has p-1 B=2, B-pow=2\\np**i=17**1 has p-1 B=2, B-pow=16'
134 >>> smoothness_p(_)
135 {3: 2, 17: 1}
137 The table of the output logic is:
139 ====== ====== ======= =======
140 | Visual
141 ------ ----------------------
142 Input True False other
143 ====== ====== ======= =======
144 dict str tuple str
145 str str tuple dict
146 tuple str tuple str
147 n str tuple tuple
148 mul str tuple tuple
149 ====== ====== ======= =======
151 See Also
152 ========
154 factorint, smoothness
155 """
157 # visual must be True, False or other (stored as None)
158 if visual in (1, 0):
159 visual = bool(visual)
160 elif visual not in (True, False):
161 visual = None
163 if isinstance(n, str):
164 if visual:
165 return n
166 d = {}
167 for li in n.splitlines():
168 k, v = [int(i) for i in
169 li.split('has')[0].split('=')[1].split('**')]
170 d[k] = v
171 if visual is not True and visual is not False:
172 return d
173 return smoothness_p(d, visual=False)
174 elif not isinstance(n, tuple):
175 facs = factorint(n, visual=False)
177 if power:
178 k = -1
179 else:
180 k = 1
181 if isinstance(n, tuple):
182 rv = n
183 else:
184 rv = (m, sorted([(f,
185 tuple([M] + list(smoothness(f + m))))
186 for f, M in list(facs.items())],
187 key=lambda x: (x[1][k], x[0])))
189 if visual is False or (visual is not True) and (type(n) in [int, Mul]):
190 return rv
191 lines = []
192 for dat in rv[1]:
193 dat = flatten(dat)
194 dat.insert(2, m)
195 lines.append('p**i=%i**%i has p%+i B=%i, B-pow=%i' % tuple(dat))
196 return '\n'.join(lines)
199def trailing(n):
200 """Count the number of trailing zero digits in the binary
201 representation of n, i.e. determine the largest power of 2
202 that divides n.
204 Examples
205 ========
207 >>> from sympy import trailing
208 >>> trailing(128)
209 7
210 >>> trailing(63)
211 0
212 """
213 n = abs(int(n))
214 if not n:
215 return 0
216 low_byte = n & 0xff
217 if low_byte:
218 return small_trailing[low_byte]
220 # 2**m is quick for z up through 2**30
221 z = bitcount(n) - 1
222 if isinstance(z, SYMPY_INTS):
223 if n == 1 << z:
224 return z
226 if z < 300:
227 # fixed 8-byte reduction
228 t = 8
229 n >>= 8
230 while not n & 0xff:
231 n >>= 8
232 t += 8
233 return t + small_trailing[n & 0xff]
235 # binary reduction important when there might be a large
236 # number of trailing 0s
237 t = 0
238 p = 8
239 while not n & 1:
240 while not n & ((1 << p) - 1):
241 n >>= p
242 t += p
243 p *= 2
244 p //= 2
245 return t
248def multiplicity(p, n):
249 """
250 Find the greatest integer m such that p**m divides n.
252 Examples
253 ========
255 >>> from sympy import multiplicity, Rational
256 >>> [multiplicity(5, n) for n in [8, 5, 25, 125, 250]]
257 [0, 1, 2, 3, 3]
258 >>> multiplicity(3, Rational(1, 9))
259 -2
261 Note: when checking for the multiplicity of a number in a
262 large factorial it is most efficient to send it as an unevaluated
263 factorial or to call ``multiplicity_in_factorial`` directly:
265 >>> from sympy.ntheory import multiplicity_in_factorial
266 >>> from sympy import factorial
267 >>> p = factorial(25)
268 >>> n = 2**100
269 >>> nfac = factorial(n, evaluate=False)
270 >>> multiplicity(p, nfac)
271 52818775009509558395695966887
272 >>> _ == multiplicity_in_factorial(p, n)
273 True
275 """
276 try:
277 p, n = as_int(p), as_int(n)
278 except ValueError:
279 from sympy.functions.combinatorial.factorials import factorial
280 if all(isinstance(i, (SYMPY_INTS, Rational)) for i in (p, n)):
281 p = Rational(p)
282 n = Rational(n)
283 if p.q == 1:
284 if n.p == 1:
285 return -multiplicity(p.p, n.q)
286 return multiplicity(p.p, n.p) - multiplicity(p.p, n.q)
287 elif p.p == 1:
288 return multiplicity(p.q, n.q)
289 else:
290 like = min(
291 multiplicity(p.p, n.p),
292 multiplicity(p.q, n.q))
293 cross = min(
294 multiplicity(p.q, n.p),
295 multiplicity(p.p, n.q))
296 return like - cross
297 elif (isinstance(p, (SYMPY_INTS, Integer)) and
298 isinstance(n, factorial) and
299 isinstance(n.args[0], Integer) and
300 n.args[0] >= 0):
301 return multiplicity_in_factorial(p, n.args[0])
302 raise ValueError('expecting ints or fractions, got %s and %s' % (p, n))
304 if n == 0:
305 raise ValueError('no such integer exists: multiplicity of %s is not-defined' %(n))
306 if p == 2:
307 return trailing(n)
308 if p < 2:
309 raise ValueError('p must be an integer, 2 or larger, but got %s' % p)
310 if p == n:
311 return 1
313 m = 0
314 n, rem = divmod(n, p)
315 while not rem:
316 m += 1
317 if m > 5:
318 # The multiplicity could be very large. Better
319 # to increment in powers of two
320 e = 2
321 while 1:
322 ppow = p**e
323 if ppow < n:
324 nnew, rem = divmod(n, ppow)
325 if not rem:
326 m += e
327 e *= 2
328 n = nnew
329 continue
330 return m + multiplicity(p, n)
331 n, rem = divmod(n, p)
332 return m
335def multiplicity_in_factorial(p, n):
336 """return the largest integer ``m`` such that ``p**m`` divides ``n!``
337 without calculating the factorial of ``n``.
340 Examples
341 ========
343 >>> from sympy.ntheory import multiplicity_in_factorial
344 >>> from sympy import factorial
346 >>> multiplicity_in_factorial(2, 3)
347 1
349 An instructive use of this is to tell how many trailing zeros
350 a given factorial has. For example, there are 6 in 25!:
352 >>> factorial(25)
353 15511210043330985984000000
354 >>> multiplicity_in_factorial(10, 25)
355 6
357 For large factorials, it is much faster/feasible to use
358 this function rather than computing the actual factorial:
360 >>> multiplicity_in_factorial(factorial(25), 2**100)
361 52818775009509558395695966887
363 """
365 p, n = as_int(p), as_int(n)
367 if p <= 0:
368 raise ValueError('expecting positive integer got %s' % p )
370 if n < 0:
371 raise ValueError('expecting non-negative integer got %s' % n )
373 factors = factorint(p)
375 # keep only the largest of a given multiplicity since those
376 # of a given multiplicity will be goverened by the behavior
377 # of the largest factor
378 test = defaultdict(int)
379 for k, v in factors.items():
380 test[v] = max(k, test[v])
381 keep = set(test.values())
382 # remove others from factors
383 for k in list(factors.keys()):
384 if k not in keep:
385 factors.pop(k)
387 mp = S.Infinity
388 for i in factors:
389 # multiplicity of i in n! is
390 mi = (n - (sum(digits(n, i)) - i))//(i - 1)
391 # multiplicity of p in n! depends on multiplicity
392 # of prime `i` in p, so we floor divide by factors[i]
393 # and keep it if smaller than the multiplicity of p
394 # seen so far
395 mp = min(mp, mi//factors[i])
397 return mp
400def perfect_power(n, candidates=None, big=True, factor=True):
401 """
402 Return ``(b, e)`` such that ``n`` == ``b**e`` if ``n`` is a unique
403 perfect power with ``e > 1``, else ``False`` (e.g. 1 is not a
404 perfect power). A ValueError is raised if ``n`` is not Rational.
406 By default, the base is recursively decomposed and the exponents
407 collected so the largest possible ``e`` is sought. If ``big=False``
408 then the smallest possible ``e`` (thus prime) will be chosen.
410 If ``factor=True`` then simultaneous factorization of ``n`` is
411 attempted since finding a factor indicates the only possible root
412 for ``n``. This is True by default since only a few small factors will
413 be tested in the course of searching for the perfect power.
415 The use of ``candidates`` is primarily for internal use; if provided,
416 False will be returned if ``n`` cannot be written as a power with one
417 of the candidates as an exponent and factoring (beyond testing for
418 a factor of 2) will not be attempted.
420 Examples
421 ========
423 >>> from sympy import perfect_power, Rational
424 >>> perfect_power(16)
425 (2, 4)
426 >>> perfect_power(16, big=False)
427 (4, 2)
429 Negative numbers can only have odd perfect powers:
431 >>> perfect_power(-4)
432 False
433 >>> perfect_power(-8)
434 (-2, 3)
436 Rationals are also recognized:
438 >>> perfect_power(Rational(1, 2)**3)
439 (1/2, 3)
440 >>> perfect_power(Rational(-3, 2)**3)
441 (-3/2, 3)
443 Notes
444 =====
446 To know whether an integer is a perfect power of 2 use
448 >>> is2pow = lambda n: bool(n and not n & (n - 1))
449 >>> [(i, is2pow(i)) for i in range(5)]
450 [(0, False), (1, True), (2, True), (3, False), (4, True)]
452 It is not necessary to provide ``candidates``. When provided
453 it will be assumed that they are ints. The first one that is
454 larger than the computed maximum possible exponent will signal
455 failure for the routine.
457 >>> perfect_power(3**8, [9])
458 False
459 >>> perfect_power(3**8, [2, 4, 8])
460 (3, 8)
461 >>> perfect_power(3**8, [4, 8], big=False)
462 (9, 4)
464 See Also
465 ========
466 sympy.core.power.integer_nthroot
467 sympy.ntheory.primetest.is_square
468 """
469 if isinstance(n, Rational) and not n.is_Integer:
470 p, q = n.as_numer_denom()
471 if p is S.One:
472 pp = perfect_power(q)
473 if pp:
474 pp = (n.func(1, pp[0]), pp[1])
475 else:
476 pp = perfect_power(p)
477 if pp:
478 num, e = pp
479 pq = perfect_power(q, [e])
480 if pq:
481 den, _ = pq
482 pp = n.func(num, den), e
483 return pp
485 n = as_int(n)
486 if n < 0:
487 pp = perfect_power(-n)
488 if pp:
489 b, e = pp
490 if e % 2:
491 return -b, e
492 return False
494 if n <= 3:
495 # no unique exponent for 0, 1
496 # 2 and 3 have exponents of 1
497 return False
498 logn = math.log(n, 2)
499 max_possible = int(logn) + 2 # only check values less than this
500 not_square = n % 10 in [2, 3, 7, 8] # squares cannot end in 2, 3, 7, 8
501 min_possible = 2 + not_square
502 if not candidates:
503 candidates = primerange(min_possible, max_possible)
504 else:
505 candidates = sorted([i for i in candidates
506 if min_possible <= i < max_possible])
507 if n%2 == 0:
508 e = trailing(n)
509 candidates = [i for i in candidates if e%i == 0]
510 if big:
511 candidates = reversed(candidates)
512 for e in candidates:
513 r, ok = integer_nthroot(n, e)
514 if ok:
515 return (r, e)
516 return False
518 def _factors():
519 rv = 2 + n % 2
520 while True:
521 yield rv
522 rv = nextprime(rv)
524 for fac, e in zip(_factors(), candidates):
525 # see if there is a factor present
526 if factor and n % fac == 0:
527 # find what the potential power is
528 if fac == 2:
529 e = trailing(n)
530 else:
531 e = multiplicity(fac, n)
532 # if it's a trivial power we are done
533 if e == 1:
534 return False
536 # maybe the e-th root of n is exact
537 r, exact = integer_nthroot(n, e)
538 if not exact:
539 # Having a factor, we know that e is the maximal
540 # possible value for a root of n.
541 # If n = fac**e*m can be written as a perfect
542 # power then see if m can be written as r**E where
543 # gcd(e, E) != 1 so n = (fac**(e//E)*r)**E
544 m = n//fac**e
545 rE = perfect_power(m, candidates=divisors(e, generator=True))
546 if not rE:
547 return False
548 else:
549 r, E = rE
550 r, e = fac**(e//E)*r, E
551 if not big:
552 e0 = primefactors(e)
553 if e0[0] != e:
554 r, e = r**(e//e0[0]), e0[0]
555 return r, e
557 # Weed out downright impossible candidates
558 if logn/e < 40:
559 b = 2.0**(logn/e)
560 if abs(int(b + 0.5) - b) > 0.01:
561 continue
563 # now see if the plausible e makes a perfect power
564 r, exact = integer_nthroot(n, e)
565 if exact:
566 if big:
567 m = perfect_power(r, big=big, factor=factor)
568 if m:
569 r, e = m[0], e*m[1]
570 return int(r), e
572 return False
575def pollard_rho(n, s=2, a=1, retries=5, seed=1234, max_steps=None, F=None):
576 r"""
577 Use Pollard's rho method to try to extract a nontrivial factor
578 of ``n``. The returned factor may be a composite number. If no
579 factor is found, ``None`` is returned.
581 The algorithm generates pseudo-random values of x with a generator
582 function, replacing x with F(x). If F is not supplied then the
583 function x**2 + ``a`` is used. The first value supplied to F(x) is ``s``.
584 Upon failure (if ``retries`` is > 0) a new ``a`` and ``s`` will be
585 supplied; the ``a`` will be ignored if F was supplied.
587 The sequence of numbers generated by such functions generally have a
588 a lead-up to some number and then loop around back to that number and
589 begin to repeat the sequence, e.g. 1, 2, 3, 4, 5, 3, 4, 5 -- this leader
590 and loop look a bit like the Greek letter rho, and thus the name, 'rho'.
592 For a given function, very different leader-loop values can be obtained
593 so it is a good idea to allow for retries:
595 >>> from sympy.ntheory.generate import cycle_length
596 >>> n = 16843009
597 >>> F = lambda x:(2048*pow(x, 2, n) + 32767) % n
598 >>> for s in range(5):
599 ... print('loop length = %4i; leader length = %3i' % next(cycle_length(F, s)))
600 ...
601 loop length = 2489; leader length = 42
602 loop length = 78; leader length = 120
603 loop length = 1482; leader length = 99
604 loop length = 1482; leader length = 285
605 loop length = 1482; leader length = 100
607 Here is an explicit example where there is a two element leadup to
608 a sequence of 3 numbers (11, 14, 4) that then repeat:
610 >>> x=2
611 >>> for i in range(9):
612 ... x=(x**2+12)%17
613 ... print(x)
614 ...
615 16
616 13
617 11
618 14
619 4
620 11
621 14
622 4
623 11
624 >>> next(cycle_length(lambda x: (x**2+12)%17, 2))
625 (3, 2)
626 >>> list(cycle_length(lambda x: (x**2+12)%17, 2, values=True))
627 [16, 13, 11, 14, 4]
629 Instead of checking the differences of all generated values for a gcd
630 with n, only the kth and 2*kth numbers are checked, e.g. 1st and 2nd,
631 2nd and 4th, 3rd and 6th until it has been detected that the loop has been
632 traversed. Loops may be many thousands of steps long before rho finds a
633 factor or reports failure. If ``max_steps`` is specified, the iteration
634 is cancelled with a failure after the specified number of steps.
636 Examples
637 ========
639 >>> from sympy import pollard_rho
640 >>> n=16843009
641 >>> F=lambda x:(2048*pow(x,2,n) + 32767) % n
642 >>> pollard_rho(n, F=F)
643 257
645 Use the default setting with a bad value of ``a`` and no retries:
647 >>> pollard_rho(n, a=n-2, retries=0)
649 If retries is > 0 then perhaps the problem will correct itself when
650 new values are generated for a:
652 >>> pollard_rho(n, a=n-2, retries=1)
653 257
655 References
656 ==========
658 .. [1] Richard Crandall & Carl Pomerance (2005), "Prime Numbers:
659 A Computational Perspective", Springer, 2nd edition, 229-231
661 """
662 n = int(n)
663 if n < 5:
664 raise ValueError('pollard_rho should receive n > 4')
665 prng = random.Random(seed + retries)
666 V = s
667 for i in range(retries + 1):
668 U = V
669 if not F:
670 F = lambda x: (pow(x, 2, n) + a) % n
671 j = 0
672 while 1:
673 if max_steps and (j > max_steps):
674 break
675 j += 1
676 U = F(U)
677 V = F(F(V)) # V is 2x further along than U
678 g = igcd(U - V, n)
679 if g == 1:
680 continue
681 if g == n:
682 break
683 return int(g)
684 V = prng.randint(0, n - 1)
685 a = prng.randint(1, n - 3) # for x**2 + a, a%n should not be 0 or -2
686 F = None
687 return None
690def pollard_pm1(n, B=10, a=2, retries=0, seed=1234):
691 """
692 Use Pollard's p-1 method to try to extract a nontrivial factor
693 of ``n``. Either a divisor (perhaps composite) or ``None`` is returned.
695 The value of ``a`` is the base that is used in the test gcd(a**M - 1, n).
696 The default is 2. If ``retries`` > 0 then if no factor is found after the
697 first attempt, a new ``a`` will be generated randomly (using the ``seed``)
698 and the process repeated.
700 Note: the value of M is lcm(1..B) = reduce(ilcm, range(2, B + 1)).
702 A search is made for factors next to even numbers having a power smoothness
703 less than ``B``. Choosing a larger B increases the likelihood of finding a
704 larger factor but takes longer. Whether a factor of n is found or not
705 depends on ``a`` and the power smoothness of the even number just less than
706 the factor p (hence the name p - 1).
708 Although some discussion of what constitutes a good ``a`` some
709 descriptions are hard to interpret. At the modular.math site referenced
710 below it is stated that if gcd(a**M - 1, n) = N then a**M % q**r is 1
711 for every prime power divisor of N. But consider the following:
713 >>> from sympy.ntheory.factor_ import smoothness_p, pollard_pm1
714 >>> n=257*1009
715 >>> smoothness_p(n)
716 (-1, [(257, (1, 2, 256)), (1009, (1, 7, 16))])
718 So we should (and can) find a root with B=16:
720 >>> pollard_pm1(n, B=16, a=3)
721 1009
723 If we attempt to increase B to 256 we find that it does not work:
725 >>> pollard_pm1(n, B=256)
726 >>>
728 But if the value of ``a`` is changed we find that only multiples of
729 257 work, e.g.:
731 >>> pollard_pm1(n, B=256, a=257)
732 1009
734 Checking different ``a`` values shows that all the ones that did not
735 work had a gcd value not equal to ``n`` but equal to one of the
736 factors:
738 >>> from sympy import ilcm, igcd, factorint, Pow
739 >>> M = 1
740 >>> for i in range(2, 256):
741 ... M = ilcm(M, i)
742 ...
743 >>> set([igcd(pow(a, M, n) - 1, n) for a in range(2, 256) if
744 ... igcd(pow(a, M, n) - 1, n) != n])
745 {1009}
747 But does aM % d for every divisor of n give 1?
749 >>> aM = pow(255, M, n)
750 >>> [(d, aM%Pow(*d.args)) for d in factorint(n, visual=True).args]
751 [(257**1, 1), (1009**1, 1)]
753 No, only one of them. So perhaps the principle is that a root will
754 be found for a given value of B provided that:
756 1) the power smoothness of the p - 1 value next to the root
757 does not exceed B
758 2) a**M % p != 1 for any of the divisors of n.
760 By trying more than one ``a`` it is possible that one of them
761 will yield a factor.
763 Examples
764 ========
766 With the default smoothness bound, this number cannot be cracked:
768 >>> from sympy.ntheory import pollard_pm1
769 >>> pollard_pm1(21477639576571)
771 Increasing the smoothness bound helps:
773 >>> pollard_pm1(21477639576571, B=2000)
774 4410317
776 Looking at the smoothness of the factors of this number we find:
778 >>> from sympy.ntheory.factor_ import smoothness_p, factorint
779 >>> print(smoothness_p(21477639576571, visual=1))
780 p**i=4410317**1 has p-1 B=1787, B-pow=1787
781 p**i=4869863**1 has p-1 B=2434931, B-pow=2434931
783 The B and B-pow are the same for the p - 1 factorizations of the divisors
784 because those factorizations had a very large prime factor:
786 >>> factorint(4410317 - 1)
787 {2: 2, 617: 1, 1787: 1}
788 >>> factorint(4869863-1)
789 {2: 1, 2434931: 1}
791 Note that until B reaches the B-pow value of 1787, the number is not cracked;
793 >>> pollard_pm1(21477639576571, B=1786)
794 >>> pollard_pm1(21477639576571, B=1787)
795 4410317
797 The B value has to do with the factors of the number next to the divisor,
798 not the divisors themselves. A worst case scenario is that the number next
799 to the factor p has a large prime divisisor or is a perfect power. If these
800 conditions apply then the power-smoothness will be about p/2 or p. The more
801 realistic is that there will be a large prime factor next to p requiring
802 a B value on the order of p/2. Although primes may have been searched for
803 up to this level, the p/2 is a factor of p - 1, something that we do not
804 know. The modular.math reference below states that 15% of numbers in the
805 range of 10**15 to 15**15 + 10**4 are 10**6 power smooth so a B of 10**6
806 will fail 85% of the time in that range. From 10**8 to 10**8 + 10**3 the
807 percentages are nearly reversed...but in that range the simple trial
808 division is quite fast.
810 References
811 ==========
813 .. [1] Richard Crandall & Carl Pomerance (2005), "Prime Numbers:
814 A Computational Perspective", Springer, 2nd edition, 236-238
815 .. [2] https://web.archive.org/web/20150716201437/http://modular.math.washington.edu/edu/2007/spring/ent/ent-html/node81.html
816 .. [3] https://www.cs.toronto.edu/~yuvalf/Factorization.pdf
817 """
819 n = int(n)
820 if n < 4 or B < 3:
821 raise ValueError('pollard_pm1 should receive n > 3 and B > 2')
822 prng = random.Random(seed + B)
824 # computing a**lcm(1,2,3,..B) % n for B > 2
825 # it looks weird, but it's right: primes run [2, B]
826 # and the answer's not right until the loop is done.
827 for i in range(retries + 1):
828 aM = a
829 for p in sieve.primerange(2, B + 1):
830 e = int(math.log(B, p))
831 aM = pow(aM, pow(p, e), n)
832 g = igcd(aM - 1, n)
833 if 1 < g < n:
834 return int(g)
836 # get a new a:
837 # since the exponent, lcm(1..B), is even, if we allow 'a' to be 'n-1'
838 # then (n - 1)**even % n will be 1 which will give a g of 0 and 1 will
839 # give a zero, too, so we set the range as [2, n-2]. Some references
840 # say 'a' should be coprime to n, but either will detect factors.
841 a = prng.randint(2, n - 2)
844def _trial(factors, n, candidates, verbose=False):
845 """
846 Helper function for integer factorization. Trial factors ``n`
847 against all integers given in the sequence ``candidates``
848 and updates the dict ``factors`` in-place. Returns the reduced
849 value of ``n`` and a flag indicating whether any factors were found.
850 """
851 if verbose:
852 factors0 = list(factors.keys())
853 nfactors = len(factors)
854 for d in candidates:
855 if n % d == 0:
856 m = multiplicity(d, n)
857 n //= d**m
858 factors[d] = m
859 if verbose:
860 for k in sorted(set(factors).difference(set(factors0))):
861 print(factor_msg % (k, factors[k]))
862 return int(n), len(factors) != nfactors
865def _check_termination(factors, n, limitp1, use_trial, use_rho, use_pm1,
866 verbose):
867 """
868 Helper function for integer factorization. Checks if ``n``
869 is a prime or a perfect power, and in those cases updates
870 the factorization and raises ``StopIteration``.
871 """
873 if verbose:
874 print('Check for termination')
876 # since we've already been factoring there is no need to do
877 # simultaneous factoring with the power check
878 p = perfect_power(n, factor=False)
879 if p is not False:
880 base, exp = p
881 if limitp1:
882 limit = limitp1 - 1
883 else:
884 limit = limitp1
885 facs = factorint(base, limit, use_trial, use_rho, use_pm1,
886 verbose=False)
887 for b, e in facs.items():
888 if verbose:
889 print(factor_msg % (b, e))
890 factors[b] = exp*e
891 raise StopIteration
893 if isprime(n):
894 factors[int(n)] = 1
895 raise StopIteration
897 if n == 1:
898 raise StopIteration
900trial_int_msg = "Trial division with ints [%i ... %i] and fail_max=%i"
901trial_msg = "Trial division with primes [%i ... %i]"
902rho_msg = "Pollard's rho with retries %i, max_steps %i and seed %i"
903pm1_msg = "Pollard's p-1 with smoothness bound %i and seed %i"
904ecm_msg = "Elliptic Curve with B1 bound %i, B2 bound %i, num_curves %i"
905factor_msg = '\t%i ** %i'
906fermat_msg = 'Close factors satisying Fermat condition found.'
907complete_msg = 'Factorization is complete.'
910def _factorint_small(factors, n, limit, fail_max):
911 """
912 Return the value of n and either a 0 (indicating that factorization up
913 to the limit was complete) or else the next near-prime that would have
914 been tested.
916 Factoring stops if there are fail_max unsuccessful tests in a row.
918 If factors of n were found they will be in the factors dictionary as
919 {factor: multiplicity} and the returned value of n will have had those
920 factors removed. The factors dictionary is modified in-place.
922 """
924 def done(n, d):
925 """return n, d if the sqrt(n) was not reached yet, else
926 n, 0 indicating that factoring is done.
927 """
928 if d*d <= n:
929 return n, d
930 return n, 0
932 d = 2
933 m = trailing(n)
934 if m:
935 factors[d] = m
936 n >>= m
937 d = 3
938 if limit < d:
939 if n > 1:
940 factors[n] = 1
941 return done(n, d)
942 # reduce
943 m = 0
944 while n % d == 0:
945 n //= d
946 m += 1
947 if m == 20:
948 mm = multiplicity(d, n)
949 m += mm
950 n //= d**mm
951 break
952 if m:
953 factors[d] = m
955 # when d*d exceeds maxx or n we are done; if limit**2 is greater
956 # than n then maxx is set to zero so the value of n will flag the finish
957 if limit*limit > n:
958 maxx = 0
959 else:
960 maxx = limit*limit
962 dd = maxx or n
963 d = 5
964 fails = 0
965 while fails < fail_max:
966 if d*d > dd:
967 break
968 # d = 6*i - 1
969 # reduce
970 m = 0
971 while n % d == 0:
972 n //= d
973 m += 1
974 if m == 20:
975 mm = multiplicity(d, n)
976 m += mm
977 n //= d**mm
978 break
979 if m:
980 factors[d] = m
981 dd = maxx or n
982 fails = 0
983 else:
984 fails += 1
985 d += 2
986 if d*d > dd:
987 break
988 # d = 6*i - 1
989 # reduce
990 m = 0
991 while n % d == 0:
992 n //= d
993 m += 1
994 if m == 20:
995 mm = multiplicity(d, n)
996 m += mm
997 n //= d**mm
998 break
999 if m:
1000 factors[d] = m
1001 dd = maxx or n
1002 fails = 0
1003 else:
1004 fails += 1
1005 # d = 6*(i + 1) - 1
1006 d += 4
1008 return done(n, d)
1011def factorint(n, limit=None, use_trial=True, use_rho=True, use_pm1=True,
1012 use_ecm=True, verbose=False, visual=None, multiple=False):
1013 r"""
1014 Given a positive integer ``n``, ``factorint(n)`` returns a dict containing
1015 the prime factors of ``n`` as keys and their respective multiplicities
1016 as values. For example:
1018 >>> from sympy.ntheory import factorint
1019 >>> factorint(2000) # 2000 = (2**4) * (5**3)
1020 {2: 4, 5: 3}
1021 >>> factorint(65537) # This number is prime
1022 {65537: 1}
1024 For input less than 2, factorint behaves as follows:
1026 - ``factorint(1)`` returns the empty factorization, ``{}``
1027 - ``factorint(0)`` returns ``{0:1}``
1028 - ``factorint(-n)`` adds ``-1:1`` to the factors and then factors ``n``
1030 Partial Factorization:
1032 If ``limit`` (> 3) is specified, the search is stopped after performing
1033 trial division up to (and including) the limit (or taking a
1034 corresponding number of rho/p-1 steps). This is useful if one has
1035 a large number and only is interested in finding small factors (if
1036 any). Note that setting a limit does not prevent larger factors
1037 from being found early; it simply means that the largest factor may
1038 be composite. Since checking for perfect power is relatively cheap, it is
1039 done regardless of the limit setting.
1041 This number, for example, has two small factors and a huge
1042 semi-prime factor that cannot be reduced easily:
1044 >>> from sympy.ntheory import isprime
1045 >>> a = 1407633717262338957430697921446883
1046 >>> f = factorint(a, limit=10000)
1047 >>> f == {991: 1, int(202916782076162456022877024859): 1, 7: 1}
1048 True
1049 >>> isprime(max(f))
1050 False
1052 This number has a small factor and a residual perfect power whose
1053 base is greater than the limit:
1055 >>> factorint(3*101**7, limit=5)
1056 {3: 1, 101: 7}
1058 List of Factors:
1060 If ``multiple`` is set to ``True`` then a list containing the
1061 prime factors including multiplicities is returned.
1063 >>> factorint(24, multiple=True)
1064 [2, 2, 2, 3]
1066 Visual Factorization:
1068 If ``visual`` is set to ``True``, then it will return a visual
1069 factorization of the integer. For example:
1071 >>> from sympy import pprint
1072 >>> pprint(factorint(4200, visual=True))
1073 3 1 2 1
1074 2 *3 *5 *7
1076 Note that this is achieved by using the evaluate=False flag in Mul
1077 and Pow. If you do other manipulations with an expression where
1078 evaluate=False, it may evaluate. Therefore, you should use the
1079 visual option only for visualization, and use the normal dictionary
1080 returned by visual=False if you want to perform operations on the
1081 factors.
1083 You can easily switch between the two forms by sending them back to
1084 factorint:
1086 >>> from sympy import Mul
1087 >>> regular = factorint(1764); regular
1088 {2: 2, 3: 2, 7: 2}
1089 >>> pprint(factorint(regular))
1090 2 2 2
1091 2 *3 *7
1093 >>> visual = factorint(1764, visual=True); pprint(visual)
1094 2 2 2
1095 2 *3 *7
1096 >>> print(factorint(visual))
1097 {2: 2, 3: 2, 7: 2}
1099 If you want to send a number to be factored in a partially factored form
1100 you can do so with a dictionary or unevaluated expression:
1102 >>> factorint(factorint({4: 2, 12: 3})) # twice to toggle to dict form
1103 {2: 10, 3: 3}
1104 >>> factorint(Mul(4, 12, evaluate=False))
1105 {2: 4, 3: 1}
1107 The table of the output logic is:
1109 ====== ====== ======= =======
1110 Visual
1111 ------ ----------------------
1112 Input True False other
1113 ====== ====== ======= =======
1114 dict mul dict mul
1115 n mul dict dict
1116 mul mul dict dict
1117 ====== ====== ======= =======
1119 Notes
1120 =====
1122 Algorithm:
1124 The function switches between multiple algorithms. Trial division
1125 quickly finds small factors (of the order 1-5 digits), and finds
1126 all large factors if given enough time. The Pollard rho and p-1
1127 algorithms are used to find large factors ahead of time; they
1128 will often find factors of the order of 10 digits within a few
1129 seconds:
1131 >>> factors = factorint(12345678910111213141516)
1132 >>> for base, exp in sorted(factors.items()):
1133 ... print('%s %s' % (base, exp))
1134 ...
1135 2 2
1136 2507191691 1
1137 1231026625769 1
1139 Any of these methods can optionally be disabled with the following
1140 boolean parameters:
1142 - ``use_trial``: Toggle use of trial division
1143 - ``use_rho``: Toggle use of Pollard's rho method
1144 - ``use_pm1``: Toggle use of Pollard's p-1 method
1146 ``factorint`` also periodically checks if the remaining part is
1147 a prime number or a perfect power, and in those cases stops.
1149 For unevaluated factorial, it uses Legendre's formula(theorem).
1152 If ``verbose`` is set to ``True``, detailed progress is printed.
1154 See Also
1155 ========
1157 smoothness, smoothness_p, divisors
1159 """
1160 if isinstance(n, Dict):
1161 n = dict(n)
1162 if multiple:
1163 fac = factorint(n, limit=limit, use_trial=use_trial,
1164 use_rho=use_rho, use_pm1=use_pm1,
1165 verbose=verbose, visual=False, multiple=False)
1166 factorlist = sum(([p] * fac[p] if fac[p] > 0 else [S.One/p]*(-fac[p])
1167 for p in sorted(fac)), [])
1168 return factorlist
1170 factordict = {}
1171 if visual and not isinstance(n, (Mul, dict)):
1172 factordict = factorint(n, limit=limit, use_trial=use_trial,
1173 use_rho=use_rho, use_pm1=use_pm1,
1174 verbose=verbose, visual=False)
1175 elif isinstance(n, Mul):
1176 factordict = {int(k): int(v) for k, v in
1177 n.as_powers_dict().items()}
1178 elif isinstance(n, dict):
1179 factordict = n
1180 if factordict and isinstance(n, (Mul, dict)):
1181 # check it
1182 for key in list(factordict.keys()):
1183 if isprime(key):
1184 continue
1185 e = factordict.pop(key)
1186 d = factorint(key, limit=limit, use_trial=use_trial, use_rho=use_rho,
1187 use_pm1=use_pm1, verbose=verbose, visual=False)
1188 for k, v in d.items():
1189 if k in factordict:
1190 factordict[k] += v*e
1191 else:
1192 factordict[k] = v*e
1193 if visual or (type(n) is dict and
1194 visual is not True and
1195 visual is not False):
1196 if factordict == {}:
1197 return S.One
1198 if -1 in factordict:
1199 factordict.pop(-1)
1200 args = [S.NegativeOne]
1201 else:
1202 args = []
1203 args.extend([Pow(*i, evaluate=False)
1204 for i in sorted(factordict.items())])
1205 return Mul(*args, evaluate=False)
1206 elif isinstance(n, (dict, Mul)):
1207 return factordict
1209 assert use_trial or use_rho or use_pm1 or use_ecm
1211 from sympy.functions.combinatorial.factorials import factorial
1212 if isinstance(n, factorial):
1213 x = as_int(n.args[0])
1214 if x >= 20:
1215 factors = {}
1216 m = 2 # to initialize the if condition below
1217 for p in sieve.primerange(2, x + 1):
1218 if m > 1:
1219 m, q = 0, x // p
1220 while q != 0:
1221 m += q
1222 q //= p
1223 factors[p] = m
1224 if factors and verbose:
1225 for k in sorted(factors):
1226 print(factor_msg % (k, factors[k]))
1227 if verbose:
1228 print(complete_msg)
1229 return factors
1230 else:
1231 # if n < 20!, direct computation is faster
1232 # since it uses a lookup table
1233 n = n.func(x)
1235 n = as_int(n)
1236 if limit:
1237 limit = int(limit)
1238 use_ecm = False
1240 # special cases
1241 if n < 0:
1242 factors = factorint(
1243 -n, limit=limit, use_trial=use_trial, use_rho=use_rho,
1244 use_pm1=use_pm1, verbose=verbose, visual=False)
1245 factors[-1] = 1
1246 return factors
1248 if limit and limit < 2:
1249 if n == 1:
1250 return {}
1251 return {n: 1}
1252 elif n < 10:
1253 # doing this we are assured of getting a limit > 2
1254 # when we have to compute it later
1255 return [{0: 1}, {}, {2: 1}, {3: 1}, {2: 2}, {5: 1},
1256 {2: 1, 3: 1}, {7: 1}, {2: 3}, {3: 2}][n]
1258 factors = {}
1260 # do simplistic factorization
1261 if verbose:
1262 sn = str(n)
1263 if len(sn) > 50:
1264 print('Factoring %s' % sn[:5] + \
1265 '..(%i other digits)..' % (len(sn) - 10) + sn[-5:])
1266 else:
1267 print('Factoring', n)
1269 if use_trial:
1270 # this is the preliminary factorization for small factors
1271 small = 2**15
1272 fail_max = 600
1273 small = min(small, limit or small)
1274 if verbose:
1275 print(trial_int_msg % (2, small, fail_max))
1276 n, next_p = _factorint_small(factors, n, small, fail_max)
1277 else:
1278 next_p = 2
1279 if factors and verbose:
1280 for k in sorted(factors):
1281 print(factor_msg % (k, factors[k]))
1282 if next_p == 0:
1283 if n > 1:
1284 factors[int(n)] = 1
1285 if verbose:
1286 print(complete_msg)
1287 return factors
1289 # continue with more advanced factorization methods
1291 # first check if the simplistic run didn't finish
1292 # because of the limit and check for a perfect
1293 # power before exiting
1294 try:
1295 if limit and next_p > limit:
1296 if verbose:
1297 print('Exceeded limit:', limit)
1299 _check_termination(factors, n, limit, use_trial, use_rho, use_pm1,
1300 verbose)
1302 if n > 1:
1303 factors[int(n)] = 1
1304 return factors
1305 else:
1306 # Before quitting (or continuing on)...
1308 # ...do a Fermat test since it's so easy and we need the
1309 # square root anyway. Finding 2 factors is easy if they are
1310 # "close enough." This is the big root equivalent of dividing by
1311 # 2, 3, 5.
1312 sqrt_n = integer_nthroot(n, 2)[0]
1313 a = sqrt_n + 1
1314 a2 = a**2
1315 b2 = a2 - n
1316 for i in range(3):
1317 b, fermat = integer_nthroot(b2, 2)
1318 if fermat:
1319 break
1320 b2 += 2*a + 1 # equiv to (a + 1)**2 - n
1321 a += 1
1322 if fermat:
1323 if verbose:
1324 print(fermat_msg)
1325 if limit:
1326 limit -= 1
1327 for r in [a - b, a + b]:
1328 facs = factorint(r, limit=limit, use_trial=use_trial,
1329 use_rho=use_rho, use_pm1=use_pm1,
1330 verbose=verbose)
1331 for k, v in facs.items():
1332 factors[k] = factors.get(k, 0) + v
1333 raise StopIteration
1335 # ...see if factorization can be terminated
1336 _check_termination(factors, n, limit, use_trial, use_rho, use_pm1,
1337 verbose)
1339 except StopIteration:
1340 if verbose:
1341 print(complete_msg)
1342 return factors
1344 # these are the limits for trial division which will
1345 # be attempted in parallel with pollard methods
1346 low, high = next_p, 2*next_p
1348 limit = limit or sqrt_n
1349 # add 1 to make sure limit is reached in primerange calls
1350 limit += 1
1351 iteration = 0
1352 while 1:
1354 try:
1355 high_ = high
1356 if limit < high_:
1357 high_ = limit
1359 # Trial division
1360 if use_trial:
1361 if verbose:
1362 print(trial_msg % (low, high_))
1363 ps = sieve.primerange(low, high_)
1364 n, found_trial = _trial(factors, n, ps, verbose)
1365 if found_trial:
1366 _check_termination(factors, n, limit, use_trial, use_rho,
1367 use_pm1, verbose)
1368 else:
1369 found_trial = False
1371 if high > limit:
1372 if verbose:
1373 print('Exceeded limit:', limit)
1374 if n > 1:
1375 factors[int(n)] = 1
1376 raise StopIteration
1378 # Only used advanced methods when no small factors were found
1379 if not found_trial:
1380 if (use_pm1 or use_rho):
1381 high_root = max(int(math.log(high_**0.7)), low, 3)
1383 # Pollard p-1
1384 if use_pm1:
1385 if verbose:
1386 print(pm1_msg % (high_root, high_))
1387 c = pollard_pm1(n, B=high_root, seed=high_)
1388 if c:
1389 # factor it and let _trial do the update
1390 ps = factorint(c, limit=limit - 1,
1391 use_trial=use_trial,
1392 use_rho=use_rho,
1393 use_pm1=use_pm1,
1394 use_ecm=use_ecm,
1395 verbose=verbose)
1396 n, _ = _trial(factors, n, ps, verbose=False)
1397 _check_termination(factors, n, limit, use_trial,
1398 use_rho, use_pm1, verbose)
1400 # Pollard rho
1401 if use_rho:
1402 max_steps = high_root
1403 if verbose:
1404 print(rho_msg % (1, max_steps, high_))
1405 c = pollard_rho(n, retries=1, max_steps=max_steps,
1406 seed=high_)
1407 if c:
1408 # factor it and let _trial do the update
1409 ps = factorint(c, limit=limit - 1,
1410 use_trial=use_trial,
1411 use_rho=use_rho,
1412 use_pm1=use_pm1,
1413 use_ecm=use_ecm,
1414 verbose=verbose)
1415 n, _ = _trial(factors, n, ps, verbose=False)
1416 _check_termination(factors, n, limit, use_trial,
1417 use_rho, use_pm1, verbose)
1419 except StopIteration:
1420 if verbose:
1421 print(complete_msg)
1422 return factors
1423 #Use subexponential algorithms if use_ecm
1424 #Use pollard algorithms for finding small factors for 3 iterations
1425 #if after small factors the number of digits of n is >= 20 then use ecm
1426 iteration += 1
1427 if use_ecm and iteration >= 3 and len(str(n)) >= 25:
1428 break
1429 low, high = high, high*2
1430 B1 = 10000
1431 B2 = 100*B1
1432 num_curves = 50
1433 while(1):
1434 if verbose:
1435 print(ecm_msg % (B1, B2, num_curves))
1436 while(1):
1437 try:
1438 factor = _ecm_one_factor(n, B1, B2, num_curves)
1439 ps = factorint(factor, limit=limit - 1,
1440 use_trial=use_trial,
1441 use_rho=use_rho,
1442 use_pm1=use_pm1,
1443 use_ecm=use_ecm,
1444 verbose=verbose)
1445 n, _ = _trial(factors, n, ps, verbose=False)
1446 _check_termination(factors, n, limit, use_trial,
1447 use_rho, use_pm1, verbose)
1448 except ValueError:
1449 break
1450 except StopIteration:
1451 if verbose:
1452 print(complete_msg)
1453 return factors
1454 B1 *= 5
1455 B2 = 100*B1
1456 num_curves *= 4
1459def factorrat(rat, limit=None, use_trial=True, use_rho=True, use_pm1=True,
1460 verbose=False, visual=None, multiple=False):
1461 r"""
1462 Given a Rational ``r``, ``factorrat(r)`` returns a dict containing
1463 the prime factors of ``r`` as keys and their respective multiplicities
1464 as values. For example:
1466 >>> from sympy import factorrat, S
1467 >>> factorrat(S(8)/9) # 8/9 = (2**3) * (3**-2)
1468 {2: 3, 3: -2}
1469 >>> factorrat(S(-1)/987) # -1/789 = -1 * (3**-1) * (7**-1) * (47**-1)
1470 {-1: 1, 3: -1, 7: -1, 47: -1}
1472 Please see the docstring for ``factorint`` for detailed explanations
1473 and examples of the following keywords:
1475 - ``limit``: Integer limit up to which trial division is done
1476 - ``use_trial``: Toggle use of trial division
1477 - ``use_rho``: Toggle use of Pollard's rho method
1478 - ``use_pm1``: Toggle use of Pollard's p-1 method
1479 - ``verbose``: Toggle detailed printing of progress
1480 - ``multiple``: Toggle returning a list of factors or dict
1481 - ``visual``: Toggle product form of output
1482 """
1483 if multiple:
1484 fac = factorrat(rat, limit=limit, use_trial=use_trial,
1485 use_rho=use_rho, use_pm1=use_pm1,
1486 verbose=verbose, visual=False, multiple=False)
1487 factorlist = sum(([p] * fac[p] if fac[p] > 0 else [S.One/p]*(-fac[p])
1488 for p, _ in sorted(fac.items(),
1489 key=lambda elem: elem[0]
1490 if elem[1] > 0
1491 else 1/elem[0])), [])
1492 return factorlist
1494 f = factorint(rat.p, limit=limit, use_trial=use_trial,
1495 use_rho=use_rho, use_pm1=use_pm1,
1496 verbose=verbose).copy()
1497 f = defaultdict(int, f)
1498 for p, e in factorint(rat.q, limit=limit,
1499 use_trial=use_trial,
1500 use_rho=use_rho,
1501 use_pm1=use_pm1,
1502 verbose=verbose).items():
1503 f[p] += -e
1505 if len(f) > 1 and 1 in f:
1506 del f[1]
1507 if not visual:
1508 return dict(f)
1509 else:
1510 if -1 in f:
1511 f.pop(-1)
1512 args = [S.NegativeOne]
1513 else:
1514 args = []
1515 args.extend([Pow(*i, evaluate=False)
1516 for i in sorted(f.items())])
1517 return Mul(*args, evaluate=False)
1521def primefactors(n, limit=None, verbose=False):
1522 """Return a sorted list of n's prime factors, ignoring multiplicity
1523 and any composite factor that remains if the limit was set too low
1524 for complete factorization. Unlike factorint(), primefactors() does
1525 not return -1 or 0.
1527 Examples
1528 ========
1530 >>> from sympy.ntheory import primefactors, factorint, isprime
1531 >>> primefactors(6)
1532 [2, 3]
1533 >>> primefactors(-5)
1534 [5]
1536 >>> sorted(factorint(123456).items())
1537 [(2, 6), (3, 1), (643, 1)]
1538 >>> primefactors(123456)
1539 [2, 3, 643]
1541 >>> sorted(factorint(10000000001, limit=200).items())
1542 [(101, 1), (99009901, 1)]
1543 >>> isprime(99009901)
1544 False
1545 >>> primefactors(10000000001, limit=300)
1546 [101]
1548 See Also
1549 ========
1551 divisors
1552 """
1553 n = int(n)
1554 factors = sorted(factorint(n, limit=limit, verbose=verbose).keys())
1555 s = [f for f in factors[:-1:] if f not in [-1, 0, 1]]
1556 if factors and isprime(factors[-1]):
1557 s += [factors[-1]]
1558 return s
1561def _divisors(n, proper=False):
1562 """Helper function for divisors which generates the divisors."""
1564 factordict = factorint(n)
1565 ps = sorted(factordict.keys())
1567 def rec_gen(n=0):
1568 if n == len(ps):
1569 yield 1
1570 else:
1571 pows = [1]
1572 for j in range(factordict[ps[n]]):
1573 pows.append(pows[-1] * ps[n])
1574 for q in rec_gen(n + 1):
1575 for p in pows:
1576 yield p * q
1578 if proper:
1579 for p in rec_gen():
1580 if p != n:
1581 yield p
1582 else:
1583 yield from rec_gen()
1586def divisors(n, generator=False, proper=False):
1587 r"""
1588 Return all divisors of n sorted from 1..n by default.
1589 If generator is ``True`` an unordered generator is returned.
1591 The number of divisors of n can be quite large if there are many
1592 prime factors (counting repeated factors). If only the number of
1593 factors is desired use divisor_count(n).
1595 Examples
1596 ========
1598 >>> from sympy import divisors, divisor_count
1599 >>> divisors(24)
1600 [1, 2, 3, 4, 6, 8, 12, 24]
1601 >>> divisor_count(24)
1602 8
1604 >>> list(divisors(120, generator=True))
1605 [1, 2, 4, 8, 3, 6, 12, 24, 5, 10, 20, 40, 15, 30, 60, 120]
1607 Notes
1608 =====
1610 This is a slightly modified version of Tim Peters referenced at:
1611 https://stackoverflow.com/questions/1010381/python-factorization
1613 See Also
1614 ========
1616 primefactors, factorint, divisor_count
1617 """
1619 n = as_int(abs(n))
1620 if isprime(n):
1621 if proper:
1622 return [1]
1623 return [1, n]
1624 if n == 1:
1625 if proper:
1626 return []
1627 return [1]
1628 if n == 0:
1629 return []
1630 rv = _divisors(n, proper)
1631 if not generator:
1632 return sorted(rv)
1633 return rv
1636def divisor_count(n, modulus=1, proper=False):
1637 """
1638 Return the number of divisors of ``n``. If ``modulus`` is not 1 then only
1639 those that are divisible by ``modulus`` are counted. If ``proper`` is True
1640 then the divisor of ``n`` will not be counted.
1642 Examples
1643 ========
1645 >>> from sympy import divisor_count
1646 >>> divisor_count(6)
1647 4
1648 >>> divisor_count(6, 2)
1649 2
1650 >>> divisor_count(6, proper=True)
1651 3
1653 See Also
1654 ========
1656 factorint, divisors, totient, proper_divisor_count
1658 """
1660 if not modulus:
1661 return 0
1662 elif modulus != 1:
1663 n, r = divmod(n, modulus)
1664 if r:
1665 return 0
1666 if n == 0:
1667 return 0
1668 n = Mul(*[v + 1 for k, v in factorint(n).items() if k > 1])
1669 if n and proper:
1670 n -= 1
1671 return n
1674def proper_divisors(n, generator=False):
1675 """
1676 Return all divisors of n except n, sorted by default.
1677 If generator is ``True`` an unordered generator is returned.
1679 Examples
1680 ========
1682 >>> from sympy import proper_divisors, proper_divisor_count
1683 >>> proper_divisors(24)
1684 [1, 2, 3, 4, 6, 8, 12]
1685 >>> proper_divisor_count(24)
1686 7
1687 >>> list(proper_divisors(120, generator=True))
1688 [1, 2, 4, 8, 3, 6, 12, 24, 5, 10, 20, 40, 15, 30, 60]
1690 See Also
1691 ========
1693 factorint, divisors, proper_divisor_count
1695 """
1696 return divisors(n, generator=generator, proper=True)
1699def proper_divisor_count(n, modulus=1):
1700 """
1701 Return the number of proper divisors of ``n``.
1703 Examples
1704 ========
1706 >>> from sympy import proper_divisor_count
1707 >>> proper_divisor_count(6)
1708 3
1709 >>> proper_divisor_count(6, modulus=2)
1710 1
1712 See Also
1713 ========
1715 divisors, proper_divisors, divisor_count
1717 """
1718 return divisor_count(n, modulus=modulus, proper=True)
1721def _udivisors(n):
1722 """Helper function for udivisors which generates the unitary divisors."""
1724 factorpows = [p**e for p, e in factorint(n).items()]
1725 for i in range(2**len(factorpows)):
1726 d, j, k = 1, i, 0
1727 while j:
1728 if (j & 1):
1729 d *= factorpows[k]
1730 j >>= 1
1731 k += 1
1732 yield d
1735def udivisors(n, generator=False):
1736 r"""
1737 Return all unitary divisors of n sorted from 1..n by default.
1738 If generator is ``True`` an unordered generator is returned.
1740 The number of unitary divisors of n can be quite large if there are many
1741 prime factors. If only the number of unitary divisors is desired use
1742 udivisor_count(n).
1744 Examples
1745 ========
1747 >>> from sympy.ntheory.factor_ import udivisors, udivisor_count
1748 >>> udivisors(15)
1749 [1, 3, 5, 15]
1750 >>> udivisor_count(15)
1751 4
1753 >>> sorted(udivisors(120, generator=True))
1754 [1, 3, 5, 8, 15, 24, 40, 120]
1756 See Also
1757 ========
1759 primefactors, factorint, divisors, divisor_count, udivisor_count
1761 References
1762 ==========
1764 .. [1] https://en.wikipedia.org/wiki/Unitary_divisor
1765 .. [2] https://mathworld.wolfram.com/UnitaryDivisor.html
1767 """
1769 n = as_int(abs(n))
1770 if isprime(n):
1771 return [1, n]
1772 if n == 1:
1773 return [1]
1774 if n == 0:
1775 return []
1776 rv = _udivisors(n)
1777 if not generator:
1778 return sorted(rv)
1779 return rv
1782def udivisor_count(n):
1783 """
1784 Return the number of unitary divisors of ``n``.
1786 Parameters
1787 ==========
1789 n : integer
1791 Examples
1792 ========
1794 >>> from sympy.ntheory.factor_ import udivisor_count
1795 >>> udivisor_count(120)
1796 8
1798 See Also
1799 ========
1801 factorint, divisors, udivisors, divisor_count, totient
1803 References
1804 ==========
1806 .. [1] https://mathworld.wolfram.com/UnitaryDivisorFunction.html
1808 """
1810 if n == 0:
1811 return 0
1812 return 2**len([p for p in factorint(n) if p > 1])
1815def _antidivisors(n):
1816 """Helper function for antidivisors which generates the antidivisors."""
1818 for d in _divisors(n):
1819 y = 2*d
1820 if n > y and n % y:
1821 yield y
1822 for d in _divisors(2*n-1):
1823 if n > d >= 2 and n % d:
1824 yield d
1825 for d in _divisors(2*n+1):
1826 if n > d >= 2 and n % d:
1827 yield d
1830def antidivisors(n, generator=False):
1831 r"""
1832 Return all antidivisors of n sorted from 1..n by default.
1834 Antidivisors [1]_ of n are numbers that do not divide n by the largest
1835 possible margin. If generator is True an unordered generator is returned.
1837 Examples
1838 ========
1840 >>> from sympy.ntheory.factor_ import antidivisors
1841 >>> antidivisors(24)
1842 [7, 16]
1844 >>> sorted(antidivisors(128, generator=True))
1845 [3, 5, 15, 17, 51, 85]
1847 See Also
1848 ========
1850 primefactors, factorint, divisors, divisor_count, antidivisor_count
1852 References
1853 ==========
1855 .. [1] definition is described in https://oeis.org/A066272/a066272a.html
1857 """
1859 n = as_int(abs(n))
1860 if n <= 2:
1861 return []
1862 rv = _antidivisors(n)
1863 if not generator:
1864 return sorted(rv)
1865 return rv
1868def antidivisor_count(n):
1869 """
1870 Return the number of antidivisors [1]_ of ``n``.
1872 Parameters
1873 ==========
1875 n : integer
1877 Examples
1878 ========
1880 >>> from sympy.ntheory.factor_ import antidivisor_count
1881 >>> antidivisor_count(13)
1882 4
1883 >>> antidivisor_count(27)
1884 5
1886 See Also
1887 ========
1889 factorint, divisors, antidivisors, divisor_count, totient
1891 References
1892 ==========
1894 .. [1] formula from https://oeis.org/A066272
1896 """
1898 n = as_int(abs(n))
1899 if n <= 2:
1900 return 0
1901 return divisor_count(2*n - 1) + divisor_count(2*n + 1) + \
1902 divisor_count(n) - divisor_count(n, 2) - 5
1905class totient(Function):
1906 r"""
1907 Calculate the Euler totient function phi(n)
1909 ``totient(n)`` or `\phi(n)` is the number of positive integers `\leq` n
1910 that are relatively prime to n.
1912 Parameters
1913 ==========
1915 n : integer
1917 Examples
1918 ========
1920 >>> from sympy.ntheory import totient
1921 >>> totient(1)
1922 1
1923 >>> totient(25)
1924 20
1925 >>> totient(45) == totient(5)*totient(9)
1926 True
1928 See Also
1929 ========
1931 divisor_count
1933 References
1934 ==========
1936 .. [1] https://en.wikipedia.org/wiki/Euler%27s_totient_function
1937 .. [2] https://mathworld.wolfram.com/TotientFunction.html
1939 """
1940 @classmethod
1941 def eval(cls, n):
1942 if n.is_Integer:
1943 if n < 1:
1944 raise ValueError("n must be a positive integer")
1945 factors = factorint(n)
1946 return cls._from_factors(factors)
1947 elif not isinstance(n, Expr) or (n.is_integer is False) or (n.is_positive is False):
1948 raise ValueError("n must be a positive integer")
1950 def _eval_is_integer(self):
1951 return fuzzy_and([self.args[0].is_integer, self.args[0].is_positive])
1953 @classmethod
1954 def _from_distinct_primes(self, *args):
1955 """Subroutine to compute totient from the list of assumed
1956 distinct primes
1958 Examples
1959 ========
1961 >>> from sympy.ntheory.factor_ import totient
1962 >>> totient._from_distinct_primes(5, 7)
1963 24
1964 """
1965 return reduce(lambda i, j: i * (j-1), args, 1)
1967 @classmethod
1968 def _from_factors(self, factors):
1969 """Subroutine to compute totient from already-computed factors
1971 Examples
1972 ========
1974 >>> from sympy.ntheory.factor_ import totient
1975 >>> totient._from_factors({5: 2})
1976 20
1977 """
1978 t = 1
1979 for p, k in factors.items():
1980 t *= (p - 1) * p**(k - 1)
1981 return t
1984class reduced_totient(Function):
1985 r"""
1986 Calculate the Carmichael reduced totient function lambda(n)
1988 ``reduced_totient(n)`` or `\lambda(n)` is the smallest m > 0 such that
1989 `k^m \equiv 1 \mod n` for all k relatively prime to n.
1991 Examples
1992 ========
1994 >>> from sympy.ntheory import reduced_totient
1995 >>> reduced_totient(1)
1996 1
1997 >>> reduced_totient(8)
1998 2
1999 >>> reduced_totient(30)
2000 4
2002 See Also
2003 ========
2005 totient
2007 References
2008 ==========
2010 .. [1] https://en.wikipedia.org/wiki/Carmichael_function
2011 .. [2] https://mathworld.wolfram.com/CarmichaelFunction.html
2013 """
2014 @classmethod
2015 def eval(cls, n):
2016 if n.is_Integer:
2017 if n < 1:
2018 raise ValueError("n must be a positive integer")
2019 factors = factorint(n)
2020 return cls._from_factors(factors)
2022 @classmethod
2023 def _from_factors(self, factors):
2024 """Subroutine to compute totient from already-computed factors
2025 """
2026 t = 1
2027 for p, k in factors.items():
2028 if p == 2 and k > 2:
2029 t = ilcm(t, 2**(k - 2))
2030 else:
2031 t = ilcm(t, (p - 1) * p**(k - 1))
2032 return t
2034 @classmethod
2035 def _from_distinct_primes(self, *args):
2036 """Subroutine to compute totient from the list of assumed
2037 distinct primes
2038 """
2039 args = [p - 1 for p in args]
2040 return ilcm(*args)
2042 def _eval_is_integer(self):
2043 return fuzzy_and([self.args[0].is_integer, self.args[0].is_positive])
2046class divisor_sigma(Function):
2047 r"""
2048 Calculate the divisor function `\sigma_k(n)` for positive integer n
2050 ``divisor_sigma(n, k)`` is equal to ``sum([x**k for x in divisors(n)])``
2052 If n's prime factorization is:
2054 .. math ::
2055 n = \prod_{i=1}^\omega p_i^{m_i},
2057 then
2059 .. math ::
2060 \sigma_k(n) = \prod_{i=1}^\omega (1+p_i^k+p_i^{2k}+\cdots
2061 + p_i^{m_ik}).
2063 Parameters
2064 ==========
2066 n : integer
2068 k : integer, optional
2069 power of divisors in the sum
2071 for k = 0, 1:
2072 ``divisor_sigma(n, 0)`` is equal to ``divisor_count(n)``
2073 ``divisor_sigma(n, 1)`` is equal to ``sum(divisors(n))``
2075 Default for k is 1.
2077 Examples
2078 ========
2080 >>> from sympy.ntheory import divisor_sigma
2081 >>> divisor_sigma(18, 0)
2082 6
2083 >>> divisor_sigma(39, 1)
2084 56
2085 >>> divisor_sigma(12, 2)
2086 210
2087 >>> divisor_sigma(37)
2088 38
2090 See Also
2091 ========
2093 divisor_count, totient, divisors, factorint
2095 References
2096 ==========
2098 .. [1] https://en.wikipedia.org/wiki/Divisor_function
2100 """
2102 @classmethod
2103 def eval(cls, n, k=S.One):
2104 k = sympify(k)
2106 if n.is_prime:
2107 return 1 + n**k
2109 if n.is_Integer:
2110 if n <= 0:
2111 raise ValueError("n must be a positive integer")
2112 elif k.is_Integer:
2113 k = int(k)
2114 return Integer(math.prod(
2115 (p**(k*(e + 1)) - 1)//(p**k - 1) if k != 0
2116 else e + 1 for p, e in factorint(n).items()))
2117 else:
2118 return Mul(*[(p**(k*(e + 1)) - 1)/(p**k - 1) if k != 0
2119 else e + 1 for p, e in factorint(n).items()])
2121 if n.is_integer: # symbolic case
2122 args = []
2123 for p, e in (_.as_base_exp() for _ in Mul.make_args(n)):
2124 if p.is_prime and e.is_positive:
2125 args.append((p**(k*(e + 1)) - 1)/(p**k - 1) if
2126 k != 0 else e + 1)
2127 else:
2128 return
2129 return Mul(*args)
2132def core(n, t=2):
2133 r"""
2134 Calculate core(n, t) = `core_t(n)` of a positive integer n
2136 ``core_2(n)`` is equal to the squarefree part of n
2138 If n's prime factorization is:
2140 .. math ::
2141 n = \prod_{i=1}^\omega p_i^{m_i},
2143 then
2145 .. math ::
2146 core_t(n) = \prod_{i=1}^\omega p_i^{m_i \mod t}.
2148 Parameters
2149 ==========
2151 n : integer
2153 t : integer
2154 core(n, t) calculates the t-th power free part of n
2156 ``core(n, 2)`` is the squarefree part of ``n``
2157 ``core(n, 3)`` is the cubefree part of ``n``
2159 Default for t is 2.
2161 Examples
2162 ========
2164 >>> from sympy.ntheory.factor_ import core
2165 >>> core(24, 2)
2166 6
2167 >>> core(9424, 3)
2168 1178
2169 >>> core(379238)
2170 379238
2171 >>> core(15**11, 10)
2172 15
2174 See Also
2175 ========
2177 factorint, sympy.solvers.diophantine.diophantine.square_factor
2179 References
2180 ==========
2182 .. [1] https://en.wikipedia.org/wiki/Square-free_integer#Squarefree_core
2184 """
2186 n = as_int(n)
2187 t = as_int(t)
2188 if n <= 0:
2189 raise ValueError("n must be a positive integer")
2190 elif t <= 1:
2191 raise ValueError("t must be >= 2")
2192 else:
2193 y = 1
2194 for p, e in factorint(n).items():
2195 y *= p**(e % t)
2196 return y
2199class udivisor_sigma(Function):
2200 r"""
2201 Calculate the unitary divisor function `\sigma_k^*(n)` for positive integer n
2203 ``udivisor_sigma(n, k)`` is equal to ``sum([x**k for x in udivisors(n)])``
2205 If n's prime factorization is:
2207 .. math ::
2208 n = \prod_{i=1}^\omega p_i^{m_i},
2210 then
2212 .. math ::
2213 \sigma_k^*(n) = \prod_{i=1}^\omega (1+ p_i^{m_ik}).
2215 Parameters
2216 ==========
2218 k : power of divisors in the sum
2220 for k = 0, 1:
2221 ``udivisor_sigma(n, 0)`` is equal to ``udivisor_count(n)``
2222 ``udivisor_sigma(n, 1)`` is equal to ``sum(udivisors(n))``
2224 Default for k is 1.
2226 Examples
2227 ========
2229 >>> from sympy.ntheory.factor_ import udivisor_sigma
2230 >>> udivisor_sigma(18, 0)
2231 4
2232 >>> udivisor_sigma(74, 1)
2233 114
2234 >>> udivisor_sigma(36, 3)
2235 47450
2236 >>> udivisor_sigma(111)
2237 152
2239 See Also
2240 ========
2242 divisor_count, totient, divisors, udivisors, udivisor_count, divisor_sigma,
2243 factorint
2245 References
2246 ==========
2248 .. [1] https://mathworld.wolfram.com/UnitaryDivisorFunction.html
2250 """
2252 @classmethod
2253 def eval(cls, n, k=S.One):
2254 k = sympify(k)
2255 if n.is_prime:
2256 return 1 + n**k
2257 if n.is_Integer:
2258 if n <= 0:
2259 raise ValueError("n must be a positive integer")
2260 else:
2261 return Mul(*[1+p**(k*e) for p, e in factorint(n).items()])
2264class primenu(Function):
2265 r"""
2266 Calculate the number of distinct prime factors for a positive integer n.
2268 If n's prime factorization is:
2270 .. math ::
2271 n = \prod_{i=1}^k p_i^{m_i},
2273 then ``primenu(n)`` or `\nu(n)` is:
2275 .. math ::
2276 \nu(n) = k.
2278 Examples
2279 ========
2281 >>> from sympy.ntheory.factor_ import primenu
2282 >>> primenu(1)
2283 0
2284 >>> primenu(30)
2285 3
2287 See Also
2288 ========
2290 factorint
2292 References
2293 ==========
2295 .. [1] https://mathworld.wolfram.com/PrimeFactor.html
2297 """
2299 @classmethod
2300 def eval(cls, n):
2301 if n.is_Integer:
2302 if n <= 0:
2303 raise ValueError("n must be a positive integer")
2304 else:
2305 return len(factorint(n).keys())
2308class primeomega(Function):
2309 r"""
2310 Calculate the number of prime factors counting multiplicities for a
2311 positive integer n.
2313 If n's prime factorization is:
2315 .. math ::
2316 n = \prod_{i=1}^k p_i^{m_i},
2318 then ``primeomega(n)`` or `\Omega(n)` is:
2320 .. math ::
2321 \Omega(n) = \sum_{i=1}^k m_i.
2323 Examples
2324 ========
2326 >>> from sympy.ntheory.factor_ import primeomega
2327 >>> primeomega(1)
2328 0
2329 >>> primeomega(20)
2330 3
2332 See Also
2333 ========
2335 factorint
2337 References
2338 ==========
2340 .. [1] https://mathworld.wolfram.com/PrimeFactor.html
2342 """
2344 @classmethod
2345 def eval(cls, n):
2346 if n.is_Integer:
2347 if n <= 0:
2348 raise ValueError("n must be a positive integer")
2349 else:
2350 return sum(factorint(n).values())
2353def mersenne_prime_exponent(nth):
2354 """Returns the exponent ``i`` for the nth Mersenne prime (which
2355 has the form `2^i - 1`).
2357 Examples
2358 ========
2360 >>> from sympy.ntheory.factor_ import mersenne_prime_exponent
2361 >>> mersenne_prime_exponent(1)
2362 2
2363 >>> mersenne_prime_exponent(20)
2364 4423
2365 """
2366 n = as_int(nth)
2367 if n < 1:
2368 raise ValueError("nth must be a positive integer; mersenne_prime_exponent(1) == 2")
2369 if n > 51:
2370 raise ValueError("There are only 51 perfect numbers; nth must be less than or equal to 51")
2371 return MERSENNE_PRIME_EXPONENTS[n - 1]
2374def is_perfect(n):
2375 """Returns True if ``n`` is a perfect number, else False.
2377 A perfect number is equal to the sum of its positive, proper divisors.
2379 Examples
2380 ========
2382 >>> from sympy.ntheory.factor_ import is_perfect, divisors, divisor_sigma
2383 >>> is_perfect(20)
2384 False
2385 >>> is_perfect(6)
2386 True
2387 >>> 6 == divisor_sigma(6) - 6 == sum(divisors(6)[:-1])
2388 True
2390 References
2391 ==========
2393 .. [1] https://mathworld.wolfram.com/PerfectNumber.html
2394 .. [2] https://en.wikipedia.org/wiki/Perfect_number
2396 """
2398 n = as_int(n)
2399 if _isperfect(n):
2400 return True
2402 # all perfect numbers for Mersenne primes with exponents
2403 # less than or equal to 43112609 are known
2404 iknow = MERSENNE_PRIME_EXPONENTS.index(43112609)
2405 if iknow <= len(PERFECT) - 1 and n <= PERFECT[iknow]:
2406 # there may be gaps between this and larger known values
2407 # so only conclude in the range for which all values
2408 # are known
2409 return False
2410 if n%2 == 0:
2411 last2 = n % 100
2412 if last2 != 28 and last2 % 10 != 6:
2413 return False
2414 r, b = integer_nthroot(1 + 8*n, 2)
2415 if not b:
2416 return False
2417 m, x = divmod(1 + r, 4)
2418 if x:
2419 return False
2420 e, b = integer_log(m, 2)
2421 if not b:
2422 return False
2423 else:
2424 if n < 10**2000: # https://www.lirmm.fr/~ochem/opn/
2425 return False
2426 if n % 105 == 0: # not divis by 105
2427 return False
2428 if not any(n%m == r for m, r in [(12, 1), (468, 117), (324, 81)]):
2429 return False
2430 # there are many criteria that the factor structure of n
2431 # must meet; since we will have to factor it to test the
2432 # structure we will have the factors and can then check
2433 # to see whether it is a perfect number or not. So we
2434 # skip the structure checks and go straight to the final
2435 # test below.
2436 rv = divisor_sigma(n) - n
2437 if rv == n:
2438 if n%2 == 0:
2439 raise ValueError(filldedent('''
2440 This even number is perfect and is associated with a
2441 Mersenne Prime, 2^%s - 1. It should be
2442 added to SymPy.''' % (e + 1)))
2443 else:
2444 raise ValueError(filldedent('''In 1888, Sylvester stated: "
2445 ...a prolonged meditation on the subject has satisfied
2446 me that the existence of any one such [odd perfect number]
2447 -- its escape, so to say, from the complex web of conditions
2448 which hem it in on all sides -- would be little short of a
2449 miracle." I guess SymPy just found that miracle and it
2450 factors like this: %s''' % factorint(n)))
2453def is_mersenne_prime(n):
2454 """Returns True if ``n`` is a Mersenne prime, else False.
2456 A Mersenne prime is a prime number having the form `2^i - 1`.
2458 Examples
2459 ========
2461 >>> from sympy.ntheory.factor_ import is_mersenne_prime
2462 >>> is_mersenne_prime(6)
2463 False
2464 >>> is_mersenne_prime(127)
2465 True
2467 References
2468 ==========
2470 .. [1] https://mathworld.wolfram.com/MersennePrime.html
2472 """
2474 n = as_int(n)
2475 if _ismersenneprime(n):
2476 return True
2477 if not isprime(n):
2478 return False
2479 r, b = integer_log(n + 1, 2)
2480 if not b:
2481 return False
2482 raise ValueError(filldedent('''
2483 This Mersenne Prime, 2^%s - 1, should
2484 be added to SymPy's known values.''' % r))
2487def abundance(n):
2488 """Returns the difference between the sum of the positive
2489 proper divisors of a number and the number.
2491 Examples
2492 ========
2494 >>> from sympy.ntheory import abundance, is_perfect, is_abundant
2495 >>> abundance(6)
2496 0
2497 >>> is_perfect(6)
2498 True
2499 >>> abundance(10)
2500 -2
2501 >>> is_abundant(10)
2502 False
2503 """
2504 return divisor_sigma(n, 1) - 2 * n
2507def is_abundant(n):
2508 """Returns True if ``n`` is an abundant number, else False.
2510 A abundant number is smaller than the sum of its positive proper divisors.
2512 Examples
2513 ========
2515 >>> from sympy.ntheory.factor_ import is_abundant
2516 >>> is_abundant(20)
2517 True
2518 >>> is_abundant(15)
2519 False
2521 References
2522 ==========
2524 .. [1] https://mathworld.wolfram.com/AbundantNumber.html
2526 """
2527 n = as_int(n)
2528 if is_perfect(n):
2529 return False
2530 return n % 6 == 0 or bool(abundance(n) > 0)
2533def is_deficient(n):
2534 """Returns True if ``n`` is a deficient number, else False.
2536 A deficient number is greater than the sum of its positive proper divisors.
2538 Examples
2539 ========
2541 >>> from sympy.ntheory.factor_ import is_deficient
2542 >>> is_deficient(20)
2543 False
2544 >>> is_deficient(15)
2545 True
2547 References
2548 ==========
2550 .. [1] https://mathworld.wolfram.com/DeficientNumber.html
2552 """
2553 n = as_int(n)
2554 if is_perfect(n):
2555 return False
2556 return bool(abundance(n) < 0)
2559def is_amicable(m, n):
2560 """Returns True if the numbers `m` and `n` are "amicable", else False.
2562 Amicable numbers are two different numbers so related that the sum
2563 of the proper divisors of each is equal to that of the other.
2565 Examples
2566 ========
2568 >>> from sympy.ntheory.factor_ import is_amicable, divisor_sigma
2569 >>> is_amicable(220, 284)
2570 True
2571 >>> divisor_sigma(220) == divisor_sigma(284)
2572 True
2574 References
2575 ==========
2577 .. [1] https://en.wikipedia.org/wiki/Amicable_numbers
2579 """
2580 if m == n:
2581 return False
2582 a, b = (divisor_sigma(i) for i in (m, n))
2583 return a == b == (m + n)
2586def dra(n, b):
2587 """
2588 Returns the additive digital root of a natural number ``n`` in base ``b``
2589 which is a single digit value obtained by an iterative process of summing
2590 digits, on each iteration using the result from the previous iteration to
2591 compute a digit sum.
2593 Examples
2594 ========
2596 >>> from sympy.ntheory.factor_ import dra
2597 >>> dra(3110, 12)
2598 8
2600 References
2601 ==========
2603 .. [1] https://en.wikipedia.org/wiki/Digital_root
2605 """
2607 num = abs(as_int(n))
2608 b = as_int(b)
2609 if b <= 1:
2610 raise ValueError("Base should be an integer greater than 1")
2612 if num == 0:
2613 return 0
2615 return (1 + (num - 1) % (b - 1))
2618def drm(n, b):
2619 """
2620 Returns the multiplicative digital root of a natural number ``n`` in a given
2621 base ``b`` which is a single digit value obtained by an iterative process of
2622 multiplying digits, on each iteration using the result from the previous
2623 iteration to compute the digit multiplication.
2625 Examples
2626 ========
2628 >>> from sympy.ntheory.factor_ import drm
2629 >>> drm(9876, 10)
2630 0
2632 >>> drm(49, 10)
2633 8
2635 References
2636 ==========
2638 .. [1] https://mathworld.wolfram.com/MultiplicativeDigitalRoot.html
2640 """
2642 n = abs(as_int(n))
2643 b = as_int(b)
2644 if b <= 1:
2645 raise ValueError("Base should be an integer greater than 1")
2646 while n > b:
2647 mul = 1
2648 while n > 1:
2649 n, r = divmod(n, b)
2650 if r == 0:
2651 return 0
2652 mul *= r
2653 n = mul
2654 return n