Number Theory¶
Generating and counting primes.

class
diofant.ntheory.generate.
Sieve
[source]¶ An infinite list of prime numbers, implemented as a dynamically growing sieve of Eratosthenes. When a lookup is requested involving an odd number that has not been sieved, the sieve is automatically extended up to that number.
>>> from array import array # this line and next for doctest only >>> sieve._list = array('l', [2, 3, 5, 7, 11, 13])
>>> 25 in sieve False >>> sieve._list array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23])

extend
(n)[source]¶ Grow the sieve to cover all primes <= n (a real number).
Examples
>>> from array import array # this line and next for doctest only >>> sieve._list = array('l', [2, 3, 5, 7, 11, 13])
>>> sieve.extend(30) >>> sieve[10] == 29 True

extend_to_no
(i)[source]¶ Extend to include the ith prime number.
i must be an integer.
The list is extended by 50% if it is too short, so it is likely that it will be longer than requested.
Examples
>>> from array import array # this line and next for doctest only >>> sieve._list = array('l', [2, 3, 5, 7, 11, 13])
>>> sieve.extend_to_no(9) >>> sieve._list array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23])


diofant.ntheory.generate.
cycle_length
(f, x0, nmax=None, values=False)[source]¶ For a given iterated sequence, return a generator that gives the length of the iterated cycle (lambda) and the length of terms before the cycle begins (mu); if
values
is True then the terms of the sequence will be returned instead. The sequence is started with valuex0
.Note: more than the first lambda + mu terms may be returned and this is the cost of cycle detection with Brent’s method; there are, however, generally less terms calculated than would have been calculated if the proper ending point were determined, e.g. by using Floyd’s method.
This will yield successive values of i <– func(i):
>>> def iter(func, i): ... while 1: ... ii = func(i) ... yield ii ... i = ii ...
A function is defined:
>>> func = lambda i: (i**2 + 1) % 51
and given a seed of 4 and the mu and lambda terms calculated:
>>> next(cycle_length(func, 4)) (6, 2)
We can see what is meant by looking at the output:
>>> n = cycle_length(func, 4, values=True) >>> list(ni for ni in n) [17, 35, 2, 5, 26, 14, 44, 50, 2, 5, 26, 14]
There are 6 repeating values after the first 2.
If a sequence is suspected of being longer than you might wish,
nmax
can be used to exit early (and mu will be returned as None):>>> next(cycle_length(func, 4, nmax = 4)) (4, None) >>> [ni for ni in cycle_length(func, 4, nmax = 4, values=True)] [17, 35, 2, 5]
References

diofant.ntheory.generate.
nextprime
(n, ith=1)[source]¶ Return the ith prime greater than n.
i must be an integer.
Notes
Potential primes are located at 6*j +/ 1. This property is used during searching.
>>> [(i, nextprime(i)) for i in range(10, 15)] [(10, 11), (11, 13), (12, 13), (13, 17), (14, 17)] >>> nextprime(2, ith=2) # the 2nd prime after 2 5
See also
prevprime()
Return the largest prime smaller than n
primerange()
Generate all primes in a given range

diofant.ntheory.generate.
prevprime
(n)[source]¶ Return the largest prime smaller than n.
Notes
Potential primes are located at 6*j +/ 1. This property is used during searching.
>>> [(i, prevprime(i)) for i in range(10, 15)] [(10, 7), (11, 7), (12, 11), (13, 11), (14, 13)]
See also
nextprime()
Return the ith prime greater than n
primerange()
Generates all primes in a given range

diofant.ntheory.generate.
prime
(nth)[source]¶ Return the nth prime, with the primes indexed as prime(1) = 2, prime(2) = 3, etc…. The nth prime is approximately n*log(n) and can never be larger than 2**n.
References
Examples
>>> prime(10) 29 >>> prime(1) 2
See also
diofant.ntheory.primetest.isprime()
Test if n is prime
primerange()
Generate all primes in a given range
primepi()
Return the number of primes less than or equal to n

diofant.ntheory.generate.
primepi
(n)[source]¶ Return the value of the prime counting function pi(n) = the number of prime numbers less than or equal to n.
Examples
>>> primepi(25) 9
See also
diofant.ntheory.primetest.isprime()
Test if n is prime
primerange()
Generate all primes in a given range
prime()
Return the nth prime

diofant.ntheory.generate.
primerange
(a, b)[source]¶ Generate a list of all prime numbers in the range [a, b).
If the range exists in the default sieve, the values will be returned from there; otherwise values will be returned but will not modify the sieve.
Notes
Some famous conjectures about the occurrence of primes in a given range are:
 Twin primes: though often not, the following will give 2 primes
an infinite number of times: primerange(6*n  1, 6*n + 2)
 Legendre’s: the following always yields at least one prime
primerange(n**2, (n+1)**2+1)
 Bertrand’s (proven): there is always a prime in the range
primerange(n, 2*n)
 Brocard’s: there are at least four primes in the range
