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 mod 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:
(x1 * x2) mod p
Now, we rewrite x2 as x2 = y1 + y2 + … + yn, so that we obtain:
(x1 * y1) mod p + (x1 * y2) mod p + … + (x1 * yn) mod p
Where each product x1 * yj for j <= n is smaller than the original x1 * x2. 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." http://www.math.uconn.edu/~kconrad/blurbs/ugradnumthy/millerrabin.pdf.
primeFactorize
, isprime
, 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.