Rcpp version of mpoly and mpoly arithmetic

DESCRIPTION

@@ -4,8 +4,10 @@ Title: Symbolic Computation and More with Multivariate Polynomials
 Version: 1.0.5.902
 URL: https://github.com/dkahle/mpoly
 BugReports: https://github.com/dkahle/mpoly/issues
-Authors@R: person("David", "Kahle", email = "david.kahle@gmail.com", role = c("aut", "cre"))
+Authors@R: c(person("David", "Kahle", email = "david.kahle@gmail.com", role = c("aut", "cre")),
+             person("Blagoje", "Ivanovic", email = "blagoje.r.ivanovic@gmail.com", role = c("aut")))
 Description: Symbolic computing with multivariate polynomials in R.
+LinkingTo: Rcpp
 Imports:
     stringr (>= 1.0.0),
     partitions,
@@ -14,10 +16,12 @@ Imports:
     ggplot2,
     polynom,
     orthopolynom,
-    tidyr
+    tidyr,
+    Rcpp
 Suggests:
     testthat,
     magrittr,
     dplyr
+SystemRequirements: C++11
 License: GPL-2
 RoxygenNote: 6.0.1

NAMESPACE

@@ -67,6 +67,8 @@ export(swap)
 export(totaldeg)
 export(tuples)
 export(vars)
+import(Rcpp)
+importFrom(Rcpp,sourceCpp)
 importFrom(ggplot2,aes)
 importFrom(ggplot2,ggplot)
 importFrom(orthopolynom,chebyshev.c.polynomials)
@@ -101,3 +103,4 @@ importFrom(stringr,str_split)
 importFrom(stringr,str_sub)
 importFrom(stringr,str_trim)
 importFrom(tidyr,gather)
+useDynLib(mpoly)

R/RcppExports.R

@@ -0,0 +1,15 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+reducePoly <- function(poly, vars) {
+    .Call('_mpoly_reducePoly', PACKAGE = 'mpoly', poly, vars)
+}
+
+mpolyMult <- function(e1, e2) {
+    .Call('_mpoly_mpolyMult', PACKAGE = 'mpoly', e1, e2)
+}
+
+mpolyPow <- function(e1, e2, vars) {
+    .Call('_mpoly_mpolyPow', PACKAGE = 'mpoly', e1, e2, vars)
+}
+

R/help.R

@@ -1,7 +1,9 @@
 #' Multivariate polynomials in R.
 #'
 #' A package for symbolic computation and more with multivariate polynomials
-#'
+#' @import Rcpp
+#' @importFrom Rcpp sourceCpp
+#' @useDynLib mpoly
 #' @importFrom stringr str_c str_detect str_sub str_replace_all str_replace
 #'   str_count str_dup str_extract str_extract_all str_locate str_locate_all
 #'   str_split str_sub<- str_trim fixed

R/mpoly.R

@@ -50,42 +50,7 @@ mpoly <- function(list, varorder){
   flatList <- flatList[names(flatList) != "coef"]
   if(!all(is.wholenumber(flatList)) || any(flatList < 0)){
     stop("degrees must be nonnegative integers.", call. = FALSE)
-  }  
-
-
-  # give terms without a coefficient a coef of 1
-  addCoefIfMissing <- function(v){
-    if(any("coef" == names(v))) return(v)
-    c(v, coef = 1)
   }
-  list <- lapply(list, addCoefIfMissing)
-
-
-
-  ## organize 
-
-  # remove terms with coef 0
-  list <- filterOutZeroTerms(list)
-
-
-  # remove 0 degrees, combine like degrees, single coef as rightmost element
-  list <- lapply(list, function(v){  	
-    # separate vardegs from coefs
-    coef_ndx <- which(names(v) == "coef")
-    coefs <- v[coef_ndx]    
-    v     <- v[-coef_ndx]
-
-  	# combine like degrees (sum)
-  	if(length(names(v)) != length(unique(names(v)))) v <- fastNamedVecTapply(v, sum)   
-
-  	# combine like coefficients (product)
-  	coefs <- c(coef = prod(coefs))
-
-  	# remove zero degree elements, combine and return  
-    c(v[v != 0], coefs)
-  })
-
-
 
   ## set intrinsic varorder - done again after 0 degrees are removed
   vars <- unique(names(flatList))
@@ -101,126 +66,10 @@ mpoly <- function(list, varorder){
     }    
     vars <- varorder
   }
