modpower: Power Function modulo m

Description

Calculates powers and orders modulo m.

Usage

modpower(n, k, m)
modorder(n, m)

Arguments

n, k, m

Natural numbers, m >= 1.

Value

Natural number.

Details

modpower calculates n to the power of k modulo m. Uses modular exponentiation, as described in the Wikipedia article.

modorder calculates the order of n in the multiplicative group module m. n and m must be coprime. Uses brute force, trick to use binary expansion and square is not more efficient in an R implementation.

Examples

Run this code

modpower(2, 100, 7)  #=> 2
modpower(3, 100, 7)  #=> 4
modorder(7, 17)      #=> 16, i.e. 7 is a primitive root mod 17

##  Gauss' table of primitive roots modulo prime numbers < 100
proots <- c(2,  2,  3,  2,  2,  6,  5, 10, 10, 10, 2,  2, 10, 17,  5,  5,
            6, 28, 10, 10, 26, 10, 10,  5, 12, 62, 5, 29, 11, 50, 30, 10)
P <- primes(100)
for (i in seq(along=P)) {
    cat(P[i], "t", modorder(proots[i], P[i]), proots[i], "t", "")
}

##  Lehmann's primality test
lehmann_test <- function(n, ntry = 25) {
    if (!is.numeric(n) || ceiling(n) != floor(n) || n < 0)
        stop("Argument 'n' must be a natural number")
    if (n >= 9e7)
        stop("Argument 'n' should be smaller than 9e7.")

    if (n < 2)                      return(FALSE)
    else if (n == 2)                return(TRUE)
    else if (n > 2 && n %% 2 == 0)  return(FALSE)

    k <- floor(ntry)
    if (k < 1) k <- 1
    if (k > n-2) a <- 2:(n-1)
    else         a <- sample(2:(n-1), k, replace = FALSE) 

    for (i in 1:length(a)) {
        m <- modpower(a[i], (n-1)/2, n)
        if (m != 1 && m != n-1) return(FALSE)
    }
    return(TRUE)
}

##  Examples
for (i in seq(1001, 1011, by = 2))
    if (lehmann_test(i)) cat(i, "\n")
# 1009
system.time(lehmann_test(27644437, 50))    # TRUE
#    user  system elapsed 
#   0.086   0.151   0.235