CoCalc Shared Filessphinx2ipynb / doc / build / ipynb / sympy / doc / src / modules / ntheory.ipynbOpen in CoCalc with one click!

1

2

class sympy.ntheory.generate.Sieve

3

4

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.

5

In [ ]:

>>> from sympy import sieve >>> from array import array # this line and next for doctest only >>> sieve._list = array('l', [2, 3, 5, 7, 11, 13])

6

In [ ]:

>>> 25 in sieve False >>> sieve._list array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23])

7

extend(n)

8

9

Grow the sieve to cover all primes <= n (a real number).

10

In [ ]:

>>> from sympy import sieve >>> from array import array # this line and next for doctest only >>> sieve._list = array('l', [2, 3, 5, 7, 11, 13])

11

In [ ]:

>>> sieve.extend(30) >>> sieve[10] == 29 True

12

extend_to_no(i)

13

14

Extend to include the ith prime number.

15

i must be an integer.

16

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

17

In [ ]:

>>> from sympy import sieve >>> from array import array # this line and next for doctest only >>> sieve._list = array('l', [2, 3, 5, 7, 11, 13])

18

In [ ]:

>>> sieve.extend_to_no(9) >>> sieve._list array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23])

19

primerange(a, b)

20

21

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

22

In [ ]:

>>> from sympy import sieve >>> print([i for i in sieve.primerange(7, 18)]) [7, 11, 13, 17]

23

search(n)

24

25

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

26

If n is prime then i == j.

27

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

28

In [ ]:

>>> from sympy import sieve >>> sieve.search(25) (9, 10) >>> sieve.search(23) (9, 9)

29

30

sympy.ntheory.generate.prime(nth)

31

32

Return the nth prime, with the primes indexed as prime(1) = 2, prime(2) = 3, etc.... The nth prime is approximately n*log(n) and can never be larger than 2**n.

33

34

In [ ]:

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

36

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

37

sympy.ntheory.generate.primepi(n)

38

39

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

40

In [ ]:

>>> from sympy import primepi >>> primepi(25) 9

41

sympy.ntheory.primetest.isprime : Test if n is prime primerange : Generate all primes in a given range prime : Return the nth prime

42

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

43

44

Return the ith prime greater than n.

45

i must be an integer.

46

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

47

In [ ]:

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

48

prevprime : Return the largest prime smaller than n primerange : Generate all primes in a given range

49

sympy.ntheory.generate.prevprime(n)

50

51

Return the largest prime smaller than n.

52

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

53

In [ ]:

>>> from sympy import prevprime >>> [(i, prevprime(i)) for i in range(10, 15)] [(10, 7), (11, 7), (12, 11), (13, 11), (14, 13)]

54

nextprime : Return the ith prime greater than n primerange : Generates all primes in a given range

55

sympy.ntheory.generate.primerange(a, b)

56

57

Generate a list of all prime numbers in the range [a, b).

58

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.

59

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

60

61

Twin primes: though often not, the following will give 2 primes

62

63

an infinite number of times:

64

65

primerange(6*n - 1, 6*n + 2)

66

67

Legendre's: the following always yields at least one prime

68

69

primerange(n**2, (n+1)**2+1)

70

71

Bertrand's (proven): there is always a prime in the range

72

73

primerange(n, 2*n)

74

75

Brocard's: there are at least four primes in the range

76

77

primerange(prime(n)**2, prime(n+1)**2)

78

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.

79

80

82

In [ ]:

