diff --git a/DESCRIPTION b/DESCRIPTION
index 07fb783..4aac1ff 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: mpoly
 Type: Package
 Title: Symbolic Computation and More with Multivariate Polynomials
-Version: 1.1.0.902
+Version: 1.1.0.903
 URL: https://github.com/dkahle/mpoly
 BugReports: https://github.com/dkahle/mpoly/issues
 Authors@R: person("David", "Kahle", email = "david@kahle.io", role = c("aut", "cre"))
diff --git a/NEWS.md b/NEWS.md
index 85bd400..b16e4bd 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,5 +2,8 @@
 
 ## Minor improvements and fixes
 
+*   __mpoly__ now is extensively tested, and many bugs have been found and fixed.
 *   `bernsteinApprox()` is now deprecated in favor of `bernstein_approx()`.
 *   `bezierFunction()` is now deprecatd in favor of `bezier_function()`.
+*   More arithmetic is available for `mpolyList` objects, such as single 
+    constant times mpolyList.
diff --git a/R/as.function.mpoly.R b/R/as.function.mpoly.R
index 0983c93..66df116 100644
--- a/R/as.function.mpoly.R
+++ b/R/as.function.mpoly.R
@@ -44,6 +44,14 @@
 #' f <- as.function(p)
 #' f(10) # -> 24
 #'
+#' bernstein(1, 2)
+#' s <- seq(0, 1, .01)
+#' as.function(bernstein(1, 2))(s)
+#' plot(
+#'   s,
+#'   as.function(bernstein(1, 2))(s)
+#' )
+#'
 #'
 #' as.function(mp("x + xx"))
 #' as.function(mp("x + xx"), squeeze = FALSE)
