isPrimeRcpp | R Documentation |
Implementation of the Miller-Rabin primality test. Based on the "mp_prime_p" function from the "factorize.c" source file found in the gmp library: https://gmplib.org.
isPrimeRcpp(v, namedVector = FALSE, nThreads = NULL)
v |
Vector of integers or numeric values. |
namedVector |
Logical flag. If |
nThreads |
Specific number of threads to be used. The default is |
The Miller-Rabin primality test is a probabilistic algorithm that makes heavy use of modular exponentiation. At the heart of modular exponentiation is the ability to accurately obtain the remainder of the product of two numbers \pmod p
.
With the gmp library, producing accurate calculations for problems like this is trivial because of the nature of the multiple precision data type. However, standard C++ does not afford this luxury and simply relying on a strict translation would have limited this algorithm to numbers less than \sqrt 2^{63} - 1
(N.B. We are taking advantage of the signed 64-bit fixed width integer from the stdint library in C++. If we were confined to base R, the limit would have been \sqrt 2^{53} - 1
). RcppAlgos::isPrimeRcpp gets around this limitation with a divide and conquer approach taking advantage of properties of arithmetic.
The problem we are trying to solve can be summarized as follows:
(x_1 * x_2) \pmod p
Now, we rewrite x_2
as x_2 = y_1 + y_2 + \dots + y_n
, so that we obtain:
(x_1 * y_1) \pmod p + (x_1 * y_2) \pmod p + \dots + (x_1 * y_n) \pmod p
Where each product (x_1 * y_j)
for j <= n
is smaller than the original x_1 * x_2
. With this approach, we are now capable of handling much larger numbers. Many details have been omitted for clarity.
For a more in depth examination of this topic see Accurate Modular Arithmetic with Double Precision.
Returns a named/unnamed logical vector. If an index is TRUE
, the number at that index is prime, otherwise the number is composite.
The maximum value for each element in v
is 2^{53} - 1
.
Conrad, Keith. "THE MILLER-RABIN TEST." https://www.math.uconn.edu/~kconrad/blurbs/ugradnumthy/millerrabin.pdf.
primeFactorize
, isprime
## check the primality of a single number
isPrimeRcpp(100)
## check the primality of every number in a vector
isPrimeRcpp(1:100)
set.seed(42)
mySamp <- sample(10^13, 10)
## return named vector for easy identification
isPrimeRcpp(mySamp, namedVector = TRUE)
## Using nThreads
system.time(isPrimeRcpp(mySamp, nThreads = 2))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.