Number Theory

Ntheory Class Reference

class sympy.ntheory.generate.Sieve(sieve_interval=1000000)[source]

A 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. Implementation details limit the number of primes to 2^32-1.

Examples

>>> from sympy import sieve
>>> sieve._reset() # this line for doctest only
>>> 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.

Examples

>>> from sympy import sieve
>>> sieve._reset() # this line for doctest only
>>> sieve.extend(30)
>>> sieve[10] == 29
True
extend_to_no(i)[source]

Extend to include the ith prime number.

Parameters:

i : integer

Examples

>>> from sympy import sieve
>>> sieve._reset() # this line for doctest only
>>> sieve.extend_to_no(9)
>>> sieve._list
array('L', [2, 3, 5, 7, 11, 13, 17, 19, 23])

Notes

The list is extended by 50% if it is too short, so it is likely that it will be longer than requested.

mobiusrange(a, b)[source]

Generate all mobius numbers for the range [a, b).

Parameters:

a : integer

First number in range

b : integer

First number outside of range

Examples

>>> from sympy import sieve
>>> print([i for i in sieve.mobiusrange(7, 18)])
[-1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1]
primerange(a, b=None)[source]

Generate all prime numbers in the range [2, a) or [a, b).

Examples

>>> from sympy import sieve, prime

All primes less than 19:

>>> print([i for i in sieve.primerange(19)])
[2, 3, 5, 7, 11, 13, 17]

All primes greater than or equal to 7 and less than 19:

>>> print([i for i in sieve.primerange(7, 19)])
[7, 11, 13, 17]

All primes through the 10th prime

>>> list(sieve.primerange(prime(10) + 1))
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
search(n)[source]

Return the indices i, j of the primes that bound n.

If n is prime then i == j.

Although n can be an expression, if ceiling cannot convert it to an integer then an n error will be raised.

Examples

>>> from sympy import sieve
>>> sieve.search(25)
(9, 10)
>>> sieve.search(23)
(9, 9)
totientrange(a, b)[source]

Generate all totient numbers for the range [a, b).

Examples

>>> from sympy import sieve
>>> print([i for i in sieve.totientrange(7, 18)])
[6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16]
class sympy.ntheory.factor_.FactorCache(maxsize: int | None = None)[source]

Provides a cache for prime factors. factor_cache is pre-prepared as an instance of FactorCache, and factorint internally references it to speed up the factorization of prime factors.

While cache is automatically added during the execution of factorint, users can also manually add prime factors independently.

>>> from sympy import factor_cache
>>> factor_cache[15] = 5

Furthermore, by customizing get_external, it is also possible to use external databases. The following is an example using http://factordb.com .

import requests
from sympy import factor_cache

def get_external(self, n: int) -> list[int] | None:
    res = requests.get("http://factordb.com/api", params={"query": str(n)})
    if res.status_code != requests.codes.ok:
        return None
    j = res.json()
    if j.get("status") in ["FF", "P"]:
        return list(int(p) for p, _ in j.get("factors"))

factor_cache.get_external = get_external

Be aware that writing this code will trigger internet access to factordb.com when calling factorint.

cache_clear() → None[source]

Clear the cache

get(n: int, default=None)[source]

Return the prime factor of n. If it does not exist in the cache, return the value of default.

property maxsize: int | None

Returns the maximum cache size; if None, it is unlimited.

Ntheory Functions Reference

sympy.ntheory.generate.prime(nth)[source]

Return the nth prime number, where primes are indexed starting from 1: prime(1) = 2, prime(2) = 3, etc.

Parameters:

nth : int

The position of the prime number to return (must be a positive integer).

Returns:

int

The nth prime number.

Examples

>>> from sympy import prime
>>> prime(10)
29
>>> prime(1)
2
>>> prime(100000)
1299709

See also

sympy.ntheory.primetest.isprime

Test if a number is prime.

primerange

Generate all primes in a given range.

primepi

Return the number of primes less than or equal to a given number.

References

sympy.ntheory.generate.primepi(n)[source]

Represents the prime counting function pi(n) = the number of prime numbers less than or equal to n.

Deprecated since version 1.13: The primepi function is deprecated. Use sympy.functions.combinatorial.numbers.primepi instead. See its documentation for more information. See Relocate symbolic functions from ntheory to functions for details.

Algorithm Description:

In sieve method, we remove all multiples of prime p except p itself.

Let phi(i,j) be the number of integers 2 <= k <= i which remain after sieving from primes less than or equal to j. Clearly, pi(n) = phi(n, sqrt(n))

If j is not a prime, phi(i,j) = phi(i, j - 1)

if j is a prime, We remove all numbers(except j) whose smallest prime factor is j.

Let \(x= j \times a\) be such a number, where \(2 \le a \le i / j\) Now, after sieving from primes \(\le j - 1\), a must remain (because x, and hence a has no prime factor \(\le j - 1\)) Clearly, there are phi(i / j, j - 1) such a which remain on sieving from primes \(\le j - 1\)

Now, if a is a prime less than equal to j - 1, \(x= j \times a\) has smallest prime factor = a, and has already been removed(by sieving from a). So, we do not need to remove it again. (Note: there will be pi(j - 1) such x)

Thus, number of x, that will be removed are: phi(i / j, j - 1) - phi(j - 1, j - 1) (Note that pi(j - 1) = phi(j - 1, j - 1))

\(\Rightarrow\) phi(i,j) = phi(i, j - 1) - phi(i / j, j - 1) + phi(j - 1, j - 1)

So,following recursion is used and implemented as dp:

phi(a, b) = phi(a, b - 1), if b is not a prime phi(a, b) = phi(a, b-1)-phi(a / b, b-1) + phi(b-1, b-1), if b is prime

Clearly a is always of the form floor(n / k), which can take at most \(2\sqrt{n}\) values. Two arrays arr1,arr2 are maintained arr1[i] = phi(i, j), arr2[i] = phi(n // i, j)

Finally the answer is arr2[1]

Examples

>>> from sympy import primepi, prime, prevprime, isprime
>>> primepi(25)
9

So there are 9 primes less than or equal to 25. Is 25 prime?

>>> isprime(25)
False

It is not. So the first prime less than 25 must be the 9th prime:

>>> prevprime(25) == prime(9)
True

See also

sympy.ntheory.primetest.isprime

Test if n is prime

primerange

Generate all primes in a given range

prime

Return the nth prime

sympy.ntheory.generate.nextprime(n, ith=1)[source]

Return the ith prime greater than n.

Parameters:

n : integer

ith : positive integer

Returns:

int : Return the ith prime greater than n

Raises:

ValueError

If ith <= 0. If n or ith is not an integer.

Notes

Potential primes are located at 6*j +/- 1. This property is used during searching.

>>> from sympy import nextprime
>>> [(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

sympy.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.

>>> from sympy import prevprime
>>> [(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

sympy.ntheory.generate.primerange(a, b=None)[source]

Generate a list of all prime numbers in the range [2, a), or [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.

Examples

>>> from sympy import primerange, prime

All primes less than 19:

>>> list(primerange(19))
[2, 3, 5, 7, 11, 13, 17]

All primes greater than or equal to 7 and less than 19:

>>> list(primerange(7, 19))
[7, 11, 13, 17]

All primes through the 10th prime

>>> list(primerange(prime(10) + 1))
[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:

>>> from sympy import sieve
>>> list(sieve.primerange(1, 30))
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

Notes

Some famous conjectures about the occurrence of primes in a given range are [1]:

  • 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) [2]; 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.

See also

prime

Return the nth prime

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.

References

sympy.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.

Note that due to implementation difficulties, the prime numbers chosen are not uniformly random. For example, there are two primes in the range [112, 128), 113 and 127, but randprime(112, 128) returns 127 with a probability of 15/17.

Examples

>>> from sympy import randprime, isprime
>>> randprime(1, 30) 
13
>>> isprime(randprime(1, 30))
True

See also

primerange

Generate all primes in a given range

References

sympy.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).

Examples

>>> from sympy.ntheory.generate import primorial, primerange
>>> from sympy import factorint, Mul, primefactors, sqrt
>>> 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

sympy.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 value x0.

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.

>>> from sympy.ntheory.generate import cycle_length

This will yield successive values of i <– func(i):

>>> def gen(func, i):
...     while 1:
...         yield i
...         i = func(i)
...

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, 3)

We can see what is meant by looking at the output:

>>> iter = cycle_length(func, 4, values=True)
>>> list(iter)
[4, 17, 35, 2, 5, 26, 14, 44, 50, 2, 5, 26, 14]

There are 6 repeating values after the first 3.

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)
>>> list(cycle_length(func, 4, nmax = 4, values=True))
[4, 17, 35, 2]
Code modified from:

https://en.wikipedia.org/wiki/Cycle_detection.

sympy.ntheory.generate.composite(nth)[source]

Return the nth composite number, with the composite numbers indexed as composite(1) = 4, composite(2) = 6, etc….

Examples

>>> from sympy import composite
>>> composite(36)
52
>>> composite(1)
4
>>> composite(17737)
20000

See also

sympy.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

prime

Return the nth prime

compositepi

Return the number of positive composite numbers less than or equal to n

sympy.ntheory.generate.compositepi(n)[source]

Return the number of positive composite numbers less than or equal to n. The first positive composite is 4, i.e. compositepi(4) = 1.

Examples

>>> from sympy import compositepi
>>> compositepi(25)
15
>>> compositepi(1000)
831

See also

sympy.ntheory.primetest.isprime

Test if n is prime

primerange

Generate all primes in a given range

prime

Return the nth prime

primepi

Return the number of primes less than or equal to n

composite

Return the nth composite number

sympy.ntheory.factor_.smoothness(n)[source]

Return the B-smooth and B-power 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.

Examples

>>> from sympy.ntheory.factor_ import smoothness
>>> smoothness(2**7*3**2)
(3, 128)
>>> smoothness(2**4*13)
(13, 16)
>>> smoothness(2)
(2, 2)

See also

factorint, smoothness_p

sympy.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:

  1. p**M is the base-p divisor of n

  2. sm(p + m) is the smoothness of p + m (m = -1 by default)

  3. 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.

>>> from sympy.ntheory.factor_ import smoothness_p, factorint
>>> 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 p-1 B=1787, B-pow=1787
p**i=4869863**1 has p-1 B=2434931, B-pow=2434931

This string can also be generated directly from a factorization dictionary and vice versa:

>>> factorint(17*9)
{3: 2, 17: 1}
>>> smoothness_p(_)
'p**i=3**2 has p-1 B=2, B-pow=2\np**i=17**1 has p-1 B=2, B-pow=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

factorint, smoothness

sympy.ntheory.factor_.multiplicity(p, n)[source]

Find the greatest integer m such that p**m divides n.

Examples

>>> from sympy import multiplicity, Rational
>>> [multiplicity(5, n) for n in [8, 5, 25, 125, 250]]
[0, 1, 2, 3, 3]
>>> multiplicity(3, Rational(1, 9))
-2

Note: when checking for the multiplicity of a number in a large factorial it is most efficient to send it as an unevaluated factorial or to call multiplicity_in_factorial directly:

>>> from sympy.ntheory import multiplicity_in_factorial
>>> from sympy import factorial
>>> p = factorial(25)
>>> n = 2**100
>>> nfac = factorial(n, evaluate=False)
>>> multiplicity(p, nfac)
52818775009509558395695966887
>>> _ == multiplicity_in_factorial(p, n)
True

See also

trailing

sympy.ntheory.factor_.perfect_power(
n,
candidates=None,
big=True,
factor=True,
)[source]

Return (b, e) such that n == b**e if n is a unique perfect power with e > 1, else False (e.g. 1 is not a perfect power). A ValueError is raised if n is not Rational.

By default, the base is recursively decomposed and the exponents collected so the largest possible e is sought. If big=False then the smallest possible e (thus prime) will be chosen.

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.

The use of candidates is primarily for internal use; if provided, False will be returned if n cannot be written as a power with one of the candidates as an exponent and factoring (beyond testing for a factor of 2) will not be attempted.

Examples

>>> from sympy import perfect_power, Rational
>>> perfect_power(16)
(2, 4)
>>> perfect_power(16, big=False)
(4, 2)

Negative numbers can only have odd perfect powers:

>>> perfect_power(-4)
False
>>> perfect_power(-8)
(-2, 3)

Rationals are also recognized:

>>> perfect_power(Rational(1, 2)**3)
(1/2, 3)
>>> perfect_power(Rational(-3, 2)**3)
(-3/2, 3)

Notes

To know whether an integer is a perfect power of 2 use

>>> is2pow = lambda n: bool(n and not n & (n - 1))
>>> [(i, is2pow(i)) for i in range(5)]
[(0, False), (1, True), (2, True), (3, False), (4, True)]

It is not necessary to provide candidates. When provided it will be assumed that they are ints. The first one that is larger than the computed maximum possible exponent will signal failure for the routine.

>>> perfect_power(3**8, [9])
False
>>> perfect_power(3**8, [2, 4, 8])
(3, 8)
>>> perfect_power(3**8, [4, 8], big=False)
(9, 4)

See also

sympy.core.intfunc.integer_nthroot, sympy.ntheory.primetest.is_square

sympy.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 pseudo-random 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) is s. Upon failure (if retries is > 0) a new a and s will be supplied; the a will be ignored if F was supplied.

The sequence of numbers generated by such functions generally have a a lead-up 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 leader-loop values can be obtained so it is a good idea to allow for retries:

>>> from sympy.ntheory.generate import cycle_length
>>> 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 =  43
loop length =   78; leader length = 121
loop length = 1482; leader length = 100
loop length = 1482; leader length = 286
loop length = 1482; leader length = 101

Here is an explicit example where there is a three element leadup to a sequence of 3 numbers (11, 14, 4) that then repeat:

>>> x=2
>>> for i in range(9):
...     print(x)
...     x=(x**2+12)%17
...
2
16
13
11
14
4
11
14
4
>>> next(cycle_length(lambda x: (x**2+12)%17, 2))
(3, 3)
>>> list(cycle_length(lambda x: (x**2+12)%17, 2, values=True))
[2, 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

>>> from sympy import pollard_rho
>>> 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=n-2, retries=0)

If retries is > 0 then perhaps the problem will correct itself when new values are generated for a:

>>> pollard_rho(n, a=n-2, retries=1)
257

References

[R654]

Richard Crandall & Carl Pomerance (2005), “Prime Numbers: A Computational Perspective”, Springer, 2nd edition, 229-231

sympy.ntheory.factor_.pollard_pm1(n, B=10, a=2, retries=0, seed=1234)[source]

Use Pollard’s p-1 method to try to extract a nontrivial factor of n. Either a divisor (perhaps composite) or None is returned.

The value of a is the base that is used in the test gcd(a**M - 1, n). The default is 2. If retries > 0 then if no factor is found after the first attempt, a new a will be generated randomly (using the seed) 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 on a and the power smoothness of the even number 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:

>>> from sympy.ntheory.factor_ import smoothness_p, pollard_pm1
>>> 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 does not 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 did not work had a gcd value not equal to n but equal to one of the factors:

>>> from sympy import ilcm, igcd, factorint, Pow
>>> M = 1
>>> for i in range(2, 256):
...     M = ilcm(M, i)
...
>>> set([igcd(pow(a, M, n) - 1, n) for a in range(2, 256) if
...      igcd(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:

  1. the power smoothness of the p - 1 value next to the root does not exceed B

  2. 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 cannot be cracked:

>>> from sympy.ntheory import pollard_pm1
>>> 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:

>>> from sympy.ntheory.factor_ import smoothness_p, factorint
>>> print(smoothness_p(21477639576571, visual=1))
p**i=4410317**1 has p-1 B=1787, B-pow=1787
p**i=4869863**1 has p-1 B=2434931, B-pow=2434931

The B and B-pow 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(4869863-1)
{2: 1, 2434931: 1}

Note that until B reaches the B-pow 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 power-smoothness 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 do not 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

[R655]

Richard Crandall & Carl Pomerance (2005), “Prime Numbers: A Computational Perspective”, Springer, 2nd edition, 236-238

sympy.ntheory.factor_.factorint(
n,
limit=None,
use_trial=True,
use_rho=True,
use_pm1=True,
use_ecm=True,
verbose=False,
visual=None,
multiple=False,
)[source]

Given a positive integer n, factorint(n) returns a dict containing the prime factors of n as keys and their respective multiplicities as values. For example:

>>> from sympy.ntheory import factorint
>>> 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) adds -1:1 to the factors and then factors n

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/p-1 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 semi-prime factor that cannot be reduced easily:

>>> from sympy.ntheory import isprime
>>> a = 1407633717262338957430697921446883
>>> f = factorint(a, limit=10000)
>>> f == {991: 1, int(202916782076162456022877024859): 1, 7: 1}
True
>>> 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}

List of Factors:

If multiple is set to True then a list containing the prime factors including multiplicities is returned.

>>> factorint(24, multiple=True)
[2, 2, 2, 3]

Visual Factorization:

If visual is set to True, then it will return a visual factorization of the integer. For example:

>>> from sympy import pprint
>>> pprint(factorint(4200, visual=True))
 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:

>>> from sympy import Mul
>>> regular = factorint(1764); regular
{2: 2, 3: 2, 7: 2}
>>> pprint(factorint(regular))
 2  2  2
2 *3 *7

>>> visual = factorint(1764, visual=True); pprint(visual)
 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

Algorithm:

The function switches between multiple algorithms. Trial division quickly finds small factors (of the order 1-5 digits), and finds all large factors if given enough time. The Pollard rho and p-1 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 division

  • use_rho: Toggle use of Pollard’s rho method

  • use_pm1: Toggle use of Pollard’s p-1 method

factorint also periodically checks if the remaining part is a prime number or a perfect power, and in those cases stops.

For unevaluated factorial, it uses Legendre’s formula(theorem).

If verbose is set to True, detailed progress is printed.

See also

smoothness, smoothness_p, divisors

sympy.ntheory.factor_.factorrat(
rat,
limit=None,
use_trial=True,
use_rho=True,
use_pm1=True,
verbose=False,
visual=None,
multiple=False,
)[source]

Given a Rational r, factorrat(r) returns a dict containing the prime factors of r as keys and their respective multiplicities as values. For example:

>>> from sympy import factorrat, S
>>> factorrat(S(8)/9)    # 8/9 = (2**3) * (3**-2)
{2: 3, 3: -2}
>>> factorrat(S(-1)/987)    # -1/789 = -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 done

  • use_trial: Toggle use of trial division

  • use_rho: Toggle use of Pollard’s rho method

  • use_pm1: Toggle use of Pollard’s p-1 method

  • verbose: Toggle detailed printing of progress

  • multiple: Toggle returning a list of factors or dict

  • visual: Toggle product form of output

sympy.ntheory.factor_.primefactors(
n,
limit=None,
verbose=False,
**kwargs,
)[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.

Parameters:

n : integer

limit, verbose, **kwargs :

Additional keyword arguments to be passed to factorint. Since kwargs is new in version 1.13, limit and verbose are retained for compatibility purposes.

Returns:

list(int) : List of prime numbers dividing n

Examples

>>> from sympy.ntheory import primefactors, factorint, isprime
>>> 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

factorint, divisors

sympy.ntheory.factor_.divisors(n, generator=False, proper=False)[source]

Return all divisors of n 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

>>> from sympy import divisors, divisor_count
>>> 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]

Notes

This is a slightly modified version of Tim Peters referenced at: https://stackoverflow.com/questions/1010381/python-factorization

See also

primefactors, factorint, divisor_count

sympy.ntheory.factor_.proper_divisors(n, generator=False)[source]

Return all divisors of n except n, sorted by default. If generator is True an unordered generator is returned.

Examples

>>> from sympy import proper_divisors, proper_divisor_count
>>> proper_divisors(24)
[1, 2, 3, 4, 6, 8, 12]
>>> proper_divisor_count(24)
7
>>> list(proper_divisors(120, generator=True))
[1, 2, 4, 8, 3, 6, 12, 24, 5, 10, 20, 40, 15, 30, 60]

See also

factorint, divisors, proper_divisor_count

sympy.ntheory.factor_.divisor_count(n, modulus=1, proper=False)[source]

Return the number of divisors of n. If modulus is not 1 then only those that are divisible by modulus are counted. If proper is True then the divisor of n will not be counted.

Examples

>>> from sympy import divisor_count
>>> divisor_count(6)
4
>>> divisor_count(6, 2)
2
>>> divisor_count(6, proper=True)
3

See also

factorint, divisors, totient, proper_divisor_count

sympy.ntheory.factor_.proper_divisor_count(n, modulus=1)[source]

Return the number of proper divisors of n.

Examples

>>> from sympy import proper_divisor_count
>>> proper_divisor_count(6)
3
>>> proper_divisor_count(6, modulus=2)
1

See also

divisors, proper_divisors, divisor_count

sympy.ntheory.factor_.udivisors(n, generator=False)[source]

Return all unitary divisors of n sorted from 1..n by default. If generator is True an unordered generator is returned.

The number of unitary divisors of n can be quite large if there are many prime factors. If only the number of unitary divisors is desired use udivisor_count(n).

Examples

>>> from sympy.ntheory.factor_ import udivisors, udivisor_count
>>> udivisors(15)
[1, 3, 5, 15]
>>> udivisor_count(15)
4

>>> sorted(udivisors(120, generator=True))
[1, 3, 5, 8, 15, 24, 40, 120]

See also

primefactors, factorint, divisors, divisor_count, udivisor_count

References

sympy.ntheory.factor_.udivisor_count(n)[source]

Return the number of unitary divisors of n.

Parameters:

n : integer

Examples

>>> from sympy.ntheory.factor_ import udivisor_count
>>> udivisor_count(120)
8

See also

factorint, divisors, udivisors, divisor_count, totient

References

sympy.ntheory.factor_.antidivisors(n, generator=False)[source]

Return all antidivisors of n sorted from 1..n by default.

Antidivisors [R661] of n are numbers that do not divide n by the largest possible margin. If generator is True an unordered generator is returned.

Examples

>>> from sympy.ntheory.factor_ import antidivisors
>>> antidivisors(24)
[7, 16]

>>> sorted(antidivisors(128, generator=True))
[3, 5, 15, 17, 51, 85]

See also

primefactors, factorint, divisors, divisor_count, antidivisor_count

References

[R661] (1,2)

definition is described in https://oeis.org/A066272/a066272a.html

sympy.ntheory.factor_.antidivisor_count(n)[source]

Return the number of antidivisors [R662] of n.

Parameters:

n : integer

Examples

>>> from sympy.ntheory.factor_ import antidivisor_count
>>> antidivisor_count(13)
4
>>> antidivisor_count(27)
5

See also

factorint, divisors, antidivisors, divisor_count, totient

References

[R662] (1,2)

formula from https://oeis.org/A066272

sympy.ntheory.factor_.totient(n)[source]

Calculate the Euler totient function phi(n)

Deprecated since version 1.13: The totient function is deprecated. Use sympy.functions.combinatorial.numbers.totient instead. See its documentation for more information. See Relocate symbolic functions from ntheory to functions for details.

totient(n) or \(\phi(n)\) is the number of positive integers \(\leq\) n that are relatively prime to n.

Parameters:

n : integer

Examples

>>> from sympy.functions.combinatorial.numbers import totient
>>> totient(1)
1
>>> totient(25)
20
>>> totient(45) == totient(5)*totient(9)
True

See also

divisor_count

References

sympy.ntheory.factor_.reduced_totient(n)[source]

Calculate the Carmichael reduced totient function lambda(n)

Deprecated since version 1.13: The reduced_totient function is deprecated. Use sympy.functions.combinatorial.numbers.reduced_totient instead. See its documentation for more information. See Relocate symbolic functions from ntheory to functions for details.

reduced_totient(n) or \(\lambda(n)\) is the smallest m > 0 such that \(k^m \equiv 1 \mod n\) for all k relatively prime to n.

Examples

>>> from sympy.functions.combinatorial.numbers import reduced_totient
>>> reduced_totient(1)
1
>>> reduced_totient(8)
2
>>> reduced_totient(30)
4

See also

totient

References

sympy.ntheory.factor_.divisor_sigma(n, k=1)[source]

Calculate the divisor function \(\sigma_k(n)\) for positive integer n

Deprecated since version 1.13: The divisor_sigma function is deprecated. Use sympy.functions.combinatorial.numbers.divisor_sigma instead. See its documentation for more information. See Relocate symbolic functions from ntheory to functions for details.

divisor_sigma(n, k) is equal to sum([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:

n : integer

k : integer, optional

power of divisors in the sum

for k = 0, 1: divisor_sigma(n, 0) is equal to divisor_count(n) divisor_sigma(n, 1) is equal to sum(divisors(n))

Default for k is 1.

Examples

>>> from sympy.functions.combinatorial.numbers import divisor_sigma
>>> divisor_sigma(18, 0)
6
>>> divisor_sigma(39, 1)
56
>>> divisor_sigma(12, 2)
210
>>> divisor_sigma(37)
38

See also

divisor_count, totient, divisors, factorint

References

sympy.ntheory.factor_.udivisor_sigma(n, k=1)[source]

Calculate the unitary divisor function \(\sigma_k^*(n)\) for positive integer n

Deprecated since version 1.13: The udivisor_sigma function is deprecated. Use sympy.functions.combinatorial.numbers.udivisor_sigma instead. See its documentation for more information. See Relocate symbolic functions from ntheory to functions for details.

udivisor_sigma(n, k) is equal to sum([x**k for x in udivisors(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^{m_ik}).\]
Parameters:

k : power of divisors in the sum

for k = 0, 1: udivisor_sigma(n, 0) is equal to udivisor_count(n) udivisor_sigma(n, 1) is equal to sum(udivisors(n))

Default for k is 1.

Examples

>>> from sympy.functions.combinatorial.numbers import udivisor_sigma
>>> udivisor_sigma(18, 0)
4
>>> udivisor_sigma(74, 1)
114
>>> udivisor_sigma(36, 3)
47450
>>> udivisor_sigma(111)
152

See also

divisor_count, totient, divisors, udivisors, udivisor_count, divisor_sigma, factorint

References

sympy.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 n

If 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:

n : integer

t : integer

core(n, t) calculates the t-th power free part of n

core(n, 2) is the squarefree part of n core(n, 3) is the cubefree part of n

Default for t is 2.

Examples

>>> from sympy.ntheory.factor_ import core
>>> core(24, 2)
6
>>> core(9424, 3)
1178
>>> core(379238)
379238
>>> core(15**11, 10)
15

See also

factorint, sympy.solvers.diophantine.diophantine.square_factor

References

sympy.ntheory.factor_.digits(n, b=10, digits=None)[source]

Return a list of the digits of n in base b. The first element in the list is b (or -b if n is negative).

Parameters:

n: integer

The number whose digits are returned.

b: integer

The base in which digits are computed.

digits: integer (or None for all digits)

The number of digits to be returned (padded with zeros, if necessary).

Examples

>>> from sympy.ntheory.digits import digits
>>> digits(35)
[10, 3, 5]

If the number is negative, the negative sign will be placed on the base (which is the first element in the returned list):

>>> digits(-35)
[-10, 3, 5]

Bases other than 10 (and greater than 1) can be selected with b:

>>> digits(27, b=2)
[2, 1, 1, 0, 1, 1]

Use the digits keyword if a certain number of digits is desired:

>>> digits(35, digits=4)
[10, 0, 0, 3, 5]

See also

sympy.core.intfunc.num_digits, count_digits

sympy.ntheory.factor_.primenu(n)[source]

Calculate the number of distinct prime factors for a positive integer n.

Deprecated since version 1.13: The primenu function is deprecated. Use sympy.functions.combinatorial.numbers.primenu instead. See its documentation for more information. See Relocate symbolic functions from ntheory to functions for details.

If n’s prime factorization is:

\[n = \prod_{i=1}^k p_i^{m_i},\]

then primenu(n) or \(\nu(n)\) is:

\[\nu(n) = k.\]

Examples

>>> from sympy.functions.combinatorial.numbers import primenu
>>> primenu(1)
0
>>> primenu(30)
3

See also

factorint

References

sympy.ntheory.factor_.primeomega(n)[source]

Calculate the number of prime factors counting multiplicities for a positive integer n.

Deprecated since version 1.13: The primeomega function is deprecated. Use sympy.functions.combinatorial.numbers.primeomega instead. See its documentation for more information. See Relocate symbolic functions from ntheory to functions for details.

If n’s prime factorization is:

\[n = \prod_{i=1}^k p_i^{m_i},\]

then primeomega(n) or \(\Omega(n)\) is:

\[\Omega(n) = \sum_{i=1}^k m_i.\]

Examples

>>> from sympy.functions.combinatorial.numbers import primeomega
>>> primeomega(1)
0
>>> primeomega(20)
3

See also

factorint

References

sympy.ntheory.factor_.mersenne_prime_exponent(nth)[source]

Returns the exponent i for the nth Mersenne prime (which has the form \(2^i - 1\)).

Examples

>>> from sympy.ntheory.factor_ import mersenne_prime_exponent
>>> mersenne_prime_exponent(1)
2
>>> mersenne_prime_exponent(20)
4423
sympy.ntheory.factor_.is_perfect(n)[source]

Returns True if n is a perfect number, else False.

A perfect number is equal to the sum of its positive, proper divisors.

Examples

>>> from sympy.functions.combinatorial.numbers import divisor_sigma
>>> from sympy.ntheory.factor_ import is_perfect, divisors
>>> is_perfect(20)
False
>>> is_perfect(6)
True
>>> 6 == divisor_sigma(6) - 6 == sum(divisors(6)[:-1])
True

References

sympy.ntheory.factor_.abundance(n)[source]

Returns the difference between the sum of the positive proper divisors of a number and the number.

Examples

>>> from sympy.ntheory import abundance, is_perfect, is_abundant
>>> abundance(6)
0
>>> is_perfect(6)
True
>>> abundance(10)
-2
>>> is_abundant(10)
False
sympy.ntheory.factor_.is_abundant(n)[source]

Returns True if n is an abundant number, else False.

A abundant number is smaller than the sum of its positive proper divisors.

Examples

>>> from sympy.ntheory.factor_ import is_abundant
>>> is_abundant(20)
True
>>> is_abundant(15)
False

References

sympy.ntheory.factor_.is_deficient(n)[source]

Returns True if n is a deficient number, else False.

A deficient number is greater than the sum of its positive proper divisors.

Examples

>>> from sympy.ntheory.factor_ import is_deficient
>>> is_deficient(20)
False
>>> is_deficient(15)
True

References

sympy.ntheory.factor_.is_amicable(m, n)[source]

Returns True if the numbers \(m\) and \(n\) are “amicable”, else False.

Amicable numbers are two different numbers so related that the sum of the proper divisors of each is equal to that of the other.

Examples

>>> from sympy.functions.combinatorial.numbers import divisor_sigma
>>> from sympy.ntheory.factor_ import is_amicable
>>> is_amicable(220, 284)
True
>>> divisor_sigma(220) == divisor_sigma(284)
True

References

sympy.ntheory.factor_.is_carmichael(n)[source]

Returns True if the numbers \(n\) is Carmichael number, else False.

Parameters:

n : Integer

References

sympy.ntheory.factor_.find_carmichael_numbers_in_range(x, y)[source]

Returns a list of the number of Carmichael in the range

See also

is_carmichael

sympy.ntheory.factor_.find_first_n_carmichaels(n)[source]

Returns the first n Carmichael numbers.

Parameters:

n : Integer

See also

is_carmichael

sympy.ntheory.modular.symmetric_residue(a, m)[source]

Return the residual mod m such that it is within half of the modulus.

>>> from sympy.ntheory.modular import symmetric_residue
>>> symmetric_residue(1, 6)
1
>>> symmetric_residue(4, 6)
-2
sympy.ntheory.modular.crt(m, v, 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 = v_i mod m_i for each pair out of v 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 co-prime 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

As an example consider a set of residues U = [49, 76, 65] and a set of moduli M = [99, 97, 95]. Then we have:

>>> from sympy.ntheory.modular import crt

>>> 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 co-prime, 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)

Note: the order of gf_crt’s arguments is reversed relative to crt, and that solve_congruence takes residue, modulus pairs.

Programmer’s note: 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

solve_congruence

sympy.polys.galoistools.gf_crt

low level crt routine used by this routine

sympy.ntheory.modular.crt1(m)[source]

First part of Chinese Remainder Theorem, for multiple application.

Examples

>>> from sympy.ntheory.modular import crt, crt1, crt2
>>> m = [99, 97, 95]
>>> v = [49, 76, 65]

The following two codes have the same result.

>>> crt(m, v)
(639985, 912285)

>>> mm, e, s = crt1(m)
>>> crt2(m, v, mm, e, s)
(639985, 912285)

However, it is faster when we want to fix m and compute for multiple v, i.e. the following cases:

>>> mm, e, s = crt1(m)
>>> vs = [[52, 21, 37], [19, 46, 76]]
>>> for v in vs:
...     print(crt2(m, v, mm, e, s))
(397042, 912285)
(803206, 912285)

See also

sympy.polys.galoistools.gf_crt1

low level crt routine used by this routine

sympy.ntheory.modular.crt, sympy.ntheory.modular.crt2

sympy.ntheory.modular.crt2(m, v, mm, e, s, symmetric=False)[source]

Second part of Chinese Remainder Theorem, for multiple application.

See crt1 for usage.

Examples

>>> from sympy.ntheory.modular import crt1, crt2
>>> mm, e, s = crt1([18, 42, 6])
>>> crt2([18, 42, 6], [0, 0, 0], mm, e, s)
(0, 4536)

See also

sympy.polys.galoistools.gf_crt2

low level crt routine used by this routine

sympy.ntheory.modular.crt, sympy.ntheory.modular.crt1

sympy.ntheory.modular.solve_congruence(*remainder_modulus_pairs, **hint)[source]

Compute the integer n that has the residual ai when it is divided by mi where the ai and mi are given as pairs to this function: ((a1, m1), (a2, m2), …). If there is no solution, return None. Otherwise return n and its modulus.

The mi values need not be co-prime. If it is known that the moduli are not co-prime then the hint check can be set to False (default=True) and the check for a quicker solution via crt() (valid when the moduli are co-prime) will be skipped.

If the hint symmetric is True (default is False), the value of n will be within 1/2 of the modulus, possibly negative.

Examples

>>> from sympy.ntheory.modular import solve_congruence

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 co-prime; 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

sympy.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

>>> from sympy.ntheory import binomial_coefficients
>>> 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}

See also

binomial_coefficients_list, multinomial_coefficients

sympy.ntheory.multinomial.binomial_coefficients_list(n)[source]

Return a list of binomial coefficients as rows of the Pascal’s triangle.

Examples

>>> from sympy.ntheory import binomial_coefficients_list
>>> binomial_coefficients_list(9)
[1, 9, 36, 84, 126, 126, 84, 36, 9, 1]

See also

binomial_coefficients, multinomial_coefficients

sympy.ntheory.multinomial.multinomial_coefficients(m, n)[source]

Return a dictionary containing pairs {(k1,k2,..,km) : C_kn} where C_kn are multinomial coefficients such that n=k1+k2+..+km.

Examples

>>> from sympy.ntheory import multinomial_coefficients
>>> multinomial_coefficients(2, 5) # indirect doctest
{(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}\]

Code contributed to Sage by Yann Laigle-Chapuy, copied with permission of the author.

See also

binomial_coefficients_list, binomial_coefficients

sympy.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.

>>> from sympy.ntheory.multinomial import multinomial_coefficients
>>> 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

>>> from sympy.ntheory.multinomial import multinomial_coefficients_iterator
>>> 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)
sympy.ntheory.partitions_.npartitions(n, verbose=False)[source]

Calculate the partition function P(n), i.e. the number of ways that n can be written as a sum of positive integers.

Deprecated since version 1.13: The npartitions function is deprecated. Use sympy.functions.combinatorial.numbers.partition instead. See its documentation for more information. See Relocate symbolic functions from ntheory to functions for details.

P(n) is computed using the Hardy-Ramanujan-Rademacher formula [R679].

The correctness of this implementation has been tested through \(10^{10}\).

Examples

>>> from sympy.functions.combinatorial.numbers import partition
>>> partition(25)
1958

References

sympy.ntheory.primetest.is_fermat_pseudoprime(n, a)[source]

Returns True if n is prime or is an odd composite integer that is coprime to a and satisfy the modular arithmetic congruence relation:

\[a^{n-1} \equiv 1 \pmod{n}\]

(where mod refers to the modulo operation).

Parameters:

n : Integer

n is a positive integer.

a : Integer

a is a positive integer. a and n should be relatively prime.

Returns:

bool : If n is prime, it always returns True.

The composite number that returns True is called an Fermat pseudoprime.

Examples

>>> from sympy.ntheory.primetest import is_fermat_pseudoprime
>>> from sympy.ntheory.factor_ import isprime
>>> for n in range(1, 1000):
...     if is_fermat_pseudoprime(n, 2) and not isprime(n):
...         print(n)
341
561
645

References

sympy.ntheory.primetest.is_euler_pseudoprime(n, a)[source]

Returns True if n is prime or is an odd composite integer that is coprime to a and satisfy the modular arithmetic congruence relation:

\[a^{(n-1)/2} \equiv \pm 1 \pmod{n}\]

(where mod refers to the modulo operation).

Parameters:

n : Integer

n is a positive integer.

a : Integer

a is a positive integer. a and n should be relatively prime.

Returns:

bool : If n is prime, it always returns True.

The composite number that returns True is called an Euler pseudoprime.

Examples

>>> from sympy.ntheory.primetest import is_euler_pseudoprime
>>> from sympy.ntheory.factor_ import isprime
>>> for n in range(1, 1000):
...     if is_euler_pseudoprime(n, 2) and not isprime(n):
...         print(n)
341
561

References

sympy.ntheory.primetest.is_euler_jacobi_pseudoprime(n, a)[source]

Returns True if n is prime or is an odd composite integer that is coprime to a and satisfy the modular arithmetic congruence relation:

\[a^{(n-1)/2} \equiv \left(\frac{a}{n}\right) \pmod{n}\]

(where mod refers to the modulo operation).

Parameters:

n : Integer

n is a positive integer.

a : Integer

a is a positive integer. a and n should be relatively prime.

Returns:

bool : If n is prime, it always returns True.

The composite number that returns True is called an Euler-Jacobi pseudoprime.

Examples

>>> from sympy.ntheory.primetest import is_euler_jacobi_pseudoprime
>>> from sympy.ntheory.factor_ import isprime
>>> for n in range(1, 1000):
...     if is_euler_jacobi_pseudoprime(n, 2) and not isprime(n):
...         print(n)
561

References

sympy.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.

Examples

>>> from sympy.ntheory.primetest import is_square
>>> is_square(25)
True
>>> is_square(2)
False

See also

sympy.core.intfunc.isqrt

References

sympy.ntheory.primetest.mr(n, bases)[source]

Perform a Miller-Rabin strong pseudoprime test on n using a given list of bases/witnesses.

Examples

>>> from sympy.ntheory.primetest import mr
>>> mr(1373651, [2, 3])
False
>>> mr(479001599, [31, 73])
True

References

[R684]

Richard Crandall & Carl Pomerance (2005), “Prime Numbers: A Computational Perspective”, Springer, 2nd edition, 135-138

A list of thresholds and the bases they require are here: https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Deterministic_variants

sympy.ntheory.primetest.is_lucas_prp(n)[source]

Standard Lucas compositeness test with Selfridge parameters. Returns False if n is definitely composite, and True if n is a Lucas probable prime.

This is typically used in combination with the Miller-Rabin test.

Examples

>>> from sympy.ntheory.primetest import isprime, is_lucas_prp
>>> for i in range(10000):
...     if is_lucas_prp(i) and not isprime(i):
...         print(i)
323
377
1159
1829
3827
5459
5777
9071
9179

References

[R685]

Robert Baillie, Samuel S. Wagstaff, Lucas Pseudoprimes, Math. Comp. Vol 35, Number 152 (1980), pp. 1391-1417, https://doi.org/10.1090%2FS0025-5718-1980-0583518-6 http://mpqs.free.fr/LucasPseudoprimes.pdf

[R686]

OEIS A217120: Lucas Pseudoprimes https://oeis.org/A217120

sympy.ntheory.primetest.is_strong_lucas_prp(n)[source]

Strong Lucas compositeness test with Selfridge parameters. Returns False if n is definitely composite, and True if n is a strong Lucas probable prime.

This is often used in combination with the Miller-Rabin test, and in particular, when combined with M-R base 2 creates the strong BPSW test.

Examples

>>> from sympy.ntheory.primetest import isprime, is_strong_lucas_prp
>>> for i in range(20000):
...     if is_strong_lucas_prp(i) and not isprime(i):
...        print(i)
5459
5777
10877
16109
18971

References

[R688]

Robert Baillie, Samuel S. Wagstaff, Lucas Pseudoprimes, Math. Comp. Vol 35, Number 152 (1980), pp. 1391-1417, https://doi.org/10.1090%2FS0025-5718-1980-0583518-6 http://mpqs.free.fr/LucasPseudoprimes.pdf

[R689]

OEIS A217255: Strong Lucas Pseudoprimes https://oeis.org/A217255

sympy.ntheory.primetest.is_extra_strong_lucas_prp(n)[source]

Extra Strong Lucas compositeness test. Returns False if n is definitely composite, and True if n is an “extra strong” Lucas probable prime.

The parameters are selected using P = 3, Q = 1, then incrementing P until (D|n) == -1. The test itself is as defined in [R692], from the Mo and Jones preprint. The parameter selection and test are the same as used in OEIS A217719, Perl’s Math::Prime::Util, and the Lucas pseudoprime page on Wikipedia.

It is 20-50% faster than the strong test.

Because of the different parameters selected, there is no relationship between the strong Lucas pseudoprimes and extra strong Lucas pseudoprimes. In particular, one is not a subset of the other.

Examples

>>> from sympy.ntheory.primetest import isprime, is_extra_strong_lucas_prp
>>> for i in range(20000):
...     if is_extra_strong_lucas_prp(i) and not isprime(i):
...        print(i)
989
3239
5777
10877

References

[R692] (1,2)

Jon Grantham, Frobenius Pseudoprimes, Math. Comp. Vol 70, Number 234 (2001), pp. 873-891, https://doi.org/10.1090%2FS0025-5718-00-01197-2

[R693]

OEIS A217719: Extra Strong Lucas Pseudoprimes https://oeis.org/A217719

sympy.ntheory.primetest.proth_test(n)[source]

Test if the Proth number \(n = k2^m + 1\) is prime. where k is a positive odd number and \(2^m > k\).

Parameters:

n : Integer

n is Proth number

Returns:

bool : If True, then n is the Proth prime

Raises:

ValueError

If n is not Proth number.

Examples

>>> from sympy.ntheory.primetest import proth_test
>>> proth_test(41)
True
>>> proth_test(57)
False

References

sympy.ntheory.primetest.is_mersenne_prime(n)[source]

Returns True if n is a Mersenne prime, else False.

A Mersenne prime is a prime number having the form \(2^i - 1\).

Examples

>>> from sympy.ntheory.factor_ import is_mersenne_prime
>>> is_mersenne_prime(6)
False
>>> is_mersenne_prime(127)
True

References

sympy.ntheory.primetest.isprime(n)[source]

Test if n is a prime number (True) or not (False). For n < 2^64 the answer is definitive; larger n values have a small probability of actually being pseudoprimes.

Negative numbers (e.g. -2) are not considered prime.

The first step is looking for trivial factors, which if found enables a quick return. Next, if the sieve is large enough, use bisection search on the sieve. For small numbers, a set of deterministic Miller-Rabin tests are performed with bases that are known to have no counterexamples in their range. Finally if the number is larger than 2^64, a strong BPSW test is performed. While this is a probable prime test and we believe counterexamples exist, there are no known counterexamples.

Examples

>>> from sympy.ntheory import isprime
>>> isprime(13)
True
>>> isprime(15)
False

Notes

This routine is intended only for integer input, not numerical expressions which may represent numbers. Floats are also rejected as input because they represent numbers of limited precision. While it is tempting to permit 7.0 to represent an integer there are errors that may “pass silently” if this is allowed:

>>> from sympy import Float, S
>>> int(1e3) == 1e3 == 10**3
True
>>> int(1e23) == 1e23
True
>>> int(1e23) == 10**23
False

>>> near_int = 1 + S(1)/10**19
>>> near_int == int(near_int)
False
>>> n = Float(near_int, 10)  # truncated by precision
>>> n % 1 == 0
True
>>> n = Float(near_int, 20)
>>> n % 1 == 0
False

See also

sympy.ntheory.generate.primerange

Generates all primes in a given range

sympy.functions.combinatorial.numbers.primepi

Return the number of primes less than or equal to n

sympy.ntheory.generate.prime

Return the nth prime

References

[R698]

Robert Baillie, Samuel S. Wagstaff, Lucas Pseudoprimes, Math. Comp. Vol 35, Number 152 (1980), pp. 1391-1417, https://doi.org/10.1090%2FS0025-5718-1980-0583518-6 http://mpqs.free.fr/LucasPseudoprimes.pdf

sympy.ntheory.primetest.is_gaussian_prime(num)[source]

Test if num is a Gaussian prime number.

References

sympy.ntheory.residue_ntheory.n_order(a, n)[source]

Returns the order of a modulo n.

Parameters:

a : integer

n : integer, n > 1. a and n should be relatively prime

Returns:

int : the order of a modulo n

Raises:

ValueError

If \(n \le 1\) or \(\gcd(a, n) \neq 1\). If a or n is not an integer.

Explanation

The order of a modulo n is the smallest integer k such that \(a^k\) leaves a remainder of 1 with n.

Examples

>>> from sympy.ntheory import n_order
>>> n_order(3, 7)
6
>>> n_order(4, 7)
3

See also

is_primitive_root

We say that a is a primitive root of n when the order of a modulo n equals totient(n)

sympy.ntheory.residue_ntheory.is_primitive_root(a, p)[source]

Returns True if a is a primitive root of p.

Parameters:

a : integer

p : integer, p > 1. a and p should be relatively prime

Returns:

bool : If True, a is the primitive root of p.

Raises:

ValueError

If \(p \le 1\) or \(\gcd(a, p) \neq 1\). If a or p is not an integer.

Explanation

a is said to be the primitive root of p if \(\gcd(a, p) = 1\) and \(\phi(p)\) is the smallest positive number s.t.

\(a^{\phi(p)} \equiv 1 \pmod{p}\).

where \(\phi(p)\) is Euler’s totient function.

The primitive root of p exist only for \(p = 2, 4, q^e, 2q^e\) (q is an odd prime). Hence, if it is not such a p, it returns False. To determine the primitive root, we need to know the prime factorization of q-1. The hardness of the determination depends on this complexity.

Examples

>>> from sympy.functions.combinatorial.numbers import totient
>>> from sympy.ntheory import is_primitive_root, n_order
>>> 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

See also

primitive_root

sympy.ntheory.residue_ntheory.primitive_root(p, smallest=True)[source]

Returns a primitive root of p or None.

Parameters:

p : integer, p > 1

smallest : if True the smallest primitive root is returned or None

Returns:

int | None :

If the primitive root exists, return the primitive root of p. If not, return None.

Raises:

ValueError

If \(p \le 1\) or p is not an integer.

Explanation

For the definition of primitive root, see the explanation of is_primitive_root.

The primitive root of p exist only for \(p = 2, 4, q^e, 2q^e\) (q is an odd prime). Now, if we know the primitive root of q, we can calculate the primitive root of \(q^e\), and if we know the primitive root of \(q^e\), we can calculate the primitive root of \(2q^e\). When there is no need to find the smallest primitive root, this property can be used to obtain a fast primitive root. On the other hand, when we want the smallest primitive root, we naively determine whether it is a primitive root or not.

Examples

>>> from sympy.ntheory.residue_ntheory import primitive_root
>>> primitive_root(19)
2
>>> primitive_root(21) is None
True
>>> primitive_root(50, smallest=False)
27

See also

is_primitive_root

References

[R701]

  1. Stein “Elementary Number Theory” (2011), page 44

[R702]

  1. Hackman “Elementary Number Theory” (2009), Chapter C

sympy.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 returned p // 2 only if it is the only root.

Use all_roots only when it is expected that all the roots fit in memory; otherwise use sqrt_mod_iter.

Examples

>>> from sympy.ntheory import sqrt_mod
>>> sqrt_mod(11, 43)
21
>>> sqrt_mod(17, 32, True)
[7, 9, 23, 25]
sympy.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 or Integer

Examples

>>> from sympy.ntheory.residue_ntheory import sqrt_mod_iter
>>> list(sqrt_mod_iter(11, 43))
[21, 22]

See also

sqrt_mod

Same functionality, but you want a sorted list or only one solution.

sympy.ntheory.residue_ntheory.quadratic_residues(p) → list[int][source]

Returns the list of quadratic residues.

Examples

>>> from sympy.ntheory.residue_ntheory import quadratic_residues
>>> quadratic_residues(7)
[0, 1, 2, 4]
sympy.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

Returns:

list[int] | int | None :

solutions to x**n = a mod p. The table of the output type is:

all_roots

has roots

Returns

True

Yes

list[int]

True

No

[]

False

Yes

int

False

No

None
Raises:

ValueError

If a, n or p is not integer. If n or p is not positive.

Examples

>>> from sympy.ntheory.residue_ntheory import nthroot_mod
>>> nthroot_mod(11, 4, 19)
8
>>> nthroot_mod(11, 4, 19, True)
[8, 11]
>>> nthroot_mod(68, 3, 109)
23

References

[R703]

  1. Hackman “Elementary Number Theory” (2009), page 76

sympy.ntheory.residue_ntheory.is_nthpow_residue(a, n, m)[source]

Returns True if x**n == a (mod m) has solutions.

References

[R704]

  1. Hackman “Elementary Number Theory” (2009), page 76

sympy.ntheory.residue_ntheory.is_quad_residue(a, p)[source]

Returns True if a (mod p) is in the set of squares mod p, i.e a % p in set([i**2 % p for i in range(p)]).

Parameters:

a : integer

p : positive integer

Returns:

bool : If True, x**2 == a (mod p) has solution.

Raises:

ValueError

If a, p is not integer. If p is not positive.

Examples

>>> from sympy.ntheory import is_quad_residue
>>> is_quad_residue(21, 100)
True

Indeed, pow(39, 2, 100) would be 21.

>>> is_quad_residue(21, 120)
False

That is, for any integer x, pow(x, 2, 120) is not 21.

If p is an odd prime, an iterative method is used to make the determination:

>>> from sympy.ntheory import is_quad_residue
>>> sorted(set([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

legendre_symbol, jacobi_symbol, sqrt_mod

sympy.ntheory.residue_ntheory.legendre_symbol(a, p)[source]

Returns the Legendre symbol \((a / p)\).

Deprecated since version 1.13: The legendre_symbol function is deprecated. Use sympy.functions.combinatorial.numbers.legendre_symbol instead. See its documentation for more information. See Relocate symbolic functions from ntheory to functions for details.

For an integer a and an odd prime p, 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

>>> from sympy.functions.combinatorial.numbers import legendre_symbol
>>> [legendre_symbol(i, 7) for i in range(7)]
[0, 1, 1, -1, 1, -1, -1]
>>> sorted(set([i**2 % 7 for i in range(7)]))
[0, 1, 2, 4]

See also

is_quad_residue, jacobi_symbol

sympy.ntheory.residue_ntheory.jacobi_symbol(m, n)[source]

Returns the Jacobi symbol \((m / n)\).

Deprecated since version 1.13: The jacobi_symbol function is deprecated. Use sympy.functions.combinatorial.numbers.jacobi_symbol instead. See its documentation for more information. See Relocate symbolic functions from ntheory to functions for details.

For any integer m and any positive odd integer n the Jacobi symbol is defined as the product of the Legendre symbols corresponding to the prime factors of n:

\[\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 modulo n.

But, unlike the Legendre symbol, if the Jacobi symbol \(\genfrac(){}{}{m}{n} = 1\) then m may or may not be a quadratic residue modulo n.

Parameters:

m : integer

n : odd positive integer

Examples

>>> from sympy.functions.combinatorial.numbers import jacobi_symbol, legendre_symbol
>>> from sympy import S
>>> jacobi_symbol(45, 77)
-1
>>> jacobi_symbol(60, 121)
1

The relationship between the jacobi_symbol and legendre_symbol can be demonstrated as follows:

>>> L = legendre_symbol
>>> S(45).factors()
{3: 2, 5: 1}
>>> jacobi_symbol(7, 45) == L(7, 3)**2 * L(7, 5)**1
True

See also

is_quad_residue, legendre_symbol

sympy.ntheory.residue_ntheory.mobius(n)[source]

Mobius function maps natural number to {-1, 0, 1}

Deprecated since version 1.13: The mobius function is deprecated. Use sympy.functions.combinatorial.numbers.mobius instead. See its documentation for more information. See Relocate symbolic functions from ntheory to functions for details.

It is defined as follows:

  1. \(1\) if \(n = 1\).

  2. \(0\) if \(n\) has a squared prime factor.

  3. \((-1)^k\) if \(n\) is a square-free 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 Mobius Function model).

Parameters:

n : positive integer

Examples

>>> from sympy.functions.combinatorial.numbers import mobius
>>> mobius(13*7)
1
>>> mobius(1)
1
>>> mobius(13*7*5)
-1
>>> mobius(13**2)
0

References

[R706]

Thomas Koshy “Elementary Number Theory with Applications”

sympy.ntheory.residue_ntheory.discrete_log(
n,
a,
b,
order=None,
prime_order=None,
)[source]

Compute the discrete logarithm of a to the base b modulo n.

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.

It employs different algorithms depending on the problem (subgroup order size, prime order or not):

  • Trial multiplication

  • Baby-step giant-step

  • Pollard’s Rho

  • Index Calculus

  • Pohlig-Hellman

Examples

>>> from sympy.ntheory import discrete_log
>>> discrete_log(41, 15, 7)
3

References

[R708]

“Handbook of applied cryptography”, Menezes, A. J., Van, O. P. C., & Vanstone, S. A. (1997).

sympy.ntheory.residue_ntheory.quadratic_congruence(a, b, c, n)[source]

Find the solutions to \(a x^2 + b x + c \equiv 0 \pmod{n}\).

Parameters:

a : int

b : int

c : int

n : int

A positive integer.

Returns:

list[int] :

A sorted list of solutions. If no solution exists, [].

Examples

>>> from sympy.ntheory.residue_ntheory import quadratic_congruence
>>> quadratic_congruence(2, 5, 3, 7) # 2x^2 + 5x + 3 = 0 (mod 7)
[2, 6]
>>> quadratic_congruence(8, 6, 4, 15) # No solution
[]

See also

polynomial_congruence

Solve the polynomial congruence

sympy.ntheory.residue_ntheory.polynomial_congruence(expr, m)[source]

Find the solutions to a polynomial congruence equation modulo m.

Parameters:

expr : integer coefficient polynomial

m : positive integer

Examples

>>> from sympy.ntheory import polynomial_congruence
>>> from sympy.abc import x
>>> expr = x**6 - 2*x**5 -35
>>> polynomial_congruence(expr, 6125)
[3257]

See also

sympy.polys.galoistools.gf_csolve

low level solving routine used by this routine

sympy.ntheory.residue_ntheory.binomial_mod(n, m, k)[source]

Compute binomial(n, m) % k.

Parameters:

n : an integer

m : an integer

k : a positive integer

Explanation

Returns binomial(n, m) % k using a generalization of Lucas’ Theorem for prime powers given by Granville [R709], in conjunction with the Chinese Remainder Theorem. The residue for each prime power is calculated in time O(log^2(n) + q^4*log(n)log(p) + q^4*p*log^3(p)).

Examples

>>> from sympy.ntheory.residue_ntheory import binomial_mod
>>> binomial_mod(10, 2, 6)  # binomial(10, 2) = 45
3
>>> binomial_mod(17, 9, 10)  # binomial(17, 9) = 24310
0

References

[R709] (1,2)

Binomial coefficients modulo prime powers, Andrew Granville, Available: https://web.archive.org/web/20170202003812/http://www.dms.umontreal.ca/~andrew/PDF/BinCoeff.pdf

sympy.ntheory.continued_fraction.continued_fraction(a) → list[source]

Return the continued fraction representation of a Rational or quadratic irrational.

Examples

>>> from sympy.ntheory.continued_fraction import continued_fraction
>>> from sympy import sqrt
>>> continued_fraction((1 + 2*sqrt(3))/5)
[0, 1, [8, 3, 34, 3]]

See also

continued_fraction_periodic, continued_fraction_reduce, continued_fraction_convergents

sympy.ntheory.continued_fraction.continued_fraction_convergents(cf)[source]

Return an iterator over the convergents of a continued fraction (cf).

The parameter should be in either of the following to forms: - A list of partial quotients, possibly with the last element being a list of repeating partial quotients, such as might be returned by continued_fraction and continued_fraction_periodic. - 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

>>> from sympy.core import pi
>>> from sympy import S
>>> from sympy.ntheory.continued_fraction import             continued_fraction_convergents, continued_fraction_iterator

>>> list(continued_fraction_convergents([0, 2, 1, 2]))
[0, 1/2, 1/3, 3/8]

>>> list(continued_fraction_convergents([1, S('1/2'), -7, S('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

>>> it = continued_fraction_convergents([1, [1, 2]])  # sqrt(3)
>>> for n in range(7):
...     print(next(it))
1
2
5/3
7/4
19/11
26/15
71/41

See also

continued_fraction_iterator, continued_fraction, continued_fraction_periodic

sympy.ntheory.continued_fraction.continued_fraction_iterator(x)[source]

Return continued fraction expansion of x as iterator.

Examples

>>> from sympy import Rational, pi
>>> from sympy.ntheory.continued_fraction import continued_fraction_iterator

>>> 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

sympy.ntheory.continued_fraction.continued_fraction_periodic(
p,
q,
d=0,
s=1,
) → list[source]

Find the periodic continued fraction expansion of a quadratic irrational.

Compute the continued fraction expansion of a rational or a quadratic irrational number, i.e. \(\frac{p + s\sqrt{d}}{q}\), where \(p\), \(q \ne 0\) and \(d \ge 0\) are integers.

Returns the continued fraction representation (canonical form) as a list of integers, optionally ending (for quadratic irrationals) with list of integers representing the repeating digits.

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

s : int, optional

the coefficient of the irrational part

Examples

>>> from sympy.ntheory.continued_fraction import continued_fraction_periodic
>>> 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]

See also

continued_fraction_iterator, continued_fraction_reduce

References

[R712]

K. Rosen. Elementary Number theory and its applications. Addison-Wesley, 3 Sub edition, pages 379-381, January 1992.

sympy.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

>>> from sympy.ntheory.continued_fraction import continued_fraction_reduce
>>> 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]).n(10)
2.718281835
>>> continued_fraction_reduce([1, 4, 2, [3, 1]])
(sqrt(21) + 287)/238
>>> continued_fraction_reduce([[1]])
(1 + sqrt(5))/2
>>> from sympy.ntheory.continued_fraction import continued_fraction_periodic
>>> continued_fraction_reduce(continued_fraction_periodic(8, 5, 13))
(sqrt(13) + 8)/5

See also

continued_fraction_periodic

sympy.ntheory.digits.count_digits(n, b=10)[source]

Return a dictionary whose keys are the digits of n in the given base, b, with keys indicating the digits appearing in the number and values indicating how many times that digit appeared.

Examples

>>> from sympy.ntheory import count_digits

>>> count_digits(1111339)
{1: 4, 3: 2, 9: 1}

The digits returned are always represented in base-10 but the number itself can be entered in any format that is understood by Python; the base of the number can also be given if it is different than 10:

>>> n = 0xFA; n
250
>>> count_digits(_)
{0: 1, 2: 1, 5: 1}
>>> count_digits(n, 16)
{10: 1, 15: 1}

The default dictionary will return a 0 for any digit that did not appear in the number. For example, which digits appear 7 times in 77!:

>>> from sympy import factorial
>>> c77 = count_digits(factorial(77))
>>> [i for i in range(10) if c77[i] == 7]
[1, 3, 7, 9]

See also

sympy.core.intfunc.num_digits, digits

sympy.ntheory.digits.digits(n, b=10, digits=None)[source]

Return a list of the digits of n in base b. The first element in the list is b (or -b if n is negative).

Parameters:

n: integer

The number whose digits are returned.

b: integer

The base in which digits are computed.

digits: integer (or None for all digits)

The number of digits to be returned (padded with zeros, if necessary).

Examples

>>> from sympy.ntheory.digits import digits
>>> digits(35)
[10, 3, 5]

If the number is negative, the negative sign will be placed on the base (which is the first element in the returned list):

>>> digits(-35)
[-10, 3, 5]

Bases other than 10 (and greater than 1) can be selected with b:

>>> digits(27, b=2)
[2, 1, 1, 0, 1, 1]

Use the digits keyword if a certain number of digits is desired:

>>> digits(35, digits=4)
[10, 0, 0, 3, 5]

See also

sympy.core.intfunc.num_digits, count_digits

sympy.ntheory.digits.is_palindromic(n, b=10)[source]

return True if n is the same when read from left to right or right to left in the given base, b.

Examples

>>> from sympy.ntheory import is_palindromic

>>> all(is_palindromic(i) for i in (-11, 1, 22, 121))
True

The second argument allows you to test numbers in other bases. For example, 88 is palindromic in base-10 but not in base-8:

>>> is_palindromic(88, 8)
False

On the other hand, a number can be palindromic in base-8 but not in base-10:

>>> 0o121, is_palindromic(0o121)
(81, False)

Or it might be palindromic in both bases:

>>> oct(121), is_palindromic(121, 8) and is_palindromic(121)
('0o171', True)
sympy.ntheory.egyptian_fraction.egyptian_fraction(r, algorithm='Greedy')[source]

Return the list of denominators of an Egyptian fraction expansion [R713] of the said rational \(r\).

Parameters:

r : Rational or (p, q)

a positive rational number, p/q.

algorithm : { “Greedy”, “Graham Jewett”, “Takenouchi”, “Golomb” }, optional

Denotes the algorithm to be used (the default is “Greedy”).

Examples

>>> from sympy import Rational
>>> from sympy.ntheory.egyptian_fraction import egyptian_fraction
>>> egyptian_fraction(Rational(3, 7))
[3, 11, 231]
>>> egyptian_fraction((3, 7), "Graham Jewett")
[7, 8, 9, 56, 57, 72, 3192]
>>> egyptian_fraction((3, 7), "Takenouchi")
[4, 7, 28]
>>> egyptian_fraction((3, 7), "Golomb")
[3, 15, 35]
>>> egyptian_fraction((11, 5), "Golomb")
[1, 2, 3, 4, 9, 234, 1118, 2580]

Notes

Currently the following algorithms are supported:

  1. Greedy Algorithm

    Also called the Fibonacci-Sylvester algorithm [R714]. At each step, extract the largest unit fraction less than the target and replace the target with the remainder.

    It has some distinct properties:

    1. 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.

    2. Uses minimal memory.

    3. The terms can blow up (standard examples of this are 5/121 and 31/311). The denominator is at most squared at each step (doubly-exponential growth) and typically exhibits singly-exponential growth.

  2. 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. See [R715].

  3. Takenouchi Algorithm

    The algorithm suggested by Takenouchi (1921). Differs from the Graham-Jewett algorithm only in the handling of duplicates. See [R715].

  4. 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). See [R716].

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.

See also

sympy.core.numbers.Rational

References

sympy.ntheory.bbp_pi.pi_hex_digits(n, prec=14)[source]

Returns a string containing prec (default 14) digits starting at the nth digit of pi in hex. Counting of digits starts at 0 and the decimal is not counted, so for n = 0 the returned value starts with 3; n = 1 corresponds to the first digit past the decimal point (which in hex is 2).

Parameters:

n : non-negative integer

prec : non-negative integer. default = 14

Returns:

str : Returns a string containing prec digits

starting at the nth digit of pi in hex. If prec = 0, returns empty string.

Raises:

ValueError

If n < 0 or prec < 0. Or n or prec is not an integer.

Examples

>>> from sympy.ntheory.bbp_pi import pi_hex_digits
>>> pi_hex_digits(0)
'3243f6a8885a30'
>>> pi_hex_digits(0, 3)
'324'

These are consistent with the following results

>>> import math
>>> hex(int(math.pi * 2**((14-1)*4)))
'0x3243f6a8885a30'
>>> hex(int(math.pi * 2**((3-1)*4)))
'0x324'

References

ECM function

The \(ecm\) function is a subexponential factoring algorithm capable of factoring numbers of around ~35 digits comfortably within few seconds. The time complexity of \(ecm\) is dependent on the smallest proper factor of the number. So even if the number is really large but its factors are comparatively smaller then \(ecm\) can easily factor them. For example we take \(N\) with 15 digit factors \(15154262241479\), \(15423094826093\), \(799333555511111\), \(809709509409109\), \(888888877777777\), \(914148152112161\). Now N is a 87 digit number. \(ECM\) takes under around 47s to factorise this.

sympy.ntheory.ecm.ecm(
n,
B1=10000,
B2=100000,
max_curve=200,
seed=1234,
)[source]

Performs factorization using Lenstra’s Elliptic curve method.

This function repeatedly calls _ecm_one_factor to compute the factors of n. First all the small factors are taken out using trial division. Then _ecm_one_factor is used to compute one factor at a time.

Parameters:

n : Number to be Factored

B1 : Stage 1 Bound. Must be an even number.

B2 : Stage 2 Bound. Must be an even number.

max_curve : Maximum number of curves generated

seed : Initialize pseudorandom generator

Examples

>>> from sympy.ntheory import ecm
>>> ecm(25645121643901801)
{5394769, 4753701529}
>>> ecm(9804659461513846513)
{4641991, 2112166839943}

Examples

>>> from sympy.ntheory import ecm
>>> ecm(7060005655815754299976961394452809, B1=100000, B2=1000000)
{6988699669998001, 1010203040506070809}
>>> ecm(122921448543883967430908091422761898618349713604256384403202282756086473494959648313841, B1=100000, B2=1000000)
{15154262241479,
15423094826093,
799333555511111,
809709509409109,
888888877777777,
914148152112161}

QS function

The \(qs\) function is a subexponential factoring algorithm, the fastest factoring algorithm for numbers within 100 digits. The time complexity of \(qs\) is dependent on the size of the number so it is used if the number contains large factors. Due to this while factoring numbers first \(ecm\) is used to get smaller factors of around ~15 digits then \(qs\) is used to get larger factors.

For factoring \(2709077133180915240135586837960864768806330782747\) which is a semi-prime number with two 25 digit factors. \(qs\) is able to factorize this in around 248s.

sympy.ntheory.qs.qs(N, prime_bound, M, ERROR_TERM=25, seed=1234)[source]

Performs factorization using Self-Initializing Quadratic Sieve. In SIQS, let N be a number to be factored, and this N should not be a perfect power. If we find two integers such that X**2 = Y**2 modN and X != +-Y modN, then \(gcd(X + Y, N)\) will reveal a proper factor of N. In order to find these integers X and Y we try to find relations of form t**2 = u modN where u is a product of small primes. If we have enough of these relations then we can form (t1*t2...ti)**2 = u1*u2...ui modN such that the right hand side is a square, thus we found a relation of X**2 = Y**2 modN.

Here, several optimizations are done like using multiple polynomials for sieving, fast changing between polynomials and using partial relations. The use of partial relations can speeds up the factoring by 2 times.

Parameters:

N : Number to be Factored

prime_bound : upper bound for primes in the factor base

M : Sieve Interval

ERROR_TERM : Error term for checking smoothness

threshold : Extra smooth relations for factorization

seed : generate pseudo prime numbers

Examples

>>> from sympy.ntheory import qs
>>> qs(25645121643901801, 2000, 10000)
{5394769, 4753701529}
>>> qs(9804659461513846513, 2000, 10000)
{4641991, 2112166839943}

References

Examples

>>> from sympy.ntheory import qs
>>> qs(5915587277*3267000013, 1000, 10000)
{3267000013, 5915587277}