A clean version of the double factorial function including unit tests.

DoubleFactorial.R

#' Double factorial function
#' The double factorial is the product of all the integers from 1 up to x 
#' that have the same parity (odd or even) as x (see https://t1p.de/byye).
#'
#' @param x Numeric vector. Will be coerced to integer.
#'
#' @note dfactorial is defined for x in [0, 300]
#' @references Gould, H. & Quaintance, J. (2012). Double Fun with Double 
#' Factorials. Mathematics Magazine, 85(3), 177-192. doi:10.4169/math.mag.85.3.177
dfactorial <- function(x) {
  x <- as.integer(x)
  if(any(x < 0, na.rm = TRUE)) stop("Double factorial not possible for values < 1")
  
  Factors <- lapply(x, 
                    function(z) {
                      if(is.na(z)) return(NA)
                      if(z != 0) return(seq(z, 1, by = -2))
                      1
                      })
  Results <- lapply(Factors, prod, na.rm = FALSE)
  
  return(unlist(Results))
}

Test_DoubleFactorial.R

library(testthat)
source("./dfactorial.R")

test_that("dfactorial", {
  
  # Special case x = 0
  expect_identical(dfactorial(0), 1 )
  
  # Double factorial for even numbers
  # taken from https://oeis.org/A000165/list
  e <- c(1, 2, 8, 48, 384, 3840,46080, 645120, 10321920,
         185794560, 3715891200, 81749606400, 1961990553600,
         51011754393600, 1428329123020800,42849873690624000,
         1371195958099968000, 46620662575398912000,
         1678343852714360832000, 63777066403145711616000)
  
  expect_equal(dfactorial(seq(0, 19)*2), e )
  
  # Double factorial for odd numbers
  # taken from https://oeis.org/A001147/list
  e <- c(1, 3, 15, 105, 945, 10395, 135135, 2027025, 34459425,
         654729075,13749310575,316234143225,7905853580625,
         213458046676875,6190283353629375,
         191898783962510625,6332659870762850625,
         221643095476699771875,8200794532637891559375)
  expect_equal(dfactorial(seq(1, 37, 2)), e )
  
  # Negative values not possible
  expect_error(dfactorial(-1), "Double factorial not possible for values < 1")
  expect_error(dfactorial(c(2, 3, 4, -3)), "Double factorial not possible for values < 1")
  expect_error(dfactorial(c(2, -9, 3, 4)), "Double factorial not possible for values < 1")
  
  # Can we handle NA?
  e <- c(1, 2, 8, 48, 384, 3840,46080, 645120, NA)
  expect_equal(dfactorial(c(seq(0, 14, 2), NA)), e)
  
  e <- c(NA, 1, 2, 8, 48, 384, 3840,46080, 645120)
  expect_equal(dfactorial(c(NA, seq(0, 14, 2))), e)
  
  # Comparison with factorial
  for(i in 1:100) {
    o <- dfactorial(i) * dfactorial(i-1)
    e <- factorial(i)
    expect_equal(o, e)
  }
  
  o <- dfactorial(1:100) * dfactorial(1:100-1)
  e <- factorial(1:100)
  expect_equal(o, e)
})