isPrimeRcpp: Vectorized Primality Test

View source: R/NumberTheory.R

isPrimeRcppR Documentation

Vectorized Primality Test


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:


isPrimeRcpp(v, namedVector = FALSE, nThreads = NULL)



Vector of integers or numeric values.


Logical flag. If TRUE, a named vector is returned. The default is FALSE.


Specific number of threads to be used. The default is NULL.


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.


See Also

primeFactorize, isprime


## check the primality of a single number

## check the primality of every number in a vector

mySamp <- sample(10^13, 10)

## return named vector for easy identification
isPrimeRcpp(mySamp, namedVector = TRUE)

## Using nThreads
system.time(isPrimeRcpp(mySamp, nThreads = 2))

RcppAlgos documentation built on Oct. 3, 2023, 1:07 a.m.