12.1. Functions

In R you can write your own functions. In the next example we define a function $sq$ that takes one argument $x$ and returns $x × x$. We can then apply this function to a vector using sapply() or we can vectorize it.

> myvec <- c(1:9)
> myvec
[1] 1 2 3 4 5 6 7 8 9
> sq <- function(x) { x*x }   # or sq <- function(x) { return(x*x) }
> sapply(myvec, sq)
[1]  1  4  9 16 25 36 49 64 81
# vectorize function
> sq <- Vectorize(sq)
> sq(c(1:9))
[1]  1  4  9 16 25 36 49 64 81
> 
 

12.1.1. Example with prime numbers

Remember that a prime number $x$ is a natural number that has only 2 divisors : 1 and $x$. Some people add the constraint that 1 and $x$ must be distinct, the consequence is that 1 is not considered as a prime number.

12.1.1.a  First version

A first version of the $isprime(x)$ function that returns TRUE if $x$ is a prime number (FALSE otherwise) can be written like a function of the C programming language with some restrictions.

We are looking for any divisor between 2 and $x-1$

# isprime version 1
# look for any divisor between 2..(x-1)
isprime <- function(x) {
	if (x <= 1) return(F)
	if (x <= 3) return(T)
	for (i in 2:(x-1)) {
		 if ((x %% i) == 0) return(F)
	}
	return(T) 
}
 
numbers <- c(1:13)
select <- sapply(numbers, isprime)
result= numbers[select]
result

If you execute this script (save it as isprime_v1.rs for example), you obtain the following result:

richer@nano:~/dev/R$ Rscript isprime_v1.rs
[1]  2  3  5  7 11 13

12.1.1.b  Second version

We can improve the search of a divisor knowing that it is not necessary to go over $⌈√x⌉$ and that if $x$ is a multiple of 2 then it is not necessary to check for even numbers. Note that $⌈x⌉$ is the ceiling function that maps $x$ to the least integer that is greater than or equal to $x$.

# isprime version 2
# look for any divisor between 3..ceiling(sqrt(x)) by increment of 2
 
isprime <- function(x) {
	if (x <= 1) return(F)
	if (x <= 3) return(T)
	if ((x %% 2) == 0) return (F) 
	for (i in seq(3, ceiling(sqrt(x)) ,2)) {
		 if ((x %% i) == 0) return(F)
	}
	return(T) 
}
 
numbers <- c(1:13)
select <- sapply(numbers, isprime)
result= numbers[select]
result

12.1.1.c  Third version

We can shorten the function using the any(condition) function of the R language, it will return true if the condition inside any() can become true. So by writing any(x %% c(3:ceiling(sqrt(x-1)),2) == 0) we say "are there any numbers in the interval $3..⌈√x⌉$ by increment of 2 that divides $x$" ? In other words, has $x$ any divisor in this interval.

# isprime version 3
# use the any function
 
isprime <- function(x) {
	if (x <= 1) return(F)
	if (x <= 3) return(T)
	if ((x %% 2) == 0) return (F)
	return( !any(x %% c(3:ceiling(sqrt(x-1)),2) == 0) )
}
 
numbers <- c(1:13)
select <- sapply(numbers, isprime)
result= numbers[select]
result