-
-
-  # sort variables in terms
-  list <- lapply(list, function(v){
-    p <- length(v) - 1L
-    if(p == 0L) return(v)
-    c( (v[1:p])[intersect(vars, names(v[1:p]))], v["coef"] )
-  })   
-
-
-
-  ## prepare to check if like terms are present
-  monomials <- vapply(list, function(v){
-  	p <- length(v) - 1 # remove coef on monomials
-    paste(names(v[1:p]), v[1:p],  sep = "", collapse = "")
-  }, character(1))  
-  # e.g. c("x1", "y1", "y1", "x2z1", "y4t3", "x1", "x1", "coef5")
-
-  unique_monomials <- unique(monomials)
-
-
-
-  ## check if like terms are present and, if so, correct  
-  if(length(monomials) != length(unique_monomials)){
-
-    matchedMonomials <- match(monomials, unique_monomials)
-    matchedMonomials <- factor(matchedMonomials, levels = 1:max(matchedMonomials))
-    ndcs2combine     <- split.default(1:length(list), matchedMonomials)
-
-    list <- lapply(ndcs2combine, function(v){
-      if(length(v) == 1) return(list[[v]])      
-      coef <- sum(vapply(list[v], `[`, double(1), length(list[[v[1]]])))
-      v <- list[[v[1]]]
-      v["coef"] <- coef
-      v
-    })
-
-    names(list) <- NULL # i.e. list <- unname(list)
-  }
-
-
-
-  ## combine constant terms
-  ## mpoly(list(c(x = 1, coef = 1), c(coef = 1), c(coef = 2)))    
-  nonConstantTerms <- fastFilter(isNotLengthOne, list)
-  constantTerms    <- fastFilter(isLengthOne, list)
-
-  if(length(constantTerms) > 0){
-    list <- c(
-      nonConstantTerms, 
-      list(Reduce(`+`, constantTerms))
-    )  
-  }
-
-
-
-  ## re-organize after like-terms combined
-  list <- filterOutZeroTerms(list)
-
 
+  list <- reducePoly(list, vars)
 
   ## return classed list
   class(list) <- "mpoly"
   list  
 }
-
-
-
-
-
-
-
-
-
-
-
-filterOutZeroTerms <- function(list){
-
-  # a term is zero if any of its coefficients are 0
-  hasNonZeroCoef <- function(v) all(v["coef"] != 0)
-
-  # remove terms with coef 0
-  list <- fastFilter(hasNonZeroCoef, list)
-
-  # if all terms have been eliminated, recreate it with a 0 coef
-  if(length(list) == 0) list <- list(c(coef = 0))  
-
-  # return
-  list  
-}
-
-
-
-
-
-
-
-isLengthOne    <- function(x) length(x) == 1L
-isNotLengthOne <- function(x) length(x) > 1L
-fastFilter     <- function(f, x) x[vapply(x, f, logical(1))]
-
-
-
-
-
-
-
-fastNamedVecTapply <- function(x, f, type = double(1)){
-  uniqueNames  <- unique(names(x))
-  matchedNames <- match(names(x), uniqueNames) # indices
-  matchedNames <- factor(matchedNames, levels = 1:max(matchedNames))
-  groupIndices <- split.default(1:length(x), matchedNames)
-  out <- vapply(groupIndices, function(ndcs) f(x[ndcs]), type)
-  names(out) <- uniqueNames
-  out
-}
-# x <- 1:10
-# names(x) <- sample(letters[1:3], 10, replace = TRUE)
-# fastNamedVecTapply(x, sum)
-# tapply(x, names(x), sum)
-
-
-