@@ -64,7 +72,7 @@ as.function.mpoly <- function(x, varorder = vars(x), vector = TRUE, silent = FAL
  
   
   # deal with constant polynomials
-  if (is.constant(x)) return( function(.) unlist(x)[["coef"]] )
+  if (is.constant(x)) return( function(.) 0*. + unlist(x)[["coef"]] )
   
   
   # print poly with stars
@@ -136,8 +144,21 @@ as.function.bernstein <- function(x, ...){
   k <- attr(x, "bernstein")$k
   n <- attr(x, "bernstein")$n
   
-  # return exp"d log function
-  function(.) exp(lchoose(n, k) + k*log(.) + (n-k)*log(1-.))
+  # return exp'd log function
+  function(.) {
+    
+    out <- vector("numeric", length = length(.))
+
+    non_pos_ndcs <- (. <= 0)
+    sup_one_ndcs <- (. >= 1)
+    if (any(non_pos_ndcs)) out[non_pos_ndcs] <- choose(n, k) * .[non_pos_ndcs]^k * (1 - .[non_pos_ndcs])^(n-k)
+    if (any(sup_one_ndcs)) out[sup_one_ndcs] <- choose(n, k) * .[sup_one_ndcs]^k * (1 - .[sup_one_ndcs])^(n-k)
+    if (!all(sup_one_ndcs | sup_one_ndcs)) out[0 < . & . < 1] <- exp(
+      lchoose(n, k)  +  k*log(.[0 < . & . < 1])  +  (n-k)*log(1-.[0 < . & . < 1])
+    )
+    
+    out
+  }
   
 }
 
@@ -147,23 +168,23 @@ as.function.bernstein <- function(x, ...){
 
 
 
-as.function.jacobi <- function(x, ...){
-  return(as.function.mpoly(x)) ## below is broken.
-  
-  # grab bernstein values
-  d <- attr(x, "jacobi")$degree
-  k <- attr(x, "jacobi")$kind
-  i <- attr(x, "jacobi")$indeterminate
-  n <- attr(x, "jacobi")$normalized
-  a <- attr(x, "jacobi")$alpha
-  b <- attr(x, "jacobi")$beta
-  
-  # return exp'd log function #
-  #http://en.wikipedia.org/wiki/Jacobi_polynomials function(.)
-  #pochhammer(a+1, d) / factorial(d) * hypergeo(-d, 1+a+b+d, a+1,
-  #(1-.)/2)
-  
-}
+# as.function.jacobi <- function(x, ...){
+#   return(as.function.mpoly(x)) ## below is broken.
+#   
+#   # grab bernstein values
+#   d <- attr(x, "jacobi")$degree
+#   k <- attr(x, "jacobi")$kind
+#   i <- attr(x, "jacobi")$indeterminate
+#   n <- attr(x, "jacobi")$normalized
+#   a <- attr(x, "jacobi")$alpha
+#   b <- attr(x, "jacobi")$beta
+#   
+#   # return exp'd log function #
+#   #http://en.wikipedia.org/wiki/Jacobi_polynomials function(.)
+#   #pochhammer(a+1, d) / factorial(d) * hypergeo(-d, 1+a+b+d, a+1,
+#   #(1-.)/2)
+#   
+# }
 
 
 
diff --git a/R/as.mpoly.R b/R/as.mpoly.R
index 420dbda..dea8919 100644
--- a/R/as.mpoly.R
+++ b/R/as.mpoly.R
@@ -66,8 +66,9 @@ as.mpoly <- function(x, ...) UseMethod("as.mpoly")
 
 
 #' @export  
-as.mpoly.default <- function(x, ...) 
-  stop("object not supported.  see ?as.mpoly for details.")
+as.mpoly.default <- function(x, ...) {
+  stop(sprintf("objects of class %s not supported by as.mpoly().", class(x)), call. = FALSE)
+}
 
 
 
@@ -116,7 +117,7 @@ parse_model_poly <- function(s) {
   inside <- str_sub(inside, 6, -2)
   inside <- str_split(inside, ", ")[[1]]
   if(!any(inside == "raw = TRUE")) {
-    stop("poly() statements currently must contain raw = TRUE.")
+    stop("poly() statements currently must contain raw = TRUE.", call. = FALSE)
   }
   vars <- inside[!(str_detect(inside, "=") | !str_detect(inside, "[:alpha:]"))]
   
diff --git a/R/bernstein.R b/R/bernstein.R
index 90574c6..bfb24ed 100644
--- a/R/bernstein.R
+++ b/R/bernstein.R
@@ -59,6 +59,8 @@
 #' @export
 bernstein <- function(k, n, indeterminate = "x"){  
   
+  stopifnot(k <= n)
+  
   ## make it possible for vector k args
   if(length(k) > 1){
     listOPolys <- lapply(k, function(.) bernstein(., n, indeterminate))
diff --git a/R/mp.R b/R/mp.R
index 8fd7dd7..6612e8b 100644
--- a/R/mp.R
+++ b/R/mp.R
@@ -69,9 +69,31 @@ mp <- function (string, varorder, stars_only = FALSE) {
   
   # deal with mpolyLists
   if (length(string) > 1) {
-    return(structure(lapply(string, mp), class = "mpolyList"))
+    
+    # do basic mpoly parsing
+    ps <- structure(
+      lapply(string, mp, stars_only = stars_only), 
+      class = "mpolyList"
+    )
+
+    # enforce varorder if present
+    if (!missing(varorder)) {
+      for (k in seq_along(ps)) {
+        ps[[k]] <- structure(
+          lapply(ps[[k]], function(term) {
+            vars_in_term <- intersect(varorder, names(term))
+            term[c(vars_in_term, "coef")]
+          }),
+          class = "mpoly"
+        )
+      }
+    }
+    
+    # return early
+    return(ps)
   } 
  
+  
   # clean spaces if needed
   if(!stars_only)  {
     # put *s in for spaces, twice for situations like "x y z"
@@ -97,9 +119,18 @@ mp <- function (string, varorder, stars_only = FALSE) {
   if (is.numeric(p)) return( structure(list(c(coef = p)), class = "mpoly") )
   
   # reorder if needed and return
-  if (!missing(varorder)) p <- reorder.mpoly(p)
-  p
+  if (!missing(varorder)) {
+    p <- structure(
+      lapply(p, function(term) {
+        vars_in_term <- intersect(varorder, names(term))
+        term[c(vars_in_term, "coef")]
+      }),
+      class = "mpoly"
+    )
+  }
   
+  # return
+  p
 }
 
 
@@ -179,18 +210,18 @@ mp <- function (string, varorder, stars_only = FALSE) {
 # parse_parenthetical_polynomial(string)
 # parse_parenthetical_polynomial("x ((x+y) + 2)")
 # parse_parenthetical_polynomial("-(x + y) + 2 x (x + y)^2 + 3 y")
-parse_parenthetical_polynomial <- function(string){
-  
-  # fix term joins
-  terms <- extract_polynomial_terms(string)
-  
-  # parse into mpolys
-  mpolys <- lapply(terms, parse_parenthetical_term)
-  
-  # add and return
-  Reduce(`+.mpoly`, mpolys)
-  
-}
+# parse_parenthetical_polynomial <- function(string){
+#   
+#   # fix term joins
+#   terms <- extract_polynomial_terms(string)
+#   
+#   # parse into mpolys
+#   mpolys <- lapply(terms, parse_parenthetical_term)
+#   
+#   # add and return
+#   Reduce(`+.mpoly`, mpolys)
+#   
+# }
 
 
 
@@ -233,47 +264,47 @@ parse_parenthetical_polynomial <- function(string){
 # parse_parenthetical_term(string)
 # parse_parenthetical_term("6 (x)")
 # parse_parenthetical_term("6.18033988749895 (x)")
-parse_parenthetical_term <- function(string){
-  
-  # short circuit if simpler
-  if(!contains_parenthetical_expression(string)) 
-    return(parse_nonparenthetical_term(string))
-  
-  # break into parenthetical pieces ("bubbles")
-  pieces <- term_parentheticals(string)
-  pieces <- pieces[pieces != ""]
-  
-  # mpoly pieces
-  mpolys <- lapply(pieces, function(piece){
-    
-    # identify expression and exponent components
-    expr <- str_extract(piece, "\\([-()\\w+ \\^.]+\\)")
-    expr <- str_sub(expr, 2, -2) # take off parens
-    
-    # check for exponent on the outer parenthetical
-    last_paren_ndx     <- nchar(piece) - str_locate(str_rev(piece), fixed(")"))[[1]] + 1L
-    string_after_paren <- str_sub(piece, last_paren_ndx + 1L) # "" or "^3"
-    
-    # if "^3", extract, otherwise 1
-    if(str_detect(string_after_paren, fixed("^"))){
-      exponent <- as.numeric(str_rev(str_extract(str_rev(string_after_paren), "[0-9]+"))) # gets first
-    } else {
-      exponent <- 1
-    }
-    
-    # parse
-    if(contains_nested_parenthetical_expression(piece)){
-      parse_parenthetical_polynomial(expr)^exponent  
-    } else {
-      parse_nonparenthetical_polynomial(expr)^exponent
-    }
-    
-  })
-  
-  # product and return
-  Reduce(`*.mpoly`, mpolys)
-  
-}
+# parse_parenthetical_term <- function(string){
+#   
+#   # short circuit if simpler
+#   if(!contains_parenthetical_expression(string)) 
+#     return(parse_nonparenthetical_term(string))
+#   
+#   # break into parenthetical pieces ("bubbles")
+#   pieces <- term_parentheticals(string)
+#   pieces <- pieces[pieces != ""]
+#   
+#   # mpoly pieces
+#   mpolys <- lapply(pieces, function(piece){
+#     
+#     # identify expression and exponent components
+#     expr <- str_extract(piece, "\\([-()\\w+ \\^.]+\\)")
+#     expr <- str_sub(expr, 2, -2) # take off parens
+#     
+#     # check for exponent on the outer parenthetical
+#     last_paren_ndx     <- nchar(piece) - str_locate(str_rev(piece), fixed(")"))[[1]] + 1L
+#     string_after_paren <- str_sub(piece, last_paren_ndx + 1L) # "" or "^3"
+#     
+#     # if "^3", extract, otherwise 1
+#     if(str_detect(string_after_paren, fixed("^"))){
+#       exponent <- as.numeric(str_rev(str_extract(str_rev(string_after_paren), "[0-9]+"))) # gets first
+#     } else {
+#       exponent <- 1
+#     }
+#     
+#     # parse
+#     if(contains_nested_parenthetical_expression(piece)){
+#       parse_parenthetical_polynomial(expr)^exponent  
+#     } else {
+#       parse_nonparenthetical_polynomial(expr)^exponent
+#     }
+#     
+#   })
+#   
+#   # product and return
+#   Reduce(`*.mpoly`, mpolys)
+#   
+# }
 
 
 
@@ -339,26 +370,26 @@ parse_parenthetical_term <- function(string){
 # parse_nonparenthetical_polynomial("x    -       y+-xy")
 # parse_nonparenthetical_polynomial("1e-2 x")
 # parse_nonparenthetical_polynomial("1e+2 x")
-parse_nonparenthetical_polynomial <- function(string){
-  
-  # check to see if it's a single term
-  if (!str_detect(string, "[+]") && !str_detect(str_sub(string, 2), "[-]")) {
-    return(parse_nonparenthetical_term(string))
-  }
-  
-  # regularize term joins (deal with minuses)
-  string <- fix_term_joins(string)  
-  
-  # split polynomial
-  terms <- str_split(string, fixed(" + "))[[1]]
-  
-  # parse terms
-  mpolyTerms <- lapply(terms, parse_nonparenthetical_term)
-  
-  # combine and return
-  Reduce(`+.mpoly`, mpolyTerms)
-  
-}
+# parse_nonparenthetical_polynomial <- function(string){
+#   
+#   # check to see if it's a single term
+#   if (!str_detect(string, "[+]") && !str_detect(str_sub(string, 2), "[-]")) {
+#     return(parse_nonparenthetical_term(string))
+#   }
+#   
+#   # regularize term joins (deal with minuses)
+#   string <- fix_term_joins(string)  
+#   
+#   # split polynomial
+#   terms <- str_split(string, fixed(" + "))[[1]]
+#   
+#   # parse terms
+#   mpolyTerms <- lapply(terms, parse_nonparenthetical_term)
+#   
+#   # combine and return
+#   Reduce(`+.mpoly`, mpolyTerms)
+#   
+# }
 
 
 
@@ -380,72 +411,72 @@ parse_nonparenthetical_polynomial <- function(string){
 # parse_nonparenthetical_term("1.5x")
 # parse_nonparenthetical_term("1.5^2x")
 # parse_nonparenthetical_term("1e-2 x") # correctly error
-parse_nonparenthetical_term <- function(string){
-  
-  # fix spaces around exponents "x ^ 2" -> "x^2"
-  string <- str_replace_all(string, " *\\^ *", "^")
-  
-  # fix spaces around minuses "x  -  2" -> "x-2"
-  string <- str_replace_all(string, " *- *", "-")
-  
-  # split based on spaces
-  parts <- str_split(string, " ")[[1]]
-  parts <- parts[nchar(parts) > 0] # for "2        -2"
-  
-  # if more than one negative provided error
-  if(str_detect(str_sub(string, 2), fixed("-")))
-    stop("Negative signs are only allowed at the beginning of terms.", call. = FALSE)
-  
-  # fix, e.g. "2x"
-  smashed_var_bool <- str_detect(parts, "^[-. ^0-9]+[a-zA-Z]")
-  if(any(smashed_var_bool)){
-    places_to_break <- str_locate(parts[smashed_var_bool], "[a-zA-Z]")[,1]
-    for(k in seq_along(places_to_break)){
-      parts[smashed_var_bool][k] <- str_c(
-        str_sub(parts[smashed_var_bool][k], 1, places_to_break[k]-1),
-        "|",
-        str_sub(parts[smashed_var_bool][k], places_to_break[k])
-      )
-    }
-    parts <- unlist(str_split(parts, fixed("|")))
-  }
-  
-  # fix, e.g. "-y"
-  minus_var_bool <- str_detect(parts, "\\-[a-zA-Z]")
-  if(any(minus_var_bool)){
-    parts[minus_var_bool] <- str_c("-1 ", str_sub(parts[minus_var_bool], 2))
-    parts <- unlist(str_split(parts, " "))
-  }
-  
-  # collect numeric elements
-  parts_with_vars <- str_detect(parts, "[a-zA-Z]")
-  if(all(parts_with_vars)){
-    coef <- 1L
-  } else {
-    coef <- prod(
-      vapply(
-        as.list(parts[which(!parts_with_vars)]), 
-        function(.) eval(parse(text = .)),
-        double(1)
-      )  
-    ) # this multiplies even, e.g., 5^2
-  }
-  
-  # if only coefs are given, return
-  if(all(parts_with_vars == FALSE)) return(mpoly(list(c(coef = coef))))
-  
-  # parse variable exponents
-  var_parts <- parts[parts_with_vars]
-  var_parts_with_exps_bool <- str_detect(var_parts, fixed("^"))
-  var_parts[!var_parts_with_exps_bool] <- str_c(var_parts[!var_parts_with_exps_bool], "^1")
-  var_parts <- str_split(var_parts, fixed("^"))
-  vars <- vapply(var_parts, `[`, character(1), 1L)
-  exps <- as.integer(vapply(var_parts, `[`, character(1), 2L))
-  names(exps) <- vars
-  
-  # mpoly and return
-  mpoly(list(c(coef = coef, exps)))  
-}
+# parse_nonparenthetical_term <- function(string){
+#   
+#   # fix spaces around exponents "x ^ 2" -> "x^2"
+#   string <- str_replace_all(string, " *\\^ *", "^")
+#   
+#   # fix spaces around minuses "x  -  2" -> "x-2"
+#   string <- str_replace_all(string, " *- *", "-")
+#   
+#   # split based on spaces
+#   parts <- str_split(string, " ")[[1]]
+#   parts <- parts[nchar(parts) > 0] # for "2        -2"
+#   
+#   # if more than one negative provided error
+#   if(str_detect(str_sub(string, 2), fixed("-")))
+#     stop("Negative signs are only allowed at the beginning of terms.", call. = FALSE)
+#   
+#   # fix, e.g. "2x"
+#   smashed_var_bool <- str_detect(parts, "^[-. ^0-9]+[a-zA-Z]")
+#   if(any(smashed_var_bool)){
+#     places_to_break <- str_locate(parts[smashed_var_bool], "[a-zA-Z]")[,1]
+#     for(k in seq_along(places_to_break)){
+#       parts[smashed_var_bool][k] <- str_c(
+#         str_sub(parts[smashed_var_bool][k], 1, places_to_break[k]-1),
+#         "|",
+#         str_sub(parts[smashed_var_bool][k], places_to_break[k])
+#       )
+#     }
+#     parts <- unlist(str_split(parts, fixed("|")))
+#   }
+#   
+#   # fix, e.g. "-y"
+#   minus_var_bool <- str_detect(parts, "\\-[a-zA-Z]")
+#   if(any(minus_var_bool)){
+#     parts[minus_var_bool] <- str_c("-1 ", str_sub(parts[minus_var_bool], 2))
+#     parts <- unlist(str_split(parts, " "))
+#   }
+#   
+#   # collect numeric elements
+#   parts_with_vars <- str_detect(parts, "[a-zA-Z]")
+#   if(all(parts_with_vars)){
+#     coef <- 1L
+#   } else {
+#     coef <- prod(
+#       vapply(
+#         as.list(parts[which(!parts_with_vars)]), 
+#         function(.) eval(parse(text = .)),
+#         double(1)
+#       )  
+#     ) # this multiplies even, e.g., 5^2
+#   }
+#   
+#   # if only coefs are given, return
+#   if(all(parts_with_vars == FALSE)) return(mpoly(list(c(coef = coef))))
+#   
+#   # parse variable exponents
+#   var_parts <- parts[parts_with_vars]
+#   var_parts_with_exps_bool <- str_detect(var_parts, fixed("^"))
+#   var_parts[!var_parts_with_exps_bool] <- str_c(var_parts[!var_parts_with_exps_bool], "^1")
+#   var_parts <- str_split(var_parts, fixed("^"))
+#   vars <- vapply(var_parts, `[`, character(1), 1L)
+#   exps <- as.integer(vapply(var_parts, `[`, character(1), 2L))
+#   names(exps) <- vars
+#   
+#   # mpoly and return
+#   mpoly(list(c(coef = coef, exps)))  
+# }
 
 
 
@@ -475,58 +506,58 @@ parse_nonparenthetical_term <- function(string){
 # fix_term_joins("1e2 x")
 # fix_term_joins("-1-1-") # error
 # fix_term_joins("-1-1+") # error
-fix_term_joins <- function(string){
-  
-  # make sure last char is not a sign
-  if(str_detect(string, "[+-]$")) stop(sprintf("Term %s does not terminate.", string), call. = FALSE)
-  
-  # zero trick for leading symbol, e.g. "-1 + x" -> "0 + -1 + x"
-  if (str_detect(string, "^[+-]")) {
-    if (str_detect(string, "^[+-]{2,}")) stop(
-      sprintf("%s cannot start an expression.", str_extract(string, "^[+-]+")), 
-      call. = FALSE
-    )
-    string <- str_c("0 + ", string)
-  }
-  
-  # fix scientific notation
-  sciRegex <- "[0-9.]+e[+-]?[0-9]+"
-  while(str_detect(string, sciRegex)){
-    stringToReplace <- str_extract(string, sciRegex)
-    replacement <- format(as.numeric(stringToReplace))
-    string <- str_replace(string, sciRegex, replacement)
-  }
-  
-  # break string into pieces of terms and joins
-  terms <- str_extract_all(string, "[\\w^.,|\\[\\]]+")[[1]]
-  joins <- str_split(string, "[\\w^.,|\\[\\]]+")[[1]]
-  if(joins[1] == "") joins <- joins[-1]
-  if(joins[length(joins)] == "") joins <- joins[-length(joins)]
-  if(length(joins) == 0L) return(string)
-  
-  # fix joins
-  pureJoins <- str_replace_all(joins, "\\s", "")
-  pureJoins[pureJoins == ""] <- "|"
-  if(any(nchar(pureJoins) > 3)) stop("Arithmetic sign sequence of more than two detected.", call. = FALSE)
-  cleanJoinMap <- c(
-    "-" = " + -1 ",  "+" = " + ", "--" = " + ", 
-    "++" = " + ", "+-" = " + -1 ", "-+" = " + -1 ", "|" = " "
-  )
-  cleanedJoins <- unname(cleanJoinMap[pureJoins]) # cbind(joins, cleanedJoins)
-  
-  # reconstruct
-  n <- length(terms) + length(joins) # n always odd, first term always a \\w
-  temp <- character(n)
-  temp[seq.int(1L, n, 2L)] <- terms
-  temp[seq.int(2L, n-1L, 2L)] <- cleanedJoins    
-  string <- str_c(temp, collapse = "")
-  
-  # strip leading "0 + " if needed
-  if(str_sub(string, 1L, 4L) == "0 + ") string <- str_sub(string, 5L)
-  
-  # return
-  string
-}
+# fix_term_joins <- function(string){
+#   
+#   # make sure last char is not a sign
+#   if(str_detect(string, "[+-]$")) stop(sprintf("Term %s does not terminate.", string), call. = FALSE)
+#   
+#   # zero trick for leading symbol, e.g. "-1 + x" -> "0 + -1 + x"
+#   if (str_detect(string, "^[+-]")) {
+#     if (str_detect(string, "^[+-]{2,}")) stop(
+#       sprintf("%s cannot start an expression.", str_extract(string, "^[+-]+")), 
+#       call. = FALSE
+#     )
+#     string <- str_c("0 + ", string)
+#   }
+#   
+#   # fix scientific notation
+#   sciRegex <- "[0-9.]+e[+-]?[0-9]+"
+#   while(str_detect(string, sciRegex)){
+#     stringToReplace <- str_extract(string, sciRegex)
+#     replacement <- format(as.numeric(stringToReplace))
+#     string <- str_replace(string, sciRegex, replacement)
+#   }
+#   
+#   # break string into pieces of terms and joins
+#   terms <- str_extract_all(string, "[\\w^.,|\\[\\]]+")[[1]]
+#   joins <- str_split(string, "[\\w^.,|\\[\\]]+")[[1]]
+#   if(joins[1] == "") joins <- joins[-1]
+#   if(joins[length(joins)] == "") joins <- joins[-length(joins)]
+#   if(length(joins) == 0L) return(string)
+#   
+#   # fix joins
+#   pureJoins <- str_replace_all(joins, "\\s", "")
+#   pureJoins[pureJoins == ""] <- "|"
+#   if(any(nchar(pureJoins) > 3)) stop("Arithmetic sign sequence of more than two detected.", call. = FALSE)
+#   cleanJoinMap <- c(
+#     "-" = " + -1 ",  "+" = " + ", "--" = " + ", 
+#     "++" = " + ", "+-" = " + -1 ", "-+" = " + -1 ", "|" = " "
+#   )
+#   cleanedJoins <- unname(cleanJoinMap[pureJoins]) # cbind(joins, cleanedJoins)
+#   
+#   # reconstruct
+#   n <- length(terms) + length(joins) # n always odd, first term always a \\w
+#   temp <- character(n)
+#   temp[seq.int(1L, n, 2L)] <- terms
+#   temp[seq.int(2L, n-1L, 2L)] <- cleanedJoins    
+#   string <- str_c(temp, collapse = "")
+#   
+#   # strip leading "0 + " if needed
+#   if(str_sub(string, 1L, 4L) == "0 + ") string <- str_sub(string, 5L)
+#   
+#   # return
+#   string
+# }
 
 
 
@@ -546,31 +577,31 @@ fix_term_joins <- function(string){
 # string <- "-1 (x + y)+ 2 x (x + y) + 3 y"
 # string <- "2 (1 + x + (x - y))+ 2 x (x + y) + 3 y"
 # extract_polynomial_terms(string)
-extract_polynomial_terms <- function(string){
-  
-  # str_split(string, " *(?<!\\([\\w ]+)[+-] *")
-  
-  # run fix_term_joins on blanked strings to get protect parentheticals
-  blanked_string <- blank_parentheticals(string, "|")
-  piped_string   <- fix_term_joins(blanked_string)
-  
-  # change +"s to *"s for breaking later
-  # they distinguish polynomial terms
-  piped_string <- str_replace_all(piped_string, fixed("+"), "*")
-  
-  # unprotect
-  string_ndcs <- str_locate_all(blanked_string, "[|]+")[[1]]
-  piped_ndcs  <- str_locate_all(piped_string,   "[|]+")[[1]]  
-  if(nrow(string_ndcs) > 0){
-    for(k in 1:nrow(string_ndcs)){
-      str_sub(piped_string, piped_ndcs[k,1], piped_ndcs[k,2]) <- 
-        str_sub(string, string_ndcs[k,1], string_ndcs[k,2])  
-    }
-  }
-  
-  # split
-  str_split(piped_string, fixed("*"))[[1]]
-}
+# extract_polynomial_terms <- function(string){
+#   
+#   # str_split(string, " *(?<!\\([\\w ]+)[+-] *")
+#   
+#   # run fix_term_joins on blanked strings to get protect parentheticals
+#   blanked_string <- blank_parentheticals(string, "|")
+#   piped_string   <- fix_term_joins(blanked_string)
+#   
+#   # change +"s to *"s for breaking later
+#   # they distinguish polynomial terms
+#   piped_string <- str_replace_all(piped_string, fixed("+"), "*")
+#   
+#   # unprotect
+#   string_ndcs <- str_locate_all(blanked_string, "[|]+")[[1]]
+#   piped_ndcs  <- str_locate_all(piped_string,   "[|]+")[[1]]  
+#   if(nrow(string_ndcs) > 0){
+#     for(k in 1:nrow(string_ndcs)){
+#       str_sub(piped_string, piped_ndcs[k,1], piped_ndcs[k,2]) <- 
+#         str_sub(string, string_ndcs[k,1], string_ndcs[k,2])  
+#     }
+#   }
+#   
+#   # split
+#   str_split(piped_string, fixed("*"))[[1]]
+# }
 
 
 
@@ -602,11 +633,11 @@ extract_polynomial_terms <- function(string){
 # extract_leftmost_inner_parenthetical("((x + 5)^10+2)^2")
 # extract_leftmost_inner_parenthetical("((x + 5)^10+2)", contents_only = TRUE)
 # extract_leftmost_inner_parenthetical("(1 + (x + 5)^10+2)^2")
-extract_leftmost_inner_parenthetical <- function(string, contents_only = FALSE){
-  string <- str_extract(string, "\\([^()]*\\)(?:\\^[0-9]+)?")
-  if(!contents_only) return(string)
-  str_extract(string, "\\(.*\\)") ->.; str_sub(., 2L, -2L) 
-}
+# extract_leftmost_inner_parenthetical <- function(string, contents_only = FALSE){
+#   string <- str_extract(string, "\\([^()]*\\)(?:\\^[0-9]+)?")
+#   if(!contents_only) return(string)
+#   str_extract(string, "\\(.*\\)") ->.; str_sub(., 2L, -2L) 
+# }
 
 
 
@@ -615,22 +646,22 @@ extract_leftmost_inner_parenthetical <- function(string, contents_only = FALSE){
 # blank_parentheticals(" -1 1 x (3 x + -1 (7 + -1 2 x))^2 7 (x + 1) -3 ")
 # blank_parentheticals(" -1 1 x (3 x + -1 (7 + -1 2 x))^2 7 (x + 1) -3 ", "*")
 # blank_parentheticals(" -1 1 x (3 x + -1 (7 + -1 2 x))^2 7 (x + 1) -3 ", "_")
-blank_parentheticals <- function(string, char = "-"){
-  # " -1 1 x (3 x + -1 (7 + -1 2 x))^2 7 (x + 1) -3 " ->
-  # " -1 1 x ------------------------- 7 ------- -3 "
-  # this blanks parentheticals from the inside out
-  # inside parentheticals are done first
-  
-  while(contains_parenthetical_expression(string)){
-    bad <- extract_leftmost_inner_parenthetical(string)
-    string <- str_replace(
-      string, 
-      "\\([^()]*\\)(?:\\^[0-9]+)?", 
-      str_dup(char, nchar(bad))
-    )
-  }
-  string
-}
+# blank_parentheticals <- function(string, char = "-"){
+#   # " -1 1 x (3 x + -1 (7 + -1 2 x))^2 7 (x + 1) -3 " ->
+#   # " -1 1 x ------------------------- 7 ------- -3 "
+#   # this blanks parentheticals from the inside out
+#   # inside parentheticals are done first
+#   
+#   while(contains_parenthetical_expression(string)){
+#     bad <- extract_leftmost_inner_parenthetical(string)
+#     string <- str_replace(
+#       string, 
+#       "\\([^()]*\\)(?:\\^[0-9]+)?", 
+#       str_dup(char, nchar(bad))
+#     )
+#   }
+#   string
+# }
 
 
 
@@ -642,21 +673,21 @@ blank_parentheticals <- function(string, char = "-"){
 # extract_nonparenthetical_elements(string)
 # string <- " (x + y)^2  (x - y)(x) "
 # extract_nonparenthetical_elements(string)
-extract_nonparenthetical_elements <- function(string){
-  
-  # remove parenthetical stuff
-  parenthetical_regex <- "\\([-+*a-zA-Z0-9.^ ()]+\\)(\\^\\d+)?"
-  nonparem_elts <- str_remove_all(string, parenthetical_regex)
-  nonparem_elts <- str_replace_all(nonparem_elts, "\\s+", " ")
-  nonparem_elts <- str_trim(nonparem_elts)
-  
-  # parenthesize and return
-  if (nonparem_elts == "") {
-    ""
-  } else {
-    str_c("(", str_trim(nonparem_elts), ")")
-  }  
-}
+# extract_nonparenthetical_elements <- function(string){
+#   
+#   # remove parenthetical stuff
+#   parenthetical_regex <- "\\([-+*a-zA-Z0-9.^ ()]+\\)(\\^\\d+)?"
+#   nonparem_elts <- str_remove_all(string, parenthetical_regex)
+#   nonparem_elts <- str_replace_all(nonparem_elts, "\\s+", " ")
+#   nonparem_elts <- str_trim(nonparem_elts)
+#   
+#   # parenthesize and return
+#   if (nonparem_elts == "") {
+#     ""
+#   } else {
+#     str_c("(", str_trim(nonparem_elts), ")")
+#   }  
+# }
 
 # string <- " -3 (x + y)^2 4 (x - y)x 4 "
 # delete_nonparenthetical_elements(string)
@@ -664,12 +695,12 @@ extract_nonparenthetical_elements <- function(string){
 # delete_nonparenthetical_elements(string)
 # string <- ".2 (x)"
 # delete_nonparenthetical_elements(string)
-extract_parenthetical_elements <- function(string){
-  
-  parenthetical_regex <- "\\([-+*a-zA-Z0-9.^ ()]+\\)(\\^\\d+)?"
-  str_extract_all(string, parenthetical_regex)[[1]]
-  
-}
+# extract_parenthetical_elements <- function(string){
+#   
+#   parenthetical_regex <- "\\([-+*a-zA-Z0-9.^ ()]+\\)(\\^\\d+)?"
+#   str_extract_all(string, parenthetical_regex)[[1]]
+#   
+# }
 
 
 
@@ -679,13 +710,13 @@ extract_parenthetical_elements <- function(string){
 # term_parentheticals(string)
 # string <- ".2 (x)"
 # term_parentheticals(string) 
-term_parentheticals <- function(string){
-  
-  nonparens <- extract_nonparenthetical_elements(string)
-  parens    <-    extract_parenthetical_elements(string)
-  c(nonparens, parens)
-  
-}
+# term_parentheticals <- function(string){
+#   
+#   nonparens <- extract_nonparenthetical_elements(string)
+#   parens    <-    extract_parenthetical_elements(string)
+#   c(nonparens, parens)
+#   
+# }
 
 
 
@@ -696,9 +727,9 @@ term_parentheticals <- function(string){
 
 
 
-contains_parenthetical_expression <- function(string){ 
-  any(str_detect(string, fixed("(")))
-}
+# contains_parenthetical_expression <- function(string){ 
+#   any(str_detect(string, fixed("(")))
+# }
 
 
 
@@ -708,36 +739,36 @@ contains_parenthetical_expression <- function(string){
 # contains_nested_parenthetical_expression("((5+5))")
 # contains_nested_parenthetical_expression("x + (5 y) + 2")
 # contains_nested_parenthetical_expression("x + ((5 y) + 2)")
-contains_nested_parenthetical_expression <- function(string){
-  only_parentheses <- str_replace_all(string, "[^()]", "")
-  str_detect(only_parentheses, fixed("(("))
-}
+# contains_nested_parenthetical_expression <- function(string){
+#   only_parentheses <- str_replace_all(string, "[^()]", "")
+#   str_detect(only_parentheses, fixed("(("))
+# }
 
 
 
-unmatched_parentheses_stop <- function(string){
-  if(contains_parenthetical_expression(string)){
-    open_paren_count <- str_count(string, fixed("("))
-    closed_paren_count <- str_count(string, fixed(")"))
-    if (open_paren_count > closed_paren_count){
-      stop("Parenthetical error: excess ('s detected.", call. = FALSE)
-    } else if(open_paren_count < closed_paren_count) {
-      stop("Parenthetical error: excess )'s detected.", call. = FALSE)
-    }
-  }
-  invisible()
-}
+# unmatched_parentheses_stop <- function(string){
+#   if(contains_parenthetical_expression(string)){
+#     open_paren_count <- str_count(string, fixed("("))
+#     closed_paren_count <- str_count(string, fixed(")"))
+#     if (open_paren_count > closed_paren_count){
+#       stop("Parenthetical error: excess ('s detected.", call. = FALSE)
+#     } else if(open_paren_count < closed_paren_count) {
+#       stop("Parenthetical error: excess )'s detected.", call. = FALSE)
+#     }
+#   }
+#   invisible()
+# }
 
 
-empty_parenthetical_stop <- function(string) {
-  if (str_detect(string, "\\( *\\)")) {
-    stop("Expression contains empty parenthetical.", call. = FALSE)
-  }
-}
+# empty_parenthetical_stop <- function(string) {
+#   if (str_detect(string, "\\( *\\)")) {
+#     stop("Expression contains empty parenthetical.", call. = FALSE)
+#   }
+# }
 
 
 
-str_rev <- function(string) str_c(rev.default(str_split(string, "")[[1]]), collapse = "")
+# str_rev <- function(string) str_c(rev.default(str_split(string, "")[[1]]), collapse = "")
 
 
 
diff --git a/R/mpolyList.R b/R/mpolyList.R
index 5e577e4..0fc059b 100644
--- a/R/mpolyList.R
+++ b/R/mpolyList.R
@@ -33,7 +33,7 @@ mpolyList <- function(...){
 
   out <- lapply(arguments, eval, parent.frame(1))
 
-  if(is.mpoly(out)) out <- list(out)
+  # if(is.mpoly(out)) out <- list(out)
   
   if(!all(vapply(out, is.mpoly, logical(1)))){
   	stop("Each argument must be of class mpoly.", call. = FALSE)
diff --git a/R/mpolyListArithmetic.R b/R/mpolyListArithmetic.R
index fd0501a..cb05a21 100644
--- a/R/mpolyListArithmetic.R
+++ b/R/mpolyListArithmetic.R
@@ -8,8 +8,8 @@
 #' @name mpolyListArithmetic
 #' @examples
 #'
-#' ( ms1 <- mp( c('x + 1', 'x + 2') ) )
-#' ( ms2 <- mp( c('x + 1', 'y + 2') ) )
+#' ( ms1 <- mp( c("x", 'y') ) )
+#' ( ms2 <- mp( c("y", '2 x^2') ) )
 #' ms1 + ms2
 #' ms1 - ms2
 #' ms1 * ms2
@@ -26,43 +26,52 @@ NULL
 #' @export
 `+.mpolyList` <- function(e1, e2){
   
-  ## argument check
-  if(is.numeric(e1) && length(e1) == 1) e1 <- mpoly(list(c(coef = e1)))
-  if(is.mpoly(e1)) e1 <- mpolyList(e1)
-  
-  if(is.numeric(e2) && length(e2) == 1) e2 <- mpoly(list(c(coef = e2)))
-  if(is.mpoly(e2)) e2 <- mpolyList(e2)
-  
-  if(!is.mpolyList(e1) || !is.mpolyList(e2)){
-    stop('e1 and e2 must be of class mpolyList.', call. = FALSE)
+  # argument check
+  if (is.numeric(e1)) {
+    e1 <- structure(
+      lapply(e1, function(.) structure(list(c("coef" = .)), class = "mpoly")),
+      class = "mpolyList"
+    )
   }
+  # if (is.mpoly(e1)) e1 <- mpolyList(e1)
   
-  if(length(e1) != length(e2)){
-    stop('e1 and e2 must have equal length.', call. = FALSE)
+  if (is.numeric(e2)) {
+    e2 <- structure(
+      lapply(e2, function(.) structure(list(c("coef" = .)), class = "mpoly")),
+      class = "mpolyList"
+    )
   }
+  # if (is.mpoly(e2)) e2 <- mpolyList(e2)
   
+  stopifnot(is.mpolyList(e1))
+  stopifnot(is.mpolyList(e2))
   
-  ## determine length, flatten, and make indices on which to add
-  n <- length(e1)
-  
-  flatList <- unlist(list(
-    unclass(e1),
-    unclass(e2)
-  ), recursive = FALSE)
+  # fix lengths
+  if (length(e1) == 1L && length(e2) != 1L) {
+    e1 <- structure(
+      replicate(length(e2), e1[[1L]], simplify = FALSE),
+      class = "mpolyList"
+    )
+  }
   
-  ndcs2add <- split(cbind(1:n, (n+1):(2*n)), 1:n)
+  if (length(e1) != 1L && length(e2) == 1L) {
+    e2 <- structure(
+      replicate(length(e1), e2[[1L]], simplify = FALSE),
+      class = "mpolyList"
+    )
+  }
   
+  if(length(e1) != length(e2)) stop("e1 and e2 must have equal length.", call. = FALSE)
+
   
-  ## sum
-  out <- lapply(ndcs2add, function(v){
-    Reduce('+', flatList[v])
-  })
-  out <- unname(out)
+  # template outcome
+  out <- vector(length = length(e1), mode = "list")
   
+  # compute sums
+  for (k in seq_along(out)) out[[k]] <- e1[[k]] + e2[[k]]
   
-  ## caste and return
-  class(out) <- 'mpolyList'
-  out
+  # return
+  structure(out, class = "mpolyList")
 }	
 
 
@@ -78,19 +87,7 @@ NULL
 
 #' @rdname mpolyListArithmetic
 #' @export
-`-.mpolyList` <- function(e1, e2){
-  ## change coefficient signs in e2
-  e2 <- lapply(e2, unclass)
-  e2 <- lapply(e2, function(l){
-    lapply(l, function(v){
-      v['coef'] <- -v['coef']
-      v
-    })	
-  })
-  class(e2) <- 'mpolyList'
-    
-  e1 + e2
-}
+`-.mpolyList` <- function(e1, e2) e1 + -1*e2
 
 
 
@@ -106,43 +103,51 @@ NULL
 #' @export
 `*.mpolyList` <- function(e1, e2){
   
-  ## argument check
-  
-  if(is.numeric(e1) && length(e1) == 1) e1 <- mpoly(list(c(coef = e1)))  
-  if(is.mpoly(e1)) e1 <- mpolyList(e1)
+  # argument check
+  if (is.numeric(e1)) {
+    e1 <- structure(
+      lapply(e1, function(.) structure(list(c("coef" = .)), class = "mpoly")),
+      class = "mpolyList"
+    )
+  }
+  # if (is.mpoly(e1)) e1 <- mpolyList(e1)
   
-  if(is.numeric(e2) && length(e2) == 1) e2 <- mpoly(list(c(coef = e2)))  
-  if(is.mpoly(e2)) e2 <- mpolyList(e2)
+  if (is.numeric(e2)) {
+    e2 <- structure(
+      lapply(e2, function(.) structure(list(c("coef" = .)), class = "mpoly")),
+      class = "mpolyList"
+    )
+  }
+  # if (is.mpoly(e2)) e2 <- mpolyList(e2)
   
   stopifnot(is.mpolyList(e1))
   stopifnot(is.mpolyList(e2))
   
-  if(length(e1) != length(e2)) stop('e1 and e2 must have equal length.', call. = FALSE)
-  
-  
-  
-  ## determine length, flatten, and make indices on which to multiply
-  n <- length(e1)
-  
-  flatList <- unlist(list(
-    unclass(e1),
-    unclass(e2)
-  ), recursive = FALSE)
-  
-  ndcs2add <- split(cbind(1:n, (n+1):(2*n)), 1:n)
-  
-  
+  # fix lengths
+  if (length(e1) == 1L && length(e2) != 1L) {
+    e1 <- structure(
+      replicate(length(e2), e1[[1L]], simplify = FALSE),
+      class = "mpolyList"
+    )
+  }
   
+  if (length(e1) != 1L && length(e2) == 1L) {
+    e2 <- structure(
+      replicate(length(e1), e2[[1L]], simplify = FALSE),
+      class = "mpolyList"
+    )
+  }
   
-  ## multiply
-  out <- lapply(ndcs2add, function(v) Reduce('*', flatList[v]))
-  out <- unname(out)
+  if(length(e1) != length(e2)) stop("e1 and e2 must have equal length.", call. = FALSE)
   
+  # template outcome
+  out <- vector(length = length(e1), mode = "list")
   
+  # compute products
+  for (k in seq_along(out)) out[[k]] <- e1[[k]] * e2[[k]]
   
-  ## caste and return
-  class(out) <- 'mpolyList'
-  out
+  # return
+  structure(out, class = "mpolyList")
 }	
 
 
diff --git a/R/reorder.mpoly.R b/R/reorder.mpoly.R
index 00ec91e..5665f45 100644
--- a/R/reorder.mpoly.R
+++ b/R/reorder.mpoly.R
@@ -49,10 +49,9 @@ reorder.mpoly <- function(x, varorder = vars(x), order = "lex", ...){
   
   if(!missing(varorder)){
     if(!all(vars %in% varorder)){
-      error <- paste(
-        "if specified, varorder must contain all computed vars - ",
-        paste(vars, collapse = ", "),
-        sep = ""
+      error <- stri_c(
+        "If specified, varorder must contain all computed vars - ",
+        paste(vars, collapse = ", ")
       )
       stop(error, call. = FALSE)
     }
@@ -69,10 +68,9 @@ reorder.mpoly <- function(x, varorder = vars(x), order = "lex", ...){
   }  
   
   if(missing(varorder) && !missing(order)){
-    message <- paste(
-      "using variable ordering - ",
-      paste(vars, collapse = ", "),
-      sep = ""
+    message <- stri_c(
+      "Using variable ordering - ",
+      paste(vars, collapse = ", ")
     )
     message(message)
   }
diff --git a/man/as.function.mpoly.Rd b/man/as.function.mpoly.Rd
index 7010c81..ca5e43d 100644
--- a/man/as.function.mpoly.Rd
+++ b/man/as.function.mpoly.Rd
@@ -61,6 +61,14 @@ p <- mp("1 2 3 4")
 f <- as.function(p)
 f(10) # -> 24
 
+bernstein(1, 2)
+s <- seq(0, 1, .01)
+as.function(bernstein(1, 2))(s)
+plot(
+  s,
+  as.function(bernstein(1, 2))(s)
+)
+
 
 as.function(mp("x + xx"))
 as.function(mp("x + xx"), squeeze = FALSE)
diff --git a/man/mpolyListArithmetic.Rd b/man/mpolyListArithmetic.Rd
index 8da59e8..4495a1c 100644
--- a/man/mpolyListArithmetic.Rd
+++ b/man/mpolyListArithmetic.Rd
@@ -26,8 +26,8 @@ Element-wise arithmetic with vectors of multivariate polynomials.
 }
 \examples{
 
-( ms1 <- mp( c('x + 1', 'x + 2') ) )
-( ms2 <- mp( c('x + 1', 'y + 2') ) )
+( ms1 <- mp( c("x", 'y') ) )
+( ms2 <- mp( c("y", '2 x^2') ) )
 ms1 + ms2
 ms1 - ms2
 ms1 * ms2
diff --git a/tests/testthat/test-as-function-mpoly.R b/tests/testthat/test-as-function-mpoly.R
new file mode 100644
index 0000000..981928d
--- /dev/null
+++ b/tests/testthat/test-as-function-mpoly.R
@@ -0,0 +1,69 @@
+context("as.function.mpoly()")
+
+
+
+test_that("univariate", {
+  
+  f <- as.function(mp("(3 - x)^2"), silent = TRUE)
+  
+  x <- 2
+  expect_equal(f(x), (3 - x)^2)
+  
+  x <- 5
+  expect_equal(f(x), (3 - x)^2)
+  
+  f <- as.function(mp("(3 - x)^2"), silent = TRUE)
+  
+  x <- 2:10
+  expect_equal(f(x), (3 - x)^2)
+
+})
+
+
+
+test_that("vector = TRUE", {
+  
+  f <- as.function(mp("3 x + y^2"), vector = TRUE, silent = TRUE)
+  
+  expect_equal(f(c(1,2)), 3*(1) + (2)^2)
+  
+})
+
+
+test_that("vector = FALSE", {
+  
+  f <- as.function(mp("3 x + y^2"), vector = FALSE, silent = TRUE)
+  
+  expect_equal(f(1, 2), 3*(1) + (2)^2)
+  
+})
+
+
+test_that("constant mpoly", {
+  
+  f <- as.function(mp("3"), silent = TRUE)
+  
+  expect_equal(f(1), 3)
+  expect_equal(f(1:4), rep(3, 4))
+  
+})
+
+
+test_that("bernstein", {
+  
+  s <- seq(0, 1, .01)
+  
+  expect_equal(
+    as.function(bernstein(1, 2))(s),
+    2*s - 2*s^2
+  )
+  
+  s <- seq(0, 1, .01)
+  
+  expect_equal(
+    as.function(bernstein(1, 2))(-3:3),
+    2*(-3:3) - 2*(-3:3)^2
+  )
+  
+})
+
diff --git a/tests/testthat/test-as-function.R b/tests/testthat/test-as-function.R
deleted file mode 100644
index e29a5ed..0000000
--- a/tests/testthat/test-as-function.R
+++ /dev/null
@@ -1,27 +0,0 @@
-context("as.function.mpoly() works properly")
-
-
-
-test_that("as.function.mpoly() basic functionality works", {
-  
-  f <- as.function(mp("(3 - x)^2"), silent = TRUE)
-  
-  x <- 2
-  expect_equal(f(x), (3 - x)^2)
-  
-  x <- 5
-  expect_equal(f(x), (3 - x)^2)
-
-})
-
-
-
-test_that("as.function.mpoly() creates a vectorized function", {
-  
-  f <- as.function(mp("(3 - x)^2"), silent = TRUE)
-  
-  x <- 2:10
-  expect_equal(f(x), (3 - x)^2)
-  
-})
-
diff --git a/tests/testthat/test-as-mpoly.R b/tests/testthat/test-as-mpoly.R
new file mode 100644
index 0000000..baa14af
--- /dev/null
+++ b/tests/testthat/test-as-mpoly.R
@@ -0,0 +1,135 @@
+context("as.mpoly()")
+
+test_that("numeric", {
+  
+  expect_equal(
+    as.mpoly(3),
+    structure(list(c("coef" = 3)), class = "mpoly")
+  )
+  
+  expect_equal(
+    as.mpoly(1:3),
+    structure(
+      list(
+        c("coef" = 1),
+        c("x" = 1, "coef" = 2),
+        c("x" = 2, "coef" = 3)
+      ), 
+      class = "mpoly"
+    )
+  )
+  
+})
+
+
+
+test_that("polynomial", {
+  
+  expect_equal(
+    as.mpoly(hermite.h.polynomials(3)[[4]]),
+    structure(
+      list(
+        c("x" = 1, "coef" = -12),
+        c("x" = 3, "coef" = 8)
+      ), 
+      class = "mpoly"
+    )
+  )
+  
+})
+
+
+
+
+test_that("lm", {
+  
+  expect_equal(
+    as.mpoly(lm(y ~ x, data = data.frame("x" = c(0,1), "y" = c(1,3)))),
+    structure(
+      list(
+        c("x" = 1, "coef" = 2),
+        c("coef" = 1)
+      ), 
+      class = "mpoly"
+    )
+  )
+  
+  
+  expect_equal(
+    as.mpoly(
+      lm(
+        y ~ poly(x, 1, raw = TRUE), 
+        data = data.frame("x" = c(0,1), "y" = c(1,3))
+      )
+    ),
+    structure(
+      list(
+        c("x" = 1, "coef" = 2),
+        c("coef" = 1)
+      ), 
+      class = "mpoly"
+    )
+  )
+  
+  
+  expect_equal(
+    as.mpoly(
+      lm(
+        y ~ poly(x, 2, raw = TRUE), 
+        data = data.frame("x" = c(-1,0,1), "y" = c(0,-1,0))
+      )
+    ),
+    structure(
+      list(
+        c("x" = 1, "coef" = 0),
+        c("x" = 2, "coef" = 1),
+        c("coef" = -1)
+      ), 
+      class = "mpoly"
+    )
+  )
+  
+  
+  expect_equal(
+    as.mpoly(
+      lm(
+        y ~ x + I(x^2), 
+        data = data.frame("x" = c(-1,0,1), "y" = c(0,-1,0))
+      )
+    ),
+    structure(
+      list(
+        c("x" = 1, "coef" = 0),
+        c("x" = 2, "coef" = 1),
+        c("coef" = -1)
+      ), 
+      class = "mpoly"
+    )
+  )
+  
+  
+  expect_error(
+    as.mpoly(
+      lm(
+        y ~ poly(x, 2, raw = FALSE), 
+        data = data.frame("x" = c(-1,0,1), "y" = c(0,-1,0))
+      )
+    ),
+    "poly() statements currently must contain raw = TRUE.",
+    fixed = TRUE
+  )
+  
+})
+
+
+
+
+test_that("errors", {
+  
+  expect_error(
+    as.mpoly("x"),
+    "objects of class character not supported by as.mpoly().",
+    fixed = TRUE
+  )
+  
+})
\ No newline at end of file
diff --git a/tests/testthat/test-burst.R b/tests/testthat/test-burst.R
index e960895..e4e188b 100644
--- a/tests/testthat/test-burst.R
+++ b/tests/testthat/test-burst.R
@@ -1,4 +1,4 @@
-context("burst() is working properly")
+context("burst()")
 
 
 test_that("basic burst works", {
diff --git a/tests/testthat/test-components.R b/tests/testthat/test-components.R
index c17a726..bcf7865 100644
--- a/tests/testthat/test-components.R
+++ b/tests/testthat/test-components.R
@@ -1,4 +1,4 @@
-context("Polynomial component extraction works properly")
+context("Polynomial component extraction")
 
 
 test_that("[.mpoly", {
diff --git a/tests/testthat/test-eq_mp.R b/tests/testthat/test-eq_mp.R
index 9ffbf38..1044b22 100644
--- a/tests/testthat/test-eq_mp.R
+++ b/tests/testthat/test-eq_mp.R
@@ -1,4 +1,4 @@
-context("eq_mp() is working properly")
+context("eq_mp()")
 
 
 test_that("eq_mp() properly deals with spaces", {
diff --git a/tests/testthat/test-homogenize.R b/tests/testthat/test-homogenize.R
index b02dc87..73aa903 100644
--- a/tests/testthat/test-homogenize.R
+++ b/tests/testthat/test-homogenize.R
@@ -1,4 +1,4 @@
-context("homogenize")
+context("homogenizing functions")
 
 
 test_that("is.homogeneous()", {
diff --git a/tests/testthat/test-mp.r b/tests/testthat/test-mp.r
index 668cf73..b306162 100644
--- a/tests/testthat/test-mp.r
+++ b/tests/testthat/test-mp.r
@@ -1,4 +1,4 @@
-context("mpoly() is working properly")
+context("mpoly()")
 
 test_that("mpoly() flags non-list arguments",{
   expect_error(mpoly(1:5), "input to mpoly must be a list.")
@@ -86,6 +86,29 @@ test_that("mp() parses character vectors properly", {
 
 
 
+
+
+
+test_that("mp() obeys varorder", {
+  
+  expect_equal(
+    mp("3 y^2 x^3", varorder = c("x", "y")),
+    structure(list(c(x = 3, y = 2, coef = 3)), class = "mpoly")
+  )
+  
+})
+
+
+
+
+
+
+
+
+
+
+
+
 # 
 # 
 # 
diff --git a/tests/testthat/test-mpolyArithmetic.R b/tests/testthat/test-mpolyArithmetic.R
index 1caee3d..7f02c84 100644
--- a/tests/testthat/test-mpolyArithmetic.R
+++ b/tests/testthat/test-mpolyArithmetic.R
@@ -1,4 +1,4 @@
-context("Arithmetic with constants works properly")
+context("mpoly arithmetic")
 
 test_that("Addition works", {
   
diff --git a/tests/testthat/test-mpolylist-arithmetic.R b/tests/testthat/test-mpolylist-arithmetic.R
new file mode 100644
index 0000000..787ae19
--- /dev/null
+++ b/tests/testthat/test-mpolylist-arithmetic.R
@@ -0,0 +1,187 @@
+context("mpolyList arithmetic")
+
+
+test_that("addition", {
+  
+  expect_equal(
+    mp(c("x", "y")) + mp(c("y^2", "x")),
+    structure(
+      list(
+        structure(
+          list(
+            c(x = 1, coef = 1), 
+            c(y = 2, coef = 1)
+          ), 
+          class = "mpoly"
+        ), 
+        structure(
+          list(
+            c(y = 1, coef = 1), 
+            c(x = 1, coef = 1)
+          ), 
+          class = "mpoly")
+        ), 
+      class = "mpolyList"
+    )
+  )
+  
+  # expect_equal(
+  #   mp("2") + mp(c("x", "y")),
+  #   mp(c("x + 2", "y + 2"))
+  # )
+  
+  expect_equal(
+    2 + mp(c("x", "y")),
+    mp(c("x + 2", "y + 2"))
+  )
+  
+  # expect_equal(
+  #   mp(c("x", "y")) + mp("2"),
+  #   mp(c("x + 2", "y + 2"))
+  # )
+  
+  expect_equal(
+    mp(c("x", "y")) + 2,
+    mp(c("x + 2", "y + 2"))
+  )
+  
+  expect_equal(
+    c(2, -3) + mp(c("x", "y")),
+    mp(c("x + 2", "y - 3"))
+  )
+  
+  expect_equal(
+    mp(c("x", "y")) + c(2, -3),
+    mp(c("x + 2", "y - 3"))
+  )
+  
+  expect_error(
+    c(2, -3) + mp(c("x", "y", "z")),
+    "e1 and e2 must have equal length.",
+    fixed = TRUE
+  )
+  
+  expect_error(
+    mp(c("x", "y", "z")) + c(2, -3),
+    "e1 and e2 must have equal length.",
+    fixed = TRUE
+  )
+    
+  
+})
+
+
+
+
+
+test_that("subtraction", {
+  
+  expect_equal(
+    mp(c("x", "y")) - mp(c("x", "x^2")),
+    structure(
+      list(
+        structure(list(c(coef = 0)), class = "mpoly"), 
+        structure(
+          list(
+            c(y = 1, coef = 1), 
+            c(x = 2, coef = -1)
+          ), 
+          class = "mpoly")
+        ), 
+      class = "mpolyList"
+    )
+  )
+  
+  # expect_equal(
+  #   mp("2") - mp(c("x", "y")),
+  #   mp(c("-1 x + 2", "-1 y + 2"))
+  # )
+  
+  expect_equal(
+    2 - mp(c("x", "y")),
+    mp(c("-1 x + 2", "-1 y + 2"))
+  )
+  
+  # expect_equal(
+  #   mp(c("x", "y")) - mp("2"),
+  #   mp(c("x - 2", "y - 2"))
+  # )
+  
+  expect_equal(
+    mp(c("x", "y")) - 2,
+    mp(c("x - 2", "y - 2"))
+  )
+  
+  expect_equal(
+    c(2, -3) - mp(c("x", "y")),
+    mp(c("-1 x + 2", "-1 y - 3"))
+  )
+  
+  expect_equal(
+    mp(c("x", "y")) - c(2, -3),
+    mp(c("x - 2", "y + 3"))
+  )
+  
+  expect_error(
+    c(2, -3) - mp(c("x", "y", "z")),
+    "e1 and e2 must have equal length.",
+    fixed = TRUE
+  )
+  
+  expect_error(
+    mp(c("x", "y", "z")) - c(2, -3),
+    "e1 and e2 must have equal length.",
+    fixed = TRUE
+  )
+  
+})
+
+
+
+
+
+test_that("multiplication", {
+  
+  expect_equal(
+    mp(c("x", "y")) * mp(c("x", "x^2")),
+    structure(
+      list(
+        structure(list(c(x = 2, coef = 1)), class = "mpoly"), 
+        structure(list(c(y = 1, x = 2, coef = 1)), class = "mpoly")
+      ), 
+      class = "mpolyList"
+    )
+  )
+  
+  # expect_equal(
+  #   mp("2") * mp(c("x", "y")),
+  #   mp(c("2 x", "2 y"))
+  # )
+  
+  expect_equal(
+    2 * mp(c("x", "y")),
+    mp(c("2 x", "2 y"))
+  )
+  
+  # expect_equal(
+  #   mp(c("x", "y")) * mp("2"),
+  #   mp(c("2 x", "2 y"))
+  # )
+  
+  expect_equal(
+    mp(c("x", "y")) * 2,
+    mp(c("2 x", "2 y"))
+  )
+  
+  expect_equal(
+    c(2, -3) * mp(c("x", "y")),
+    mp(c("2 x", "-3 y"))
+  )
+  
+  expect_error(
+    c(2, -3) * mp(c("x", "y", "z")),
+    "e1 and e2 must have equal length.",
+    fixed = TRUE
+  )
+  
+})
diff --git a/tests/testthat/test-permutations.R b/tests/testthat/test-permutations.R
index 726db1c..5532ba6 100644
--- a/tests/testthat/test-permutations.R
+++ b/tests/testthat/test-permutations.R
@@ -1,4 +1,4 @@
-context("permutations() is working properly")
+context("permutations()")
 
 test_that("permutations of a single integer works", {
   
diff --git a/tests/testthat/test-plug.R b/tests/testthat/test-plug.R
index 61022a7..68022db 100644
--- a/tests/testthat/test-plug.R
+++ b/tests/testthat/test-plug.R
@@ -1,4 +1,4 @@
-context("plug() works properly")
+context("plug()")
 
 
 test_that("plug() inserts numbers properly", {
diff --git a/tests/testthat/test-predicates.R b/tests/testthat/test-predicates.R
index fe93511..47127b4 100644
--- a/tests/testthat/test-predicates.R
+++ b/tests/testthat/test-predicates.R
@@ -1,4 +1,4 @@
-context("predicate functions are working properly")
+context("predicate functions")
 
 
 test_that("is.constant() works - single", {
diff --git a/tests/testthat/test-print-mpoly.R b/tests/testthat/test-print-mpoly.R
index 7f87fdf..7f50692 100644
--- a/tests/testthat/test-print-mpoly.R
+++ b/tests/testthat/test-print-mpoly.R
@@ -1,4 +1,4 @@
-context("print_term() works properly")
+context("print_term()")
 
 
 
@@ -120,3 +120,124 @@ test_that("print_term() prints multivariate terms properly", {
 
 
 
+
+
+
+
+context("print.mpoly()")
+
+
+test_that("constants", {
+  
+  test_coef <- function(number) {
+    p <- structure(list(c("coef" = number)), class = "mpoly")
+    expect_equal(
+      print.mpoly(p, silent = TRUE),
+      as.character(number)
+    )
+  }
+  
+  test_coef(0)
+  test_coef(1)
+  test_coef(-1)
+  test_coef(1.1)
+  test_coef(-1.1) 
+  test_coef(1e-3)
+  test_coef(-1e-3)
+  
+})
+
+
+
+
+
+test_that("sums: no coefs", {
+  
+  p <- structure(
+    list(
+      c("x" = 1, "coef" = 1),
+      c("y" = 2, "coef" = 1)
+    ), 
+    class = "mpoly"
+  )
+  
+  expect_equal(
+    print.mpoly(p, silent = TRUE),
+    "x  +  y^2"
+  )
+  
+})
+
+
+
+test_that("sums: + coefs", {
+  
+  p <- structure(
+    list(
+      c("x" = 2, "coef" = 2),
+      c("y" = 3, "coef" = 3)
+    ), 
+    class = "mpoly"
+  )
+  
+  expect_equal(
+    print.mpoly(p, silent = TRUE),
+    "2 x^2  +  3 y^3"
+  )
+  
+})
+
+
+test_that("subtraction: ordinary", {
+  
+  p <- structure(
+    list(
+      c("x" = 1, "coef" = 2),
+      c("y" = 3, "coef" = -3)
+    ), 
+    class = "mpoly"
+  )
+  
+  expect_equal(
+    print.mpoly(p, silent = TRUE),
+    "2 x  -  3 y^3"
+  )
+  
+})
+
+
+test_that("subtraction: 2 x - y^2", {
+  
+  p <- structure(
+    list(
+      c("x" = 1, "coef" = 2),
+      c("y" = 2, "coef" = -1)
+    ), 
+    class = "mpoly"
+  )
+  
+  expect_equal(
+    print.mpoly(p, silent = TRUE),
+    "2 x  -  y^2"
+  )
+  
+})
+
+
+
+
+
+test_that("errors", {
+  
+  expect_error(
+    print.mpoly(
+      structure(list(c("x" = 1, "y" = 2, "coef" = 1)), class = "mpoly"),
+      varorder = "x",
+      silent = TRUE
+    ),
+    "if specified, varorder must contain all computed vars - x, y",
+    fixed = TRUE
+  )
+  
+})
+
diff --git a/tests/testthat/test-print-mpolylist.R b/tests/testthat/test-print-mpolylist.R
index 03af98e..dd0e60c 100644
--- a/tests/testthat/test-print-mpolylist.R
+++ b/tests/testthat/test-print-mpolylist.R
@@ -1,4 +1,4 @@
-context("print.mpolyList() works properly")
+context("print.mpolyList()")
 
 
 test_that("basic print.mpolyList() functionality", {
diff --git a/tests/testthat/test-reorder-mpoly.R b/tests/testthat/test-reorder-mpoly.R
index da84854..ad781a3 100644
--- a/tests/testthat/test-reorder-mpoly.R
+++ b/tests/testthat/test-reorder-mpoly.R
@@ -1,4 +1,4 @@
-context("reorder.mpoly() is working properly")
+context("reorder.mpoly()")
 
 
 # examples from cox little and o'shea, 3e, p.59
@@ -57,3 +57,30 @@ test_that("grlex works", {
   )
   
 })
+
+
+
+test_that("message when no order specified but varorder not", {
+  
+  expect_message(
+    reorder.mpoly(mp("x + y^2"), order = "glex"),
+    "Using variable ordering - x, y",
+    fixed = TRUE
+  )
+  
+})
+
+
+
+test_that("all vars in varorder", {
+  
+  expect_error(
+    reorder.mpoly(mp("x + y^2"), varorder = "x"),
+    "If specified, varorder must contain all computed vars - x, y",
+    fixed = TRUE
+  )
+  
+})
+
+
+
diff --git a/tests/testthat/test-solve-unipoly.R b/tests/testthat/test-solve-unipoly.R
index 01b1adc..ab715c1 100644
--- a/tests/testthat/test-solve-unipoly.R
+++ b/tests/testthat/test-solve-unipoly.R
@@ -1,4 +1,4 @@
-context("solve_unipoly() is working properly")
+context("solve_unipoly()")
 
 
 test_that("solve_unipoly() works with real roots", {
diff --git a/tests/testthat/test-tuples.R b/tests/testthat/test-tuples.R
index b380620..71a57c6 100644
--- a/tests/testthat/test-tuples.R
+++ b/tests/testthat/test-tuples.R
@@ -1,4 +1,4 @@
-context("tuples() works properly")
+context("tuples()")
 
 
 test_that("basic tuples() functionality works", {