primerange(prime(n)**2, prime(n+1)**2)
The average gap between primes is log(n); the gap between primes can be arbitrarily large since sequences of composite numbers are arbitrarily large, e.g. the numbers in the sequence n! + 2, n! + 3 … n! + n are all composite.
References
Examples
>>> print([i for i in primerange(1, 30)]) [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
The Sieve method, primerange, is generally faster but it will occupy more memory as the sieve stores values. The default instance of Sieve, named sieve, can be used:
>>> list(sieve.primerange(1, 30)) [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
See also
nextprime()
Return the ith prime greater than n
prevprime()
Return the largest prime smaller than n
randprime()
Returns a random prime in a given range
primorial()
Returns the product of primes based on condition
Sieve.primerange()
return range from already computed primes or extend the sieve to contain the requested range.

diofant.ntheory.generate.
primorial
(n, nth=True)[source]¶ Returns the product of the first n primes (default) or the primes less than or equal to n (when
nth=False
).>>> primorial(4) # the first 4 primes are 2, 3, 5, 7 210 >>> primorial(4, nth=False) # primes <= 4 are 2 and 3 6 >>> primorial(1) 2 >>> primorial(1, nth=False) 1 >>> primorial(sqrt(101), nth=False) 210
One can argue that the primes are infinite since if you take a set of primes and multiply them together (e.g. the primorial) and then add or subtract 1, the result cannot be divided by any of the original factors, hence either 1 or more new primes must divide this product of primes.
In this case, the number itself is a new prime:
>>> factorint(primorial(4) + 1) {211: 1}
In this case two new primes are the factors:
>>> factorint(primorial(4)  1) {11: 1, 19: 1}
Here, some primes smaller and larger than the primes multiplied together are obtained:
>>> p = list(primerange(10, 20)) >>> sorted(set(primefactors(Mul(*p) + 1)).difference(set(p))) [2, 5, 31, 149]
See also
primerange()
Generate all primes in a given range

diofant.ntheory.generate.
randprime
(a, b)[source]¶ Return a random prime number in the range [a, b).
Bertrand’s postulate assures that randprime(a, 2*a) will always succeed for a > 1.
References
Examples
>>> randprime(1, 30) #doctest: +SKIP 13 >>> isprime(randprime(1, 30)) True
See also
primerange()
Generate all primes in a given range
Integer factorization

diofant.ntheory.factor_.
antidivisor_count
(n)[source]¶ Return the number of antidivisors of
n
.References
formula from https://oeis.org/A066272
Examples
>>> antidivisor_count(13) 4 >>> antidivisor_count(27) 5
See also
factorint()
,divisors()
,antidivisors()
,divisor_count()
,totient()

diofant.ntheory.factor_.
antidivisors
(n, generator=False)[source]¶ Return all antidivisors of n sorted from 1..n by default.
Antidivisors of n are numbers that do not divide n by the largest possible margin. If generator is True an unordered generator is returned.
References
definition is described in http://oeis.org/A066272/a066272a.html
Examples
>>> antidivisors(24) [7, 16]
>>> sorted(antidivisors(128, generator=True)) [3, 5, 15, 17, 51, 85]

diofant.ntheory.factor_.
core
(n, t=2)[source]¶ Calculate core(n,t) = \(core_t(n)\) of a positive integer n.
core_2(n)
is equal to the squarefree part of nIf n’s prime factorization is:
\[n = \prod_{i=1}^\omega p_i^{m_i},\]then
\[core_t(n) = \prod_{i=1}^\omega p_i^{m_i \mod t}.\] Parameters
t (core(n,t) calculates the tth power free part of n) –
core(n, 2)
is the squarefree part ofn
core(n, 3)
is the cubefree part ofn
Default for t is 2.
References
Examples
>>> from diofant.ntheory.factor_ import core >>> core(24, 2) 6 >>> core(9424, 3) 1178 >>> core(379238) 379238 >>> core(15**11, 10) 15

diofant.ntheory.factor_.
divisor_count
(n, modulus=1)[source]¶ Return the number of divisors of
n
.If
modulus
is not 1 then only those that are divisible bymodulus
are counted.References
Examples
>>> divisor_count(6) 4
See also

class
diofant.ntheory.factor_.
divisor_sigma
[source]¶ Calculate the divisor function \(\sigma_k(n)\) for positive integer n.
divisor_sigma(n, k)
is equal tosum(x**k for x in divisors(n))
If n’s prime factorization is:
\[n = \prod_{i=1}^\omega p_i^{m_i},\]then
\[\sigma_k(n) = \prod_{i=1}^\omega (1+p_i^k+p_i^{2k}+\cdots + p_i^{m_ik}).\] Parameters
k (power of divisors in the sum) – for k = 0, 1:
divisor_sigma(n, 0)
is equal todivisor_count(n)
divisor_sigma(n, 1)
is equal tosum(divisors(n))
Default for k is 1.
References
Examples
>>> divisor_sigma(18, 0) 6 >>> divisor_sigma(39, 1) 56 >>> divisor_sigma(12, 2) 210 >>> divisor_sigma(37) 38
See also

diofant.ntheory.factor_.
divisors
(n, generator=False)[source]¶ Return all divisors of n.
Divisors are sorted from 1..n by default. If generator is True an unordered generator is returned.
The number of divisors of n can be quite large if there are many prime factors (counting repeated factors). If only the number of factors is desired use divisor_count(n).
Examples
>>> divisors(24) [1, 2, 3, 4, 6, 8, 12, 24] >>> divisor_count(24) 8
>>> list(divisors(120, generator=True)) [1, 2, 4, 8, 3, 6, 12, 24, 5, 10, 20, 40, 15, 30, 60, 120]
See also
References

diofant.ntheory.factor_.
factorint
(n, limit=None, use_trial=True, use_rho=True, use_pm1=True, verbose=False, visual=None)[source]¶ Given a positive integer
n
,factorint(n)
returns a dict containing the prime factors ofn
as keys and their respective multiplicities as values. For example:>>> factorint(2000) # 2000 = (2**4) * (5**3) {2: 4, 5: 3} >>> factorint(65537) # This number is prime {65537: 1}
For input less than 2, factorint behaves as follows:
factorint(1)
returns the empty factorization,{}
factorint(0)
returns{0:1}
factorint(n)
adds1:1
to the factors and then factorsn
Partial Factorization:
If
limit
(> 3) is specified, the search is stopped after performing trial division up to (and including) the limit (or taking a corresponding number of rho/p1 steps). This is useful if one has a large number and only is interested in finding small factors (if any). Note that setting a limit does not prevent larger factors from being found early; it simply means that the largest factor may be composite. Since checking for perfect power is relatively cheap, it is done regardless of the limit setting.This number, for example, has two small factors and a huge semiprime factor that cannot be reduced easily:
>>> a = 1407633717262338957430697921446883 >>> f = factorint(a, limit=10000) >>> f {7: 1, 991: 1, 202916782076162456022877024859: 1} >>> isprime(max(f)) False
This number has a small factor and a residual perfect power whose base is greater than the limit:
>>> factorint(3*101**7, limit=5) {3: 1, 101: 7}
Visual Factorization:
If
visual
is set toTrue
, then it will return a visual factorization of the integer. For example:>>> pprint(factorint(4200, visual=True), use_unicode=False) 3 1 2 1 2 *3 *5 *7
Note that this is achieved by using the evaluate=False flag in Mul and Pow. If you do other manipulations with an expression where evaluate=False, it may evaluate. Therefore, you should use the visual option only for visualization, and use the normal dictionary returned by visual=False if you want to perform operations on the factors.
You can easily switch between the two forms by sending them back to factorint:
>>> regular = factorint(1764); regular {2: 2, 3: 2, 7: 2} >>> pprint(factorint(regular), use_unicode=False) 2 2 2 2 *3 *7
>>> visual = factorint(1764, visual=True); pprint(visual, use_unicode=False) 2 2 2 2 *3 *7 >>> print(factorint(visual)) {2: 2, 3: 2, 7: 2}
If you want to send a number to be factored in a partially factored form you can do so with a dictionary or unevaluated expression:
>>> factorint(factorint({4: 2, 12: 3})) # twice to toggle to dict form {2: 10, 3: 3} >>> factorint(Mul(4, 12, evaluate=False)) {2: 4, 3: 1}
The table of the output logic is:
Input
True
False
other
dict
mul
dict
mul
n
mul
dict
dict
mul
mul
dict
dict
Notes
The function switches between multiple algorithms. Trial division quickly finds small factors (of the order 15 digits), and finds all large factors if given enough time. The Pollard rho and p1 algorithms are used to find large factors ahead of time; they will often find factors of the order of 10 digits within a few seconds:
>>> factors = factorint(12345678910111213141516) >>> for base, exp in sorted(factors.items()): ... print('%s %s' % (base, exp)) ... 2 2 2507191691 1 1231026625769 1
Any of these methods can optionally be disabled with the following boolean parameters:
use_trial
: Toggle use of trial divisionuse_rho
: Toggle use of Pollard’s rho methoduse_pm1
: Toggle use of Pollard’s p1 method
factorint
also periodically checks if the remaining part is a prime number or a perfect power, and in those cases stops.If
verbose
is set toTrue
, detailed progress is printed.See also

diofant.ntheory.factor_.
factorrat
(rat, limit=None, use_trial=True, use_rho=True, use_pm1=True, verbose=False, visual=None)[source]¶ Given a Rational
r
,factorrat(r)
returns a dict containing the prime factors ofr
as keys and their respective multiplicities as values. For example:>>> factorrat(Rational(8, 9)) # = (2**3) * (3**2) {2: 3, 3: 2} >>> factorrat(Rational(1, 987)) # = 1 * (3**1) * (7**1) * (47**1) {1: 1, 3: 1, 7: 1, 47: 1}
Please see the docstring for
factorint
for detailed explanations and examples of the following keywords:limit
: Integer limit up to which trial division is doneuse_trial
: Toggle use of trial divisionuse_rho
: Toggle use of Pollard’s rho methoduse_pm1
: Toggle use of Pollard’s p1 methodverbose
: Toggle detailed printing of progressvisual
: Toggle product form of output

diofant.ntheory.factor_.
multiplicity
(p, n)[source]¶ Find the greatest integer m such that p**m divides n.
Examples
>>> [multiplicity(5, n) for n in [8, 5, 25, 125, 250]] [0, 1, 2, 3, 3] >>> multiplicity(3, Rational(1, 9)) 2

diofant.ntheory.factor_.
perfect_power
(n, candidates=None, big=True, factor=True)[source]¶ Return
(b, e)
such thatn
==b**e
ifn
is a perfect power; otherwise returnFalse
.By default, the base is recursively decomposed and the exponents collected so the largest possible
e
is sought. Ifbig=False
then the smallest possiblee
(thus prime) will be chosen.If
candidates
for exponents are given, they are assumed to be sorted and the first one that is larger than the computed maximum will signal failure for the routine.If
factor=True
then simultaneous factorization of n is attempted since finding a factor indicates the only possible root for n. This is True by default since only a few small factors will be tested in the course of searching for the perfect power.Examples
>>> perfect_power(16) (2, 4) >>> perfect_power(16, big = False) (4, 2)

diofant.ntheory.factor_.
pollard_pm1
(n, B=10, a=2, retries=0, seed=1234)[source]¶ Use Pollard’s p1 method to try to extract a nontrivial factor of
n
. Either a divisor (perhaps composite) orNone
is returned.The value of
a
is the base that is used in the test gcd(a**M  1, n). The default is 2. Ifretries
> 0 then if no factor is found after the first attempt, a newa
will be generated randomly (using theseed
) and the process repeated.Note: the value of M is lcm(1..B) = reduce(ilcm, range(2, B + 1)).
A search is made for factors next to even numbers having a power smoothness less than
B
. Choosing a larger B increases the likelihood of finding a larger factor but takes longer. Whether a factor of n is found or not depends ona
and the power smoothness of the even mumber just less than the factor p (hence the name p  1).Although some discussion of what constitutes a good
a
some descriptions are hard to interpret. At the modular.math site referenced below it is stated that if gcd(a**M  1, n) = N then a**M % q**r is 1 for every prime power divisor of N. But consider the following:>>> n = 257*1009 >>> smoothness_p(n) (1, [(257, (1, 2, 256)), (1009, (1, 7, 16))])
So we should (and can) find a root with B=16:
>>> pollard_pm1(n, B=16, a=3) 1009
If we attempt to increase B to 256 we find that it doesn’t work:
>>> pollard_pm1(n, B=256) >>>
But if the value of
a
is changed we find that only multiples of 257 work, e.g.:>>> pollard_pm1(n, B=256, a=257) 1009
Checking different
a
values shows that all the ones that didn’t work had a gcd value not equal ton
but equal to one of the factors:>>> M = 1 >>> for i in range(2, 256): ... M = ilcm(M, i) ... >>> {math.gcd(pow(a, M, n)  1, n) for a in range(2, 256) if ... math.gcd(pow(a, M, n)  1, n) != n} {1009}
But does aM % d for every divisor of n give 1?
>>> aM = pow(255, M, n) >>> [(d, aM%Pow(*d.args)) for d in factorint(n, visual=True).args] [(257**1, 1), (1009**1, 1)]
No, only one of them. So perhaps the principle is that a root will be found for a given value of B provided that:
the power smoothness of the p  1 value next to the root does not exceed B
a**M % p != 1 for any of the divisors of n.
By trying more than one
a
it is possible that one of them will yield a factor.Examples
With the default smoothness bound, this number can’t be cracked:
>>> pollard_pm1(21477639576571)
Increasing the smoothness bound helps:
>>> pollard_pm1(21477639576571, B=2000) 4410317
Looking at the smoothness of the factors of this number we find:
>>> print(smoothness_p(21477639576571, visual=1)) p**i=4410317**1 has p1 B=1787, Bpow=1787 p**i=4869863**1 has p1 B=2434931, Bpow=2434931
The B and Bpow are the same for the p  1 factorizations of the divisors because those factorizations had a very large prime factor:
>>> factorint(4410317  1) {2: 2, 617: 1, 1787: 1} >>> factorint(48698631) {2: 1, 2434931: 1}
Note that until B reaches the Bpow value of 1787, the number is not cracked;
>>> pollard_pm1(21477639576571, B=1786) >>> pollard_pm1(21477639576571, B=1787) 4410317
The B value has to do with the factors of the number next to the divisor, not the divisors themselves. A worst case scenario is that the number next to the factor p has a large prime divisisor or is a perfect power. If these conditions apply then the powersmoothness will be about p/2 or p. The more realistic is that there will be a large prime factor next to p requiring a B value on the order of p/2. Although primes may have been searched for up to this level, the p/2 is a factor of p  1, something that we don’t know. The modular.math reference below states that 15% of numbers in the range of 10**15 to 15**15 + 10**4 are 10**6 power smooth so a B of 10**6 will fail 85% of the time in that range. From 10**8 to 10**8 + 10**3 the percentages are nearly reversed…but in that range the simple trial division is quite fast.
References
Richard Crandall & Carl Pomerance (2005), “Prime Numbers: A Computational Perspective”, Springer, 2nd edition, 236238
https://web.archive.org/web/20170830055619/http://www.cs.toronto.edu/~yuvalf/Factorization.pdf

diofant.ntheory.factor_.
pollard_rho
(n, s=2, a=1, retries=5, seed=1234, max_steps=None, F=None)[source]¶ Use Pollard’s rho method to try to extract a nontrivial factor of
n
. The returned factor may be a composite number. If no factor is found,None
is returned.The algorithm generates pseudorandom values of x with a generator function, replacing x with F(x). If F is not supplied then the function x**2 +
a
is used. The first value supplied to F(x) iss
. Upon failure (ifretries
is > 0) a newa
ands
will be supplied; thea
will be ignored if F was supplied.The sequence of numbers generated by such functions generally have a a leadup to some number and then loop around back to that number and begin to repeat the sequence, e.g. 1, 2, 3, 4, 5, 3, 4, 5 – this leader and loop look a bit like the Greek letter rho, and thus the name, ‘rho’.
For a given function, very different leaderloop values can be obtained so it is a good idea to allow for retries:
>>> n = 16843009 >>> F = lambda x: (2048*pow(x, 2, n) + 32767)%n >>> for s in range(5): ... print('loop length = %4i; leader length = %3i' % next(cycle_length(F, s))) ... loop length = 2489; leader length = 42 loop length = 78; leader length = 120 loop length = 1482; leader length = 99 loop length = 1482; leader length = 285 loop length = 1482; leader length = 100
Here is an explicit example where there is a two element leadup to a sequence of 3 numbers (11, 14, 4) that then repeat:
>>> x = 2 >>> for i in range(9): ... x = (x**2 + 12)%17 ... print(x) ... 16 13 11 14 4 11 14 4 11 >>> next(cycle_length(lambda x: (x**2+12)%17, 2)) (3, 2) >>> list(cycle_length(lambda x: (x**2+12)%17, 2, values=True)) [16, 13, 11, 14, 4]
Instead of checking the differences of all generated values for a gcd with n, only the kth and 2*kth numbers are checked, e.g. 1st and 2nd, 2nd and 4th, 3rd and 6th until it has been detected that the loop has been traversed. Loops may be many thousands of steps long before rho finds a factor or reports failure. If
max_steps
is specified, the iteration is cancelled with a failure after the specified number of steps.Examples
>>> n = 16843009 >>> F = lambda x: (2048*pow(x, 2, n) + 32767) % n >>> pollard_rho(n, F=F) 257
Use the default setting with a bad value of
a
and no retries:>>> pollard_rho(n, a=n2, retries=0)
If retries is > 0 then perhaps the problem will correct itself when new values are generated for a:
>>> pollard_rho(n, a=n2, retries=1) 257
References
Richard Crandall & Carl Pomerance (2005), “Prime Numbers: A Computational Perspective”, Springer, 2nd edition, 229231

diofant.ntheory.factor_.
primefactors
(n, limit=None, verbose=False)[source]¶ Return a sorted list of n’s prime factors, ignoring multiplicity and any composite factor that remains if the limit was set too low for complete factorization. Unlike factorint(), primefactors() does not return 1 or 0.
Examples
>>> primefactors(6) [2, 3] >>> primefactors(5) [5]
>>> sorted(factorint(123456).items()) [(2, 6), (3, 1), (643, 1)] >>> primefactors(123456) [2, 3, 643]
>>> sorted(factorint(10000000001, limit=200).items()) [(101, 1), (99009901, 1)] >>> isprime(99009901) False >>> primefactors(10000000001, limit=300) [101]
See also

diofant.ntheory.factor_.
smoothness
(n)[source]¶ Return the Bsmooth and Bpower smooth values of n.
The smoothness of n is the largest prime factor of n; the power smoothness is the largest divisor raised to its multiplicity.
>>> smoothness(2**7*3**2) (3, 128) >>> smoothness(2**4*13) (13, 16) >>> smoothness(2) (2, 2)
See also

diofant.ntheory.factor_.
smoothness_p
(n, m=1, power=0, visual=None)[source]¶ Return a list of [m, (p, (M, sm(p + m), psm(p + m)))…] where:
p**M is the basep divisor of n
sm(p + m) is the smoothness of p + m (m = 1 by default)
psm(p + m) is the power smoothness of p + m
The list is sorted according to smoothness (default) or by power smoothness if power=1.
The smoothness of the numbers to the left (m = 1) or right (m = 1) of a factor govern the results that are obtained from the p +/ 1 type factoring methods.
>>> smoothness_p(10431, m=1) (1, [(3, (2, 2, 4)), (19, (1, 5, 5)), (61, (1, 31, 31))]) >>> smoothness_p(10431) (1, [(3, (2, 2, 2)), (19, (1, 3, 9)), (61, (1, 5, 5))]) >>> smoothness_p(10431, power=1) (1, [(3, (2, 2, 2)), (61, (1, 5, 5)), (19, (1, 3, 9))])
If visual=True then an annotated string will be returned:
>>> print(smoothness_p(21477639576571, visual=1)) p**i=4410317**1 has p1 B=1787, Bpow=1787 p**i=4869863**1 has p1 B=2434931, Bpow=2434931
This string can also be generated directly from a factorization dictionary and vice versa:
>>> f = factorint(17*9) >>> f {3: 2, 17: 1} >>> smoothness_p(f) 'p**i=3**2 has p1 B=2, Bpow=2\np**i=17**1 has p1 B=2, Bpow=16' >>> smoothness_p(_) {3: 2, 17: 1}
The table of the output logic is:
Visual
Input
True
False
other
dict
str
tuple
str
str
str
tuple
dict
tuple
str
tuple
str
n
str
tuple
tuple
mul
str
tuple
tuple
See also

class
diofant.ntheory.factor_.
totient
[source]¶ Calculate the Euler totient function phi(n)
>>> totient(1) 1 >>> totient(25) 20
See also

diofant.ntheory.factor_.
trailing
(n)[source]¶ Count the number of trailing zero digits in the binary representation of n, i.e. determine the largest power of 2 that divides n.
Examples
>>> trailing(128) 7 >>> trailing(63) 0

diofant.ntheory.modular.
crt
(M, U, symmetric=False, check=True)[source]¶ Chinese Remainder Theorem.
The moduli in M are assumed to be pairwise coprime. The output is then an integer f, such that f = u_i mod m_i for each pair out of U and M. If
symmetric
is False a positive integer will be returned, else f will be less than or equal to the LCM of the moduli, and thus f may be negative.If the moduli are not coprime the correct result will be returned if/when the test of the result is found to be incorrect. This result will be None if there is no solution.
The keyword
check
can be set to False if it is known that the moduli are coprime.Examples
>>> crt([99, 97, 95], [49, 76, 65]) (639985, 912285)
This is the correct result because:
>>> [639985 % m for m in [99, 97, 95]] [49, 76, 65]
If the moduli are not coprime, you may receive an incorrect result if you use
check=False
:>>> crt([12, 6, 17], [3, 4, 2], check=False) (954, 1224) >>> [954 % m for m in [12, 6, 17]] [6, 0, 2] >>> crt([12, 6, 17], [3, 4, 2]) is None True >>> crt([3, 6], [2, 5]) (5, 6)
Notes
Rather than checking that all pairs of moduli share no GCD (an O(n**2) test) and rather than factoring all moduli and seeing that there is no factor in common, a check that the result gives the indicated residuals is performed – an O(n) operation.
See also

diofant.ntheory.modular.
crt1
(M)[source]¶ First part of Chinese Remainder Theorem, for multiple application.
Examples
>>> crt1([18, 42, 6]) (4536, [252, 108, 756], [0, 2, 0])

diofant.ntheory.modular.
crt2
(M, U, p, E, S, symmetric=False)[source]¶ Second part of Chinese Remainder Theorem, for multiple application.
Examples
>>> mm, e, s = crt1([18, 42, 6]) >>> crt2([18, 42, 6], [0, 0, 0], mm, e, s) (0, 4536)

diofant.ntheory.modular.
integer_rational_reconstruction
(c, m, domain)[source]¶ Reconstruct a rational number \(\frac a b\) from
\[c = \frac a b \; \mathrm{mod} \, m,\]where \(c\) and \(m\) are integers.
The algorithm is based on the Euclidean Algorithm. In general, \(m\) is not a prime number, so it is possible that \(b\) is not invertible modulo \(m\). In that case
None
is returned. Parameters
c (Integer) – \(c = \frac a b \; \mathrm{mod} \, m\)
m (Integer) – modulus, not necessarily prime
domain (IntegerRing) – \(a, b, c\) are elements of
domain
 Returns
frac (Rational) – either \(\frac a b\) in \(\mathbb Q\) or
None
References

diofant.ntheory.modular.
solve_congruence
(*remainder_modulus_pairs, **hint)[source]¶ Compute the integer
n
that has the residualai
when it is divided bymi
where theai
andmi
are given as pairs to this function: ((a1, m1), (a2, m2), …). If there is no solution, return. Otherwise returnn
and its modulus.The
mi
values need not be coprime. If it is known that the moduli are not coprime then the hintcheck
can be set to False (default=True) and the check for a quicker solution via crt() (valid when the moduli are coprime) will be skipped.If the hint
symmetric
is True (default is False), the value ofn
will be within 1/2 of the modulus, possibly negative.Examples
What number is 2 mod 3, 3 mod 5 and 2 mod 7?
>>> solve_congruence((2, 3), (3, 5), (2, 7)) (23, 105) >>> [23 % m for m in [3, 5, 7]] [2, 3, 2]
If you prefer to work with all remainder in one list and all moduli in another, send the arguments like this:
>>> solve_congruence(*zip((2, 3, 2), (3, 5, 7))) (23, 105)
The moduli need not be coprime; in this case there may or may not be a solution:
>>> solve_congruence((2, 3), (4, 6)) is None True
>>> solve_congruence((2, 3), (5, 6)) (5, 6)
The symmetric flag will make the result be within 1/2 of the modulus:
>>> solve_congruence((2, 3), (5, 6), symmetric=True) (1, 6)
See also
crt()
high level routine implementing the Chinese Remainder Theorem

diofant.ntheory.modular.
symmetric_residue
(a, m)[source]¶ Return the residual mod m such that it is within half of the modulus.
>>> symmetric_residue(1, 6) 1 >>> symmetric_residue(4, 6) 2

diofant.ntheory.multinomial.
binomial_coefficients
(n)[source]¶ Return a dictionary containing pairs \({(k1,k2) : C_{kn}}\) where \(C_{kn}\) are binomial coefficients and \(n=k1+k2\).
Examples
>>> binomial_coefficients(9) {(0, 9): 1, (1, 8): 9, (2, 7): 36, (3, 6): 84, (4, 5): 126, (5, 4): 126, (6, 3): 84, (7, 2): 36, (8, 1): 9, (9, 0): 1}

diofant.ntheory.multinomial.
binomial_coefficients_list
(n)[source]¶ Return a list of binomial coefficients as rows of the Pascal’s triangle.
Examples
>>> binomial_coefficients_list(9) [1, 9, 36, 84, 126, 126, 84, 36, 9, 1]

diofant.ntheory.multinomial.
multinomial_coefficients
(m, n)[source]¶ Return a dictionary containing pairs
{(k1,k2,..,km) : C_kn}
whereC_kn
are multinomial coefficients such thatn=k1+k2+..+km
.Examples
>>> multinomial_coefficients(2, 5) {(0, 5): 1, (1, 4): 5, (2, 3): 10, (3, 2): 10, (4, 1): 5, (5, 0): 1}
Notes
The algorithm is based on the following result:
\[\binom{n}{k_1, \ldots, k_m} = \frac{k_1 + 1}{n  k_1} \sum_{i=2}^m \binom{n}{k_1 + 1, \ldots, k_i  1, \ldots}\]

diofant.ntheory.multinomial.
multinomial_coefficients_iterator
(m, n, _tuple=<class 'tuple'>)[source]¶ multinomial coefficient iterator
This routine has been optimized for \(m\) large with respect to \(n\) by taking advantage of the fact that when the monomial tuples \(t\) are stripped of zeros, their coefficient is the same as that of the monomial tuples from
multinomial_coefficients(n, n)
. Therefore, the latter coefficients are precomputed to save memory and time.>>> m53, m33 = multinomial_coefficients(5, 3), multinomial_coefficients(3, 3) >>> m53[(0, 0, 0, 1, 2)] == m53[(0, 0, 1, 0, 2)] == m53[(1, 0, 2, 0, 0)] == m33[(0, 1, 2)] True
Examples
>>> it = multinomial_coefficients_iterator(20, 3) >>> next(it) ((3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 1)

diofant.ntheory.partitions_.
npartitions
(n)[source]¶ Calculate the partition function P(n), i.e. the number of ways that n can be written as a sum of positive integers.
P(n) is computed using the HardyRamanujanRademacher formula.
The correctness of this implementation has been tested for 10**n up to n = 8.
Examples
>>> npartitions(25) 1958
References
Primality testing

diofant.ntheory.primetest.
is_square
(n, prep=True)[source]¶ Return True if n == a * a for some integer a, else False.
If n is suspected of not being a square then this is a quick method of confirming that it is not.
References
See also

diofant.ntheory.primetest.
isprime
(n)[source]¶ Test if n is a prime number (True) or not (False). For n < 10**16 the answer is accurate; greater n values have a small probability of actually being pseudoprimes.
Negative primes (e.g. 2) are not considered prime.
The function first looks for trivial factors, and if none is found, performs a safe MillerRabin strong pseudoprime test with bases that are known to prove a number prime. Finally, a general MillerRabin test is done with the first k bases which will report a pseudoprime as a prime with an error of about 4**k. The current value of k is 46 so the error is about 2 x 10**28.
Examples
>>> isprime(13) True >>> isprime(15) False
See also
diofant.ntheory.generate.primerange()
Generates all primes in a given range
diofant.ntheory.generate.primepi()
Return the number of primes less than or equal to n
diofant.ntheory.generate.prime()
Return the nth prime

diofant.ntheory.primetest.
mr
(n, bases)[source]¶ Perform a MillerRabin strong pseudoprime test on n using a given list of bases/witnesses.
References
Richard Crandall & Carl Pomerance (2005), “Prime Numbers: A Computational Perspective”, Springer, 2nd edition, 135138
A list of thresholds and the bases they require are here: https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Deterministic_variants
Examples
>>> mr(1373651, [2, 3]) False >>> mr(479001599, [31, 73]) True

diofant.ntheory.residue_ntheory.
discrete_log
(n, a, b, order=None, prime_order=None)[source]¶ Compute the discrete logarithm of
a
to the baseb
modulon
.This is a recursive function to reduce the discrete logarithm problem in cyclic groups of composite order to the problem in cyclic groups of prime order.
Notes
It employs different algorithms depending on the problem (subgroup order size, prime order or not):
Trial multiplication
Babystep giantstep
PohligHellman
References
Examples
>>> discrete_log(41, 15, 7) 3

diofant.ntheory.residue_ntheory.
is_nthpow_residue
(a, n, m)[source]¶ Returns True if
x**n == a (mod m)
has solutions.References
Hackman “Elementary Number Theory” (2009), page 76

diofant.ntheory.residue_ntheory.
is_primitive_root
(a, p)[source]¶ Returns True if
a
is a primitive root ofp
a
is said to be the primitive root ofp
if gcd(a, p) == 1 and totient(p) is the smallest positive number s.t.:a**totient(p) cong 1 mod(p)
Examples
>>> is_primitive_root(3, 10) True >>> is_primitive_root(9, 10) False >>> n_order(3, 10) == totient(10) True >>> n_order(9, 10) == totient(10) False

diofant.ntheory.residue_ntheory.
is_quad_residue
(a, p)[source]¶ Returns True if
a
(modp
) is in the set of squares modp
, i.e a % p in {i**2 % p for i in range(p)}. Ifp
is an odd prime, an iterative method is used to make the determination:>>> sorted({i**2 % 7 for i in range(7)}) [0, 1, 2, 4] >>> [j for j in range(7) if is_quad_residue(j, 7)] [0, 1, 2, 4]
See also

diofant.ntheory.residue_ntheory.
jacobi_symbol
(m, n)[source]¶ Returns the Jacobi symbol \((m / n)\).
For any integer
m
and any positive odd integern
the Jacobi symbol is defined as the product of the Legendre symbols corresponding to the prime factors ofn
:\[\genfrac(){}{}{m}{n} = \genfrac(){}{}{m}{p^{1}}^{\alpha_1} \genfrac(){}{}{m}{p^{2}}^{\alpha_2} ... \genfrac(){}{}{m}{p^{k}}^{\alpha_k} \text{ where } n = p_1^{\alpha_1} p_2^{\alpha_2} ... p_k^{\alpha_k}\]Like the Legendre symbol, if the Jacobi symbol \(\genfrac(){}{}{m}{n} = 1\) then
m
is a quadratic nonresidue modulon
.But, unlike the Legendre symbol, if the Jacobi symbol \(\genfrac(){}{}{m}{n} = 1\) then
m
may or may not be a quadratic residue modulon
. Parameters
m (integer)
n (odd positive integer)
Examples
>>> jacobi_symbol(45, 77) 1 >>> jacobi_symbol(60, 121) 1
The relationship between the
jacobi_symbol
andlegendre_symbol
can be demonstrated as follows:>>> L = legendre_symbol >>> Integer(45).factors() {3: 2, 5: 1} >>> jacobi_symbol(7, 45) == L(7, 3)**2 * L(7, 5)**1 True
See also

diofant.ntheory.residue_ntheory.
legendre_symbol
(a, p)[source]¶ Returns the Legendre symbol \((a / p)\).
For an integer
a
and an odd primep
, the Legendre symbol is defined as\[\begin{split}\genfrac(){}{}{a}{p} = \begin{cases} 0 & \text{if } p \text{ divides } a\\ 1 & \text{if } a \text{ is a quadratic residue modulo } p\\ 1 & \text{if } a \text{ is a quadratic nonresidue modulo } p \end{cases}\end{split}\] Parameters
a (integer)
p (odd prime)
Examples
>>> [legendre_symbol(i, 7) for i in range(7)] [0, 1, 1, 1, 1, 1, 1] >>> sorted({i**2 % 7 for i in range(7)}) [0, 1, 2, 4]
See also
References

class
diofant.ntheory.residue_ntheory.
mobius
[source]¶ Möbius function maps natural number to {1, 0, 1}
 It is defined as follows:
\(1\) if \(n = 1\).
\(0\) if \(n\) has a squared prime factor.
\((1)^k\) if \(n\) is a squarefree positive integer with \(k\) number of prime factors.
It is an important multiplicative function in number theory and combinatorics. It has applications in mathematical series, algebraic number theory and also physics (Fermion operator has very concrete realization with Möbius Function model).
 Parameters
n (positive integer)
Examples
>>> mobius(13*7) 1 >>> mobius(1) 1 >>> mobius(13*7*5) 1 >>> mobius(13**2) 0
References
Thomas Koshy “Elementary Number Theory with Applications”

diofant.ntheory.residue_ntheory.
n_order
(a, n)[source]¶ Returns the order of
a
modulon
.The order of
a
modulon
is the smallest integerk
such thata**k
leaves a remainder of 1 withn
.Examples
>>> n_order(3, 7) 6 >>> n_order(4, 7) 3

diofant.ntheory.residue_ntheory.
nthroot_mod
(a, n, p, all_roots=False)[source]¶ Find the solutions to
x**n = a mod p
. Parameters
a (integer)
n (positive integer)
p (positive integer)
all_roots (if False returns the smallest root, else the list of roots)
Examples
>>> nthroot_mod(11, 4, 19) 8 >>> nthroot_mod(11, 4, 19, True) [8, 11] >>> nthroot_mod(68, 3, 109) 23

diofant.ntheory.residue_ntheory.
primitive_root
(p)[source]¶ Returns the smallest primitive root or None.
References
Stein “Elementary Number Theory” (2011), page 44
Hackman “Elementary Number Theory” (2009), Chapter C
 Parameters
p (positive integer)
Examples
>>> primitive_root(19) 2

diofant.ntheory.residue_ntheory.
quadratic_residues
(p)[source]¶ Returns the list of quadratic residues.
Examples
>>> quadratic_residues(7) [0, 1, 2, 4]

diofant.ntheory.residue_ntheory.
sqrt_mod
(a, p, all_roots=False)[source]¶ Find a root of
x**2 = a mod p
. Parameters
a (integer)
p (positive integer)
all_roots (if True the list of roots is returned or None)
Notes
If there is no root it is returned None; else the returned root is less or equal to
p // 2
; in general is not the smallest one. It is returnedp // 2
only if it is the only root.Use
all_roots
only when it is expected that all the roots fit in memory; otherwise usesqrt_mod_iter
.Examples
>>> sqrt_mod(11, 43) 21 >>> sqrt_mod(17, 32, True) [7, 9, 23, 25]

diofant.ntheory.residue_ntheory.
sqrt_mod_iter
(a, p, domain=<class 'int'>)[source]¶ Iterate over solutions to
x**2 = a mod p
. Parameters
a (integer)
p (positive integer)
domain (integer domain,
int
,ZZ
orInteger
)
Examples
>>> list(sqrt_mod_iter(11, 43)) [21, 22]

diofant.ntheory.residue_ntheory.
discrete_log
(n, a, b, order=None, prime_order=None)[source] Compute the discrete logarithm of
a
to the baseb
modulon
.This is a recursive function to reduce the discrete logarithm problem in cyclic groups of composite order to the problem in cyclic groups of prime order.
Notes
It employs different algorithms depending on the problem (subgroup order size, prime order or not):
Trial multiplication
Babystep giantstep
PohligHellman
References
Examples
>>> discrete_log(41, 15, 7) 3

diofant.ntheory.continued_fraction.
continued_fraction_convergents
(cf)[source]¶ Return an iterator over the convergents of a continued fraction.
The parameter should be an iterable returning successive partial quotients of the continued fraction, such as might be returned by continued_fraction_iterator. In computing the convergents, the continued fraction need not be strictly in canonical form (all integers, all but the first positive). Rational and negative elements may be present in the expansion.
Examples
>>> list(continued_fraction_convergents([0, 2, 1, 2])) [0, 1/2, 1/3, 3/8]
>>> list(continued_fraction_convergents([1, Rational(1, 2), 7, Rational(1, 4)])) [1, 3, 19/5, 7]
>>> it = continued_fraction_convergents(continued_fraction_iterator(pi)) >>> for n in range(7): ... print(next(it)) 3 22/7 333/106 355/113 103993/33102 104348/33215 208341/66317
See also

diofant.ntheory.continued_fraction.
continued_fraction_iterator
(x)[source]¶ Return continued fraction expansion of x as iterator.
Examples
>>> list(continued_fraction_iterator(Rational(3, 8))) [0, 2, 1, 2] >>> list(continued_fraction_iterator(Rational(3, 8))) [1, 1, 1, 1, 2]
>>> for i, v in enumerate(continued_fraction_iterator(pi)): ... if i > 7: ... break ... print(v) 3 7 15 1 292 1 1 1
References

diofant.ntheory.continued_fraction.
continued_fraction_periodic
(p, q, d=0)[source]¶ Find the periodic continued fraction expansion.
Compute the continued fraction expansion of a rational or a quadratic surd, i.e. \(\frac{p + \sqrt{d}}{q}\), where \(p\), \(q\) and \(d \ge 0\) are integers.
 Returns
list – the continued fraction representation (canonical form) as a list of integers, optionally ending (for quadratic irrationals) with repeating block as the last term of this list.
 Parameters
p (int) – the rational part of the number’s numerator
q (int) – the denominator of the number
d (int, optional) – the irrational part (discriminator) of the number’s numerator
Examples
>>> continued_fraction_periodic(3, 2, 7) [2, [1, 4, 1, 1]]
Golden ratio has the simplest continued fraction expansion:
>>> continued_fraction_periodic(1, 2, 5) [[1]]
If the discriminator is zero or a perfect square then the number will be a rational number:
>>> continued_fraction_periodic(4, 3, 0) [1, 3] >>> continued_fraction_periodic(4, 3, 49) [3, 1, 2]
References
K. Rosen. Elementary Number theory and its applications. AddisonWesley, 3 Sub edition, pages 379381, January 1992.

diofant.ntheory.continued_fraction.
continued_fraction_reduce
(cf)[source]¶ Reduce a continued fraction to a rational or quadratic irrational.
Compute the rational or quadratic irrational number from its terminating or periodic continued fraction expansion. The continued fraction expansion (cf) should be supplied as a terminating iterator supplying the terms of the expansion. For terminating continued fractions, this is equivalent to
list(continued_fraction_convergents(cf))[1]
, only a little more efficient. If the expansion has a repeating part, a list of the repeating terms should be returned as the last element from the iterator. This is the format returned by continued_fraction_periodic.For quadratic irrationals, returns the largest solution found, which is generally the one sought, if the fraction is in canonical form (all terms positive except possibly the first).
Examples
>>> continued_fraction_reduce([1, 2, 3, 4, 5]) 225/157 >>> continued_fraction_reduce([2, 1, 9, 7, 1, 2]) 256/233 >>> continued_fraction_reduce([2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8]).evalf(10) 2.718281835 >>> continued_fraction_reduce([1, 4, 2, [3, 1]]) (sqrt(21) + 287)/238 >>> continued_fraction_reduce([[1]]) 1/2 + sqrt(5)/2 >>> continued_fraction_reduce(continued_fraction_periodic(8, 5, 13)) (sqrt(13) + 8)/5
See also

diofant.ntheory.egyptian_fraction.
egyptian_fraction
(r, algorithm='Greedy')[source]¶ Compute an Egyptian fraction of the rational \(r\).
 Returns
list – The list of denominators of an Egyptian fraction expansion.
 Parameters
r (Rational) – a positive rational number.
algorithm ({ “Greedy”, “Graham Jewett”, “Takenouchi”, “Golomb” }, optional) – Denotes the algorithm to be used (the default is “Greedy”).
Examples
>>> egyptian_fraction(Rational(3, 7)) [3, 11, 231] >>> egyptian_fraction(Rational(3, 7), "Graham Jewett") [7, 8, 9, 56, 57, 72, 3192] >>> egyptian_fraction(Rational(3, 7), "Takenouchi") [4, 7, 28] >>> egyptian_fraction(Rational(3, 7), "Golomb") [3, 15, 35] >>> egyptian_fraction(Rational(11, 5), "Golomb") [1, 2, 3, 4, 9, 234, 1118, 2580]
See also
Notes
Currently the following algorithms are supported:
Greedy Algorithm
Also called the FibonacciSylvester algorithm. At each step, extract the largest unit fraction less than the target and replace the target with the remainder.
It has some distinct properties:
Given \(p/q\) in lowest terms, generates an expansion of maximum length \(p\). Even as the numerators get large, the number of terms is seldom more than a handful.
Uses minimal memory.
The terms can blow up (standard examples of this are 5/121 and 31/311). The denominator is at most squared at each step (doublyexponential growth) and typically exhibits singlyexponential growth.
Graham Jewett Algorithm
The algorithm suggested by the result of Graham and Jewett. Note that this has a tendency to blow up: the length of the resulting expansion is always
2**(x/gcd(x, y))  1
.Takenouchi Algorithm
The algorithm suggested by Takenouchi (1921). Differs from the GrahamJewett algorithm only in the handling of duplicates.
Golomb’s Algorithm
A method given by Golumb (1962), using modular arithmetic and inverses. It yields the same results as a method using continued fractions proposed by Bleicher (1972).
If the given rational is greater than or equal to 1, a greedy algorithm of summing the harmonic sequence 1/1 + 1/2 + 1/3 + … is used, taking all the unit fractions of this sequence until adding one more would be greater than the given number. This list of denominators is prefixed to the result from the requested algorithm used on the remainder. For example, if r is 8/3, using the Greedy algorithm, we get [1, 2, 3, 4, 5, 6, 7, 14, 420], where the beginning of the sequence, [1, 2, 3, 4, 5, 6, 7] is part of the harmonic sequence summing to 363/140, leaving a remainder of 31/420, which yields [14, 420] by the Greedy algorithm. The result of egyptian_fraction(Rational(8, 3), “Golomb”) is [1, 2, 3, 4, 5, 6, 7, 14, 574, 2788, 6460, 11590, 33062, 113820], and so on.
References