>>> from sympy import primerange, sieve >>> print([i for i in primerange(1, 30)]) [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

84

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:

85

In [ ]:

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

86

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

87

88

or extend the sieve to contain the requested range.

89

sympy.ntheory.generate.randprime(a, b)

90

91

Return a random prime number in the range [a, b).

92

Bertrand's postulate assures that randprime(a, 2*a) will always succeed for a > 1.

93

94

In [ ]:

>>> from sympy import randprime, isprime >>> randprime(1, 30) #doctest: +SKIP 13 >>> isprime(randprime(1, 30)) True

96

primerange : Generate all primes in a given range

97

sympy.ntheory.generate.primorial(n, nth=True)

98

99

Returns the product of the first n primes (default) or the primes less than or equal to n (when "nth=False").

100

In [ ]:

>>> from sympy.ntheory.generate import primorial, randprime, 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

101

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.

102

In this case, the number itself is a new prime:

103

In [ ]:

>>> factorint(primorial(4) + 1) {211: 1}

104

In this case two new primes are the factors:

105

In [ ]:

>>> factorint(primorial(4) - 1) {11: 1, 19: 1}

106

Here, some primes smaller and larger than the primes multiplied together are obtained:

107

In [ ]:

>>> p = list(primerange(10, 20)) >>> sorted(set(primefactors(Mul(*p) + 1)).difference(set(p))) [2, 5, 31, 149]

108

primerange : Generate all primes in a given range

109

sympy.ntheory.generate.cycle_length(f, x0, nmax=None, values=False)

110

111

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

112

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.

113

In [ ]:

>>> from sympy.ntheory.generate import cycle_length

114

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

115

In [ ]:

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

116

A function is defined:

117

In [ ]:

>>> func = lambda i: (i**2 + 1) % 51

118

and given a seed of 4 and the mu and lambda terms calculated:

119

In [ ]:

>>> next(cycle_length(func, 4)) (6, 2)

120

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

121

In [ ]:

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

122

There are 6 repeating values after the first 2.

123

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

124

In [ ]:

>>> next(cycle_length(func, 4, nmax = 4)) (4, None) >>> [ni for ni in cycle_length(func, 4, nmax = 4, values=True)] [17, 35, 2, 5]

125

Code modified from:

126

127

sympy.ntheory.factor_.smoothness(n)

129

130

Return the B-smooth and B-power smooth values of n.

131

The smoothness of n is the largest prime factor of n; the power- smoothness is the largest divisor raised to its multiplicity.

132

In [ ]:

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

133

factorint, smoothness_p

134

sympy.ntheory.factor_.smoothness_p(n, m=-1, power=0, visual=None)

135

136

Return a list of [m, (p, (M, sm(p + m), psm(p + m)))...] where:

137

138

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

139

140

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

141

142

psm(p + m) is the power smoothness of p + m

143

The list is sorted according to smoothness (default) or by power smoothness if power=1.

144

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.

145

In [ ]:

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

146

If visual=True then an annotated string will be returned:

147

In [ ]:

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

148

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

149

In [ ]:

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

150

The table of the output logic is:

151

152

153

In [ ]:

154

Visual

155

Input

156

True

157

False

158

other

159

dict

160

str

161

tuple

162

str

163

str

164

str

165

tuple

166

dict

167

tuple

168

str

169

tuple

170

str

171

n

172

str

173

tuple

174

tuple

175

mul

176

str

177

tuple

178

tuple

179

factorint, smoothness

180

sympy.ntheory.factor_.trailing(n)

181

182

Count the number of trailing zero digits in the binary representation of n, i.e. determine the largest power of 2 that divides n.

183

In [ ]:

>>> from sympy import trailing >>> trailing(128) 7 >>> trailing(63) 0

184

sympy.ntheory.factor_.multiplicity(p, n)

185

186

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

187

In [ ]:

>>> 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] >>> multiplicity(3, R(1, 9)) -2

188

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

189

190

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

191

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.

192

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.

193

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.

194

In [ ]:

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

195

sympy.ntheory.factor_.pollard_rho(n, s=2, a=1, retries=5, seed=1234, max_steps=None, F=None)

196

197

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.

198

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.

199

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

200

For a given function, very different leader-loop values can be obtained so it is a good idea to allow for retries:

201

In [ ]:

>>> 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 = 42 loop length = 78; leader length = 120 loop length = 1482; leader length = 99 loop length = 1482; leader length = 285 loop length = 1482; leader length = 100

202

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

203

In [ ]:

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

204

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.

205

In [ ]:

>>> from sympy import pollard_rho >>> n=16843009 >>> F=lambda x:(2048*pow(x,2,n) + 32767) % n >>> pollard_rho(n, F=F) 257

206

Use the default setting with a bad value of "a" and no retries:

207

In [ ]:

>>> pollard_rho(n, a=n-2, retries=0)

208

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

209

In [ ]:

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

210

211

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

212

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

213

214

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

215

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.

216

Note: the value of M is lcm(1..B) = reduce(ilcm, range(2, B + 1)).

217

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 mumber just less than the factor p (hence the name p - 1).

218

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:

219

In [ ]:

>>> from sympy.ntheory.factor_ import smoothness_p, pollard_pm1 >>> n=257*1009 >>> smoothness_p(n) (-1, [(257, (1, 2, 256)), (1009, (1, 7, 16))])

220

So we should (and can) find a root with B=16:

221

In [ ]:

>>> pollard_pm1(n, B=16, a=3) 1009

222

If we attempt to increase B to 256 we find that it doesn't work:

223

In [ ]:

>>> pollard_pm1(n, B=256) >>>

224

But if the value of "a" is changed we find that only multiples of 257 work, e.g.:

225

In [ ]:

>>> pollard_pm1(n, B=256, a=257) 1009

226

Checking different "a" values shows that all the ones that didn't work had a gcd value not equal to "n" but equal to one of the factors:

227

In [ ]:

>>> from sympy.core.numbers import ilcm, igcd >>> from sympy import 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]) set([1009])

228

But does aM % d for every divisor of n give 1?

229

In [ ]:

>>> aM = pow(255, M, n) >>> [(d, aM%Pow(*d.args)) for d in factorint(n, visual=True).args] [(257**1, 1), (1009**1, 1)]

230

No, only one of them. So perhaps the principle is that a root will be found for a given value of B provided that:

231

232

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

233

234

a**M % p != 1 for any of the divisors of n.

235

By trying more than one "a" it is possible that one of them will yield a factor.

236

With the default smoothness bound, this number can't be cracked:

237

In [ ]:

>>> from sympy.ntheory import pollard_pm1, primefactors >>> pollard_pm1(21477639576571)

238

Increasing the smoothness bound helps:

239

In [ ]:

>>> pollard_pm1(21477639576571, B=2000) 4410317

240

Looking at the smoothness of the factors of this number we find:

241

In [ ]:

>>> from sympy.utilities import flatten >>> 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

242

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:

243

In [ ]:

>>> factorint(4410317 - 1) {2: 2, 617: 1, 1787: 1} >>> factorint(4869863-1) {2: 1, 2434931: 1}

244

Note that until B reaches the B-pow value of 1787, the number is not cracked;

245

In [ ]:

>>> pollard_pm1(21477639576571, B=1786) >>> pollard_pm1(21477639576571, B=1787) 4410317

246

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 don't
know. The modular.math reference below states that 15% of numbers in the
range of 10**15 to 15**15 + 10**4 are 10**6 power smooth so a B of 10**6
will fail 85% of the time in that range. From 10**8 to 10**8 + 10**3 the
percentages are nearly reversed...but in that range the simple trial
division is quite fast.

247

248

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

249

250

252

sympy.ntheory.factor_.factorint(n, limit=None, use_trial=True, use_rho=True, use_pm1=True, verbose=False, visual=None)

254

255

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:

256

In [ ]:

>>> from sympy.ntheory import factorint >>> factorint(2000) # 2000 = (2**4) * (5**3) {2: 4, 5: 3} >>> factorint(65537) # This number is prime {65537: 1}

257

For input less than 2, factorint behaves as follows:

258

259

260

"factorint(1)" returns the empty factorization, "{}"

261

262

"factorint(0)" returns "{0:1}"

263

264

"factorint(-n)" adds "-1:1" to the factors and then factors "n"

265

Partial Factorization:

266

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.

267

This number, for example, has two small factors and a huge semi-prime factor that cannot be reduced easily:

268

In [ ]:

>>> from sympy.ntheory import isprime >>> from sympy.core.compatibility import long >>> a = 1407633717262338957430697921446883 >>> f = factorint(a, limit=10000) >>> f == {991: 1, long(202916782076162456022877024859): 1, 7: 1} True >>> isprime(max(f)) False

269

This number has a small factor and a residual perfect power whose base is greater than the limit:

270

In [ ]:

>>> factorint(3*101**7, limit=5) {3: 1, 101: 7}

271

Visual Factorization:

272

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

273

In [ ]:

>>> from sympy import pprint >>> pprint(factorint(4200, visual=True)) 3 1 2 1 2 *3 *5 *7

274

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.

275

You can easily switch between the two forms by sending them back to factorint:

276

In [ ]:

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

277

In [ ]:

>>> visual = factorint(1764, visual=True); pprint(visual) 2 2 2 2 *3 *7 >>> print(factorint(visual)) {2: 2, 3: 2, 7: 2}

278

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:

279

In [ ]:

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

280

The table of the output logic is:

281

282

283

Input

284

True

285

False

286

other

287

dict

288

mul

289

dict

290

mul

291

n

292

mul

293

dict

294

dict

295

mul

296

mul

297

dict

298

dict

299

Algorithm:

300

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:

301

In [ ]:

>>> factors = factorint(12345678910111213141516) >>> for base, exp in sorted(factors.items()): ... print('%s %s' % (base, exp)) ... 2 2 2507191691 1 1231026625769 1

302

Any of these methods can optionally be disabled with the following boolean parameters:

303

304

305

"use_trial": Toggle use of trial division

306

307

"use_rho": Toggle use of Pollard's rho method

308

309

"use_pm1": Toggle use of Pollard's p-1 method

310

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

311

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

312

smoothness, smoothness_p, divisors

313

sympy.ntheory.factor_.primefactors(n, limit=None, verbose=False)

314

315

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.

316

In [ ]:

>>> from sympy.ntheory import primefactors, factorint, isprime >>> primefactors(6) [2, 3] >>> primefactors(-5) [5]

317

In [ ]:

>>> sorted(factorint(123456).items()) [(2, 6), (3, 1), (643, 1)] >>> primefactors(123456) [2, 3, 643]

318

In [ ]:

>>> sorted(factorint(10000000001, limit=200).items()) [(101, 1), (99009901, 1)] >>> isprime(99009901) False >>> primefactors(10000000001, limit=300) [101]

319

divisors

320

sympy.ntheory.factor_.divisors(n, generator=False)

321

322

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

323

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

324

In [ ]:

>>> from sympy import divisors, divisor_count >>> divisors(24) [1, 2, 3, 4, 6, 8, 12, 24] >>> divisor_count(24) 8

325

In [ ]:

>>> list(divisors(120, generator=True)) [1, 2, 4, 8, 3, 6, 12, 24, 5, 10, 20, 40, 15, 30, 60, 120]

326

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

327

primefactors, factorint, divisor_count

328

sympy.ntheory.factor_.divisor_count(n, modulus=1)

329

330

Return the number of divisors of "n". If "modulus" is not 1 then only those that are divisible by "modulus" are counted.

331

332

In [ ]:

>>> from sympy import divisor_count >>> divisor_count(6) 4

334

factorint, divisors, totient

335

sympy.ntheory.factor_.totient(*args, **kw_args)

336

337

Calculate the Euler totient function phi(n)

338

In [ ]:

>>> from sympy.ntheory import totient >>> totient(1) 1 >>> totient(25) 20

339

divisor_count

340

sympy.ntheory.modular.symmetric_residue(a, m)

341

342

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

343

In [ ]:

>>> from sympy.ntheory.modular import symmetric_residue >>> symmetric_residue(1, 6) 1 >>> symmetric_residue(4, 6) -2

344

sympy.ntheory.modular.crt(m, v, symmetric=False, check=True)

345

346

Chinese Remainder Theorem.

347

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.

348

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.

349

The keyword "check" can be set to False if it is known that the moduli are coprime.

350

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

351

In [ ]:

>>> from sympy.ntheory.modular import crt, solve_congruence >>> crt([99, 97, 95], [49, 76, 65]) (639985, 912285)

352

This is the correct result because:

353

In [ ]:

>>> [639985 % m for m in [99, 97, 95]] [49, 76, 65]

354

If the moduli are not co-prime, you may receive an incorrect result if you use "check=False":

355

In [ ]:

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

356

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

357

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.

358

solve_congruence sympy.polys.galoistools.gf_crt : low level crt routine used by this routine

359

sympy.ntheory.modular.crt1(m)

360

361

First part of Chinese Remainder Theorem, for multiple application.

362

In [ ]:

>>> from sympy.ntheory.modular import crt1 >>> crt1([18, 42, 6]) (4536, [252, 108, 756], [0, 2, 0])

363

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

364

365

Second part of Chinese Remainder Theorem, for multiple application.

366

In [ ]:

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

367

sympy.ntheory.modular.solve_congruence(*remainder_modulus_pairs, **hint)

368

369

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.

370

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.

371

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

372

In [ ]:

>>> from sympy.ntheory.modular import solve_congruence

373

What number is 2 mod 3, 3 mod 5 and 2 mod 7?

374

In [ ]:

>>> solve_congruence((2, 3), (3, 5), (2, 7)) (23, 105) >>> [23 % m for m in [3, 5, 7]] [2, 3, 2]

375

If you prefer to work with all remainder in one list and all moduli in another, send the arguments like this:

376

In [ ]:

>>> solve_congruence(*zip((2, 3, 2), (3, 5, 7))) (23, 105)

377

The moduli need not be co-prime; in this case there may or may not be a solution:

378

In [ ]:

>>> solve_congruence((2, 3), (4, 6)) is None True

379

In [ ]:

>>> solve_congruence((2, 3), (5, 6)) (5, 6)

380

The symmetric flag will make the result be within 1/2 of the modulus:

381

In [ ]:

>>> solve_congruence((2, 3), (5, 6), symmetric=True) (-1, 6)

382

crt : high level routine implementing the Chinese Remainder Theorem

383

sympy.ntheory.multinomial.binomial_coefficients(n)

384

385

386

In [ ]:

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

387

binomial_coefficients_list, multinomial_coefficients

388

sympy.ntheory.multinomial.binomial_coefficients_list(n)

389

390

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

391

In [ ]:

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

392

binomial_coefficients, multinomial_coefficients

393

sympy.ntheory.multinomial.multinomial_coefficients(m, n)

394

395

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

396

For example:

397

In [ ]:

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

398

The algorithm is based on the following result:

399

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

400

binomial_coefficients_list, binomial_coefficients

401

sympy.ntheory.multinomial.multinomial_coefficients_iterator(m, n, _tuple=<type 'tuple'>)

402

403

multinomial coefficient iterator

404

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.

405

In [ ]:

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

406

In [ ]:

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

407

sympy.ntheory.partitions_.npartitions(n, verbose=False)

408

409

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

410

P(n) is computed using the Hardy-Ramanujan-Rademacher formula, described e.g. at http://mathworld.wolfram.com/PartitionFunctionP.html

411

The correctness of this implementation has been tested for 10**n up to n = 8.

412

In [ ]:

>>> from sympy.ntheory import npartitions >>> npartitions(25) 1958

413

sympy.ntheory.primetest.mr(n, bases)

414

415

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

416

417

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

418

A list of thresholds and the bases they require are here: http://en.wikipedia.org/wiki/Millerâ€“Rabin_primality_test#Deterministic_variants_of_the_test

419

In [ ]:

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

420

sympy.ntheory.primetest.isprime(n)

421

422

Test if n is a prime number (True) or not (False). For n < 10**16 the answer is accurate; greater n values have a small probability of actually being pseudoprimes.

423

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

424

The function first looks for trivial factors, and if none is found, performs a safe Miller-Rabin strong pseudoprime test with bases that are known to prove a number prime. Finally, a general Miller-Rabin test is done with the first k bases which will report a pseudoprime as a prime with an error of about 4**-k. The current value of k is 46 so the error is about 2 x 10**-28.

425

In [ ]:

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

426

sympy.ntheory.generate.primerange : Generates all primes in a given range sympy.ntheory.generate.primepi : Return the number of primes less than or equal to n sympy.ntheory.generate.prime : Return the nth prime

427

sympy.ntheory.residue_ntheory.n_order(a, n)

428

429

Returns the order of "a" modulo "n".

430

The order of "a" modulo "n" is the smallest integer "k" such that "a**k" leaves a remainder of 1 with "n".

431

In [ ]:

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

432

sympy.ntheory.residue_ntheory.is_primitive_root(a, p)

433

434

Returns True if "a" is a primitive root of "p"

435

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

436

437

a**totient(p) cong 1 mod(p)

438

In [ ]:

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

439

sympy.ntheory.residue_ntheory.primitive_root(p)

440

441

Returns the smallest primitive root or None

442

[1] W. Stein "Elementary Number Theory" (2011), page 44 [2] P. Hackman "Elementary Number Theory" (2009), Chapter C

443

p : positive integer

444

In [ ]:

>>> from sympy.ntheory.residue_ntheory import primitive_root >>> primitive_root(19) 2

445

sympy.ntheory.residue_ntheory.sqrt_mod(a, p, all_roots=False)

446

447

find a root of "x**2 = a mod p"

448

a : integer p : positive integer all_roots : if True the list of roots is returned or None

449

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.

450

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

451

In [ ]:

>>> from sympy.ntheory import sqrt_mod >>> sqrt_mod(11, 43) 21 >>> sqrt_mod(17, 32, True) [7, 9, 23, 25]

452

sympy.ntheory.residue_ntheory.quadratic_residues(p)

453

454

Returns the list of quadratic residues.

455

In [ ]:

>>> from sympy.ntheory.residue_ntheory import quadratic_residues >>> quadratic_residues(7) [0, 1, 2, 4]

456

sympy.ntheory.residue_ntheory.nthroot_mod(a, n, p, all_roots=False)

457

458

find the solutions to "x**n = a mod p"

459

a : integer n : positive integer p : positive integer all_roots : if False returns the smallest root, else the list of roots

460

In [ ]:

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

461

sympy.ntheory.residue_ntheory.is_nthpow_residue(a, n, m)

462

463

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

464

465

Hackman "Elementary Number Theory" (2009), page 76

466

sympy.ntheory.residue_ntheory.is_quad_residue(a, p)

467

468

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)]). If "p" is an odd prime, an iterative method is used to make the determination:

469

In [ ]:

>>> from sympy.ntheory import is_quad_residue >>> list(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]

470

legendre_symbol, jacobi_symbol

471

sympy.ntheory.residue_ntheory.legendre_symbol(a, p)

472

473

474

0 if a is multiple of p

475

476

1 if a is a quadratic residue of p

477

478

-1 otherwise

479

p should be an odd prime by definition

480

In [ ]:

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

481

is_quad_residue, jacobi_symbol

482

sympy.ntheory.residue_ntheory.jacobi_symbol(m, n)

483

484

Returns the product of the legendre_symbol(m, p) for all the prime factors, p, of n.

485

486

0 if m cong 0 mod(n)

487

488

1 if x**2 cong m mod(n) has a solution

489

490

-1 otherwise

491

In [ ]:

>>> from sympy.ntheory import jacobi_symbol, legendre_symbol >>> from sympy import Mul, S >>> jacobi_symbol(45, 77) -1 >>> jacobi_symbol(60, 121) 1

492

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

493

In [ ]:

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

494

is_quad_residue, legendre_symbol

495