12.37 math.prime - Prime numbers

Module: math.prime ¶

This module provides utilities related to prime numbers.

Sequence of prime numbers

Variable: *primes* ¶

{math.prime} An infinite lazy sequence of primes.

;; show 10 prime numbers from 100-th one.
(take (drop *primes* 100) 10)
 ⇒ (547 557 563 569 571 577 587 593 599 601)

Function: reset-primes ¶

{math.prime} Once you take a very large prime out of *primes*, all primes before that has been calculated remains in memory, since the head of sequence is held in *primes*. Sometimes you know you need no more prime numbers and you wish those calculated ones to be garbage-collected. Calling reset-primes rebinds *primes* to unrealized lazy sequence, allowing the previously realized primes to be GCed.

Function: primes ¶

{math.prime} Returns a fresh lazy sequence of primes. It is useful when you need certain primes in a short period of time—if you don’t keep a reference to the head of the returned sequence, it will be garbage collected after you’ve done with the primes. (Note that calculation of a prime number needs the sequence of primes from the beginning, so even if your code only keep a reference in the middle of the sequence, the entire sequence will be kept in the thunk within the lazy sequence—you have to release all references in order to make the sequence GCed.)

On the other hand, each sequence returned by primes are realized individually, duplicating calculation.

The rule of thumb is—if you use primes repeatedly throughout the program, just use *primes* and you’ll save calculation. If you need primes one-shot, call primes and abandon it and you’ll save space.

Testing primality

Function: small-prime? n ¶

{math.prime} For relatively small positive integers (below *small-prime-bound*, to be specific), this procedure determines if the input is prime or not, quickly and deterministically. If n is on or above the bound, this procedure returns #f.

This can be used to quickly filter out known primes; it never returns #t on composite numbers (while it may return #f on large prime numbers). Miller-Rabin test below can tell if the input is composite for sure, but it may return #t on some composite numbers.

Variable: *small-prime-bound* ¶

{math.prime} For all positive integers below this value (slightly above 3.4e14 in the current implementation), small-prime? can determines whether it is a prime or not.

Function: miller-rabin-prime? n :key num-tests random-integer ¶

{math.prime} Check if an exact integer n is a prime number, using probabilistic Miller-Rabin algorithm (n must be greater than 1). If this procedure returns #f, n is a composite number. If this procedure returns #t, n is likely a prime, but there’s a small probability that it is a false positive.

Note that if n is smaller than a certain number (*small-prime-bound*), the algorithm is deterministic; if it returns #t, n is certainly a prime.

If n is greater than or equal to *small-prime-bound*, we use a probabilistic test. We choosing random base integer to perform Miller-Rabin test up to 7 times by default. You can change the number of tests by the keyword argument num-tests. The error probability (to return #t for a composite number) is at most (expt 4 (- num-tests)).

For a probabilistic test, miller-rabin-prime? uses its own fixed random seed by default. We chose fixed seed so that the behavior can be reproducible. To change the random sequence, you can provide your own random integer generator to the random-integer keyword argument. It must be a procedure that takes a positive integer k and returns a random integer from 0 to k-1, including.

Function: bpsw-prime? n ¶

{math.prime} Check if an exact integer n is a prime number, using Baillie-PSW primality test (http://www.trnicely.net/misc/bpsw.html). It is deterministic, and returns the definitive answer below 2^64 (around 1.8e19). For larger integers this can return #t on a composite number, although such number hasn’t been found yet. This never returns #f on a prime number.

This is slower than Miller-Rabin but fast enough for casual use, so it is handy when you want a definitive answer below the above range.

Factorization

Function: naive-factorize n :optional divisor-limit ¶

{math.prime} Factorize a positive exact integer n by trying to divide it with all primes up to (sqrt n). Returns a list of prime factors (each of which is equal to or greater than 2), smaller ones first.

(naive-factorize 142857)
  ⇒ (3 3 3 11 13 37)

Note that (naive-factorize 1) is ().

Although this is pretty naive method, this works well as far as any of n’s factors are up to the order of around 1e7. For example, the following example runs in about 0.4sec on 2.4GHz Core2 machine. (The first time will take about 1.3sec to realize lazy prime sequences.)

(naive-factorize 3644357367494986671013))
  ⇒ (10670053 10670053 32010157)

Of course, if n includes any factors above that order, the performance becomes abysmal. So it is better to use this procedure below 1e14 or so.

Alternatively, you can give divisor-limit argument that specifies the upper bound of the prime number to be tried. If it is given, naive-factorize leaves a factor f as is if it can’t be divided by any primes less than or equal to divisor-limit. So, the last element of the returned list may be composite number. This is handy to exclude trivial factors before applying more sophisticated factorizing algorithms.

(naive-factorize 825877877739 1000)
  ⇒ (3 43 6402154091)

;; whereas
(naive-factorize 825877877739)
  ⇒ (3 43 4591 1394501)

The procedure also memoizes the results on smaller n to make things faster.

Function: mc-factorize n ¶

{math.prime} Factorize a positive exact integer n using the algorithm described in R. P. Brent, An improved Monte Carlo factorization algorithm, BIT 20 (1980), 176-184. http://maths-people.anu.edu.au/~brent/pub/pub051.html.

This one is capable to handle much larger range than naive-factorize, somewhere around 1e20 or so.

Since this method is probabilistic, the execution time may vary on the same n. But it will always return the definitive results as far as every prime factor of n is smaller than 2^64.

At this moment, if n contains a prime factor greater than 2^64, this routine would keep trying factorizing it forever. Practical applications should have some means to interrupt the function and give it up after some time bounds. This will be addressed once we have deterministic primality test.

Miscellaneous

Function: jacobi a n ¶

{math.prime} Calculates Jacobi symbol (a/n) (http://en.wikipedia.org/wiki/Jacobi_symbol).

Function: totient n ¶

{math.prime} Euler’s totient function of nonnegative integer n.

The current implementation relies on mc-factorize above, so it may take very long if n contains large prime factors.