## Number Theory with Sympy's Sieve

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.


```python
from sympy import sieve
```


```python
sieve._reset()
25 in sieve
```




    False




```python
sieve._list
```




    array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23])




```python
# Grow the sieve to cover all primes <= n
sieve._reset()
sieve.extend(30)
sieve[10] == 28
```




    False




```python
sieve[10] == 29
```




    True




```python
sieve[10] == 23
```




    False




```python
# Extend to include the ith prime number
sieve._reset()
sieve.extend_to_no(9)
sieve._list
```




    array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23])




```python
# $primerange(a,b)$

print([i for i in sieve.primerange(7, 23)])
```

    [7, 11, 13, 17, 19]



```python
# Search = returns the indice i, j of the primes that bound n
#if n is prime then i = j

sieve.search(25)
```




    (9, 10)




```python
sieve.search(23)
```




    (9, 9)




```python
# Prime
# Return the nth prime, with the primes indexed as prime(1) = 2, prime(2) = 3, etc…. 
# The nth prime is approximately n*log(n).
from sympy import prime
prime(10)
```




    29



## Primes 


```python
prime(1)
```




    2




```python
%time
prime(1000000)
```

    CPU times: user 3 µs, sys: 0 ns, total: 3 µs
    Wall time: 6.91 µs





    15485863




```python
# primepi(n) - gives n number of primes

from sympy import primepi
primepi(25)
```




    9




```python
%time
primepi(1000000)
```

    CPU times: user 3 µs, sys: 0 ns, total: 3 µs
    Wall time: 6.91 µs





    78498




```python
from sympy import nextprime
[(i, nextprime(i)) for i in range(10, 15)]
```




    [(10, 11), (11, 13), (12, 13), (13, 17), (14, 17)]




```python
from sympy import prevprime
[(i, prevprime(i)) for i in range(10, 15)]
```




    [(10, 7), (11, 7), (12, 11), (13, 11), (14, 13)]



## Prime Ranges

Some famous conjectures about the occurence 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.


```python
from sympy import primerange, sieve
print([i for i in primerange(1, 30)])
```

    [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]



```python
list(sieve.primerange(1, 30))
```




    [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]




```python
# randprime
from sympy import randprime, isprime
randprime(1, 30)
```




    19




```python
isprime(randprime(1, 30))
```




    True




```python
# This returns the product of the first n primes or teh primes <= n (when nth=False)
from sympy.ntheory.generate import primorial, randprime, primerange
from sympy import factorint, Mul, primefactors, sqrt
primorial(5) # product of 2, 3, 5, 7, 11
```




    2310




```python
2*3*5*7*11
```




    2310




```python
primorial(2)
```




    6




```python
primorial(3, nth=False) # primes <= 3 are 2 and 3
```




    6




```python
primorial(5, nth=False) # product of 2*3*5
```




    30




```python
primorial(sqrt(100), nth=False)
```




    210




```python
# Adding or subtracting by 1 of a primorial product gives you a prime
factorint(primorial(5) - 1)
```




    {2309: 1}




```python
# here we get two new primes that are factors
factorint(primorial(7) - 1)
```




    {61: 1, 8369: 1}




```python
# Some primes smaller and larger than the primes multiplied together
p = list(primerange(10, 20))
sorted(set(primefactors(Mul(*p) + 1)).difference(set(p)))
```




    [2, 5, 31, 149]



### cycle_length

$cycle_length(f, x0, nmax=None, values=False)$

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.


```python
from sympy.ntheory.generate import cycle_length # will give succesive values of i <- func(i)

```


```python
def iter(func, i):
    while 1:
        ii = func(1)
        yield ii
        i = ii
# give a seed of 4 and the mu and lambda terms
func = lambda i: (i**2 + 1) % 51
```


```python
next(cycle_length(func, 4))
```




    (6, 2)




```python
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]



### composite(nth)

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


```python
from sympy import composite
composite(24)
```




    36




```python
composite(1)
```




    4



### compositepi


```python
from sympy import compositepi
compositepi(20)
```




    11



### smoothness_p(n, m=-1, power=0, visual=None)

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.


```python
from sympy.ntheory.factor_ import smoothness_p, factorint
smoothness_p(10345, m=1)
```




    (1, [(5, (1, 3, 3)), (2069, (1, 23, 23))])




```python
smoothness_p(10345)
```




    (-1, [(5, (1, 2, 4)), (2069, (1, 47, 47))])




```python
smoothness_p(10345, power=1)
```




    (-1, [(5, (1, 2, 4)), (2069, (1, 47, 47))])




```python
print(smoothness_p(344556576677878, visual=1))
```

    p**i=2**1 has p-1 B=1, B-pow=1
    p**i=172278288338939**1 has p-1 B=18836462753, B-pow=18836462753



```python
factorint(15*11)
```




    {3: 1, 5: 1, 11: 1}




```python
smoothness_p(_)
```




    'p**i=3**1 has p-1 B=2, B-pow=2\np**i=5**1 has p-1 B=2, B-pow=4\np**i=11**1 has p-1 B=5, B-pow=5'




```python
smoothness_p(_)
```




    {3: 1, 5: 1, 11: 1}



### Table for output logic is like this

#### Visual

| Input         | True        | False       | Other       |
|---------------|-------------|-------------|-------------|
| ``str``      |``str``      |``tuple``    |``str``      |
| ``str``      |``str``      |``tuple``    |``dict``      |
| ``tuple``    |``str``      |``tuple``    |``str``      |
| ``n``        |``str``      |``tuple``    |``tuple``      |
| ``mul``      |``str``      |``tuple``    |``tuple``      |


```python
# training(n)
# Count the number of trailing zero digits in the binary representation of n, 
# i.e. determine the largest power of 2 that divides n.
from sympy import trailing
trailing(128)
```




    7




```python
trailing(51)
```




    0




```python
# multiplicity
# Find the greatest integer m such that p**m divides n.
from sympy.ntheory import multiplicity
from sympy.core.numbers import Rational as R
[multiplicity(5, n) for n in [8, 5, 25, 125, 250]]
```




    [0, 1, 2, 3, 3]




```python
multiplicity(3, R(1, 9))
```




    -2



###  sympy.ntheory.factor_.perfect_power

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

Return `(b, e)` such that `n == b**e` if n is a perfect power; otherwise return False.

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




```python
from sympy import perfect_power
perfect_power(16)
```




    (2, 4)




```python
perfect_power(25, big=False)
```




    (5, 2)



### Pollard_rho

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:


```python
from sympy.ntheory.generate import cycle_length
n = 14345656
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 =  660; leader length =  19
    loop length =  660; leader length =  19
    loop length =  120; leader length =  26
    loop length =  660; leader length =   6
    loop length =  660; leader length =  17



```python
# An explicit example where there is a two element leadup to a seq of 3 numbers
x = 2
for i in range(9):
    x = (x**2 + 12)%17
    print(x)
```

    16
    13
    11
    14
    4
    11
    14
    4
    11



```python
next(cycle_length(lambda x:(x**2+12)%17, 2))
```




    (3, 2)




```python
list(cycle_length(lambda x: (x**2+12)%17, 2, values=True))
```




    [16, 13, 11, 14, 4]



### Note

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.


```python
from sympy import pollard_rho
n = 14345656
F = lambda x:(2048*pow(x,2,n) + 32767) % n
pollard_rho(n, F=F)
```




    8




```python
pollard_rho(n, a=n-2, retries=1)
```




    8




```python

```

943b80ad53e4270fe76872e0e554f52cf135600f
ipynb
Jupyter Notebook
sympy/Number Theory.ipynb
DanielMabadeje/Artificial-Intelligence-Deep-Learning-Machine-Learning-Tutorials
7adab3877fc1d3f1d5f57e6c1743dae8f76f72c5
Apache-2.0
2017-08-06T16:51:46.000Z
2022-03-30T07:34:24.000Z
sympy/Number Theory.ipynb
nuhaltinsoy/Artificial-Intelligence-Deep-Learning-Machine-Learning-Tutorials
6017441f2d476f9c6c568dd886da43c6c0fd89bd
Apache-2.0
2017-08-28T14:59:36.000Z
2022-03-11T23:21:35.000Z
sympy/Number Theory.ipynb
nuhaltinsoy/Artificial-Intelligence-Deep-Learning-Machine-Learning-Tutorials
6017441f2d476f9c6c568dd886da43c6c0fd89bd
Apache-2.0
2017-08-06T17:40:59.000Z
2022-03-31T12:03:24.000Z
Qwen/Qwen-72B
1. YES
2. YES