R/mpolyArithmetic.R

@@ -63,14 +63,25 @@ NULL
 #' @rdname mpolyArithmetic
 #' @export
 `-.mpoly` <- function(e1, e2){
-
+
+  ## if either is constant, mpoly it
+  if(!is.list(e1)){
+    stopifnot(is.numeric(e1) && length(e1) == 1)
+    e1 <- list(c(coef = e1))
+  }
+
+  if(!is.list(e2)){
+    stopifnot(is.numeric(e2) && length(e2) == 1)
+    e2 <- list(c(coef = e2))
+  }
+
   ## flip coefficients of each term
   e2 <- lapply(e2, function(v){
     v["coef"] <- -v["coef"]
     v
   })
   class(e2) <- "mpoly"
-  
+
   ## add
   e1 + e2	
 }
@@ -83,36 +94,24 @@ NULL
 #' @rdname mpolyArithmetic
 #' @export
 `*.mpoly` <- function(e1, e2){
-
   ## allow for multiplication by a constant
   if(is.numeric(e1) && length(e1) == 1) e1 <- mpoly(list(c(coef = e1)))
   if(is.numeric(e2) && length(e2) == 1) e2 <- mpoly(list(c(coef = e2)))
-
-
-
+
   ## argument check
   stopifnot(is.mpoly(e1))
   stopifnot(is.mpoly(e2))
-
-
-
+
   ## multiply
-  list <- lapply(e1, function(v){
-    lapply(e2, function(z){
-      c(v, z)  
-    })	
-  })
-
-
-
+  list <- mpolyMult(e1, e2)
   ## return
-  mpoly( unlist(list, recursive = FALSE) )
+  mpoly( list )
 }	
 
 
-l <- list(
-  c(x = 1, coef = 1, x = 1, coef = 1)  
-)
+# l <- list(
+#   c(x = 1, coef = 1, x = 1, coef = 1)
+# )
 
 
 
@@ -125,16 +124,14 @@ l <- list(
 #' @rdname mpolyArithmetic
 #' @export
 `^.mpoly` <- function(e1, e2){
-  
+
   if(!is.mpoly(e1)) stop('e1 must be of class mpoly.', call. = FALSE)
-  
+
   if(!is.wholenumber(e2) || e2 < 0) stop('exponent must be a positive integer.', call. = FALSE)
-
-  out <- mpoly(list(c(coef = 1)))
-
+
   if(e2 == 0) return(out)
-    
-  for(k in 1:e2) out <- out * e1
-  
-  out
+
+  out <- mpolyPow(e1, e2, vars(e1))
+
+  mpoly(out)
 }

R/mpolyEqual.R

@@ -51,5 +51,7 @@
 }
 
 
-
+all.equal.mpoly <- function(e1, e2) {
+  e1 == e2
+}
 

inst/include/helpers.h

@@ -0,0 +1,24 @@
+#ifndef mpoly_HELPERS_H
+#define mpoly_HELPERS_H
+
+using namespace Rcpp;
+
+// concatenate string vectors
+std::vector<std::string> scat(std::vector<std::string> x,
+                              std::vector<std::string> y);
+
+// concatenate two numeric vectors (as c() in R)
+NumericVector vcat(const NumericVector& x, const NumericVector& y);
+
+// helper function to convert types to string
+template <typename T>
+std::string to_string(T const& value) {
+  std::stringstream sstr;
+  sstr << value;
+  return sstr.str();
+};
+
+// Vectorized equality operator for string vs Character vector
+LogicalVector vecstreq(const CharacterVector& vec, std::string str);
+
+#endif

inst/include/reducepoly.h

@@ -0,0 +1,8 @@
+#ifndef mpoly_REDUCEPOLY_H
+#define mpoly_REDUCEPOLY_H
+
+using namespace Rcpp;
+
+List reducePoly(const List& poly, const CharacterVector& vars);
+
+#endif

src/Makevars

@@ -0,0 +1 @@
+PKG_CPPFLAGS =  -I../inst/include/