diff --git a/R/HeumanLambda.R b/R/HeumanLambda.R
index b6bdee8..52ca385 100644
--- a/R/HeumanLambda.R
+++ b/R/HeumanLambda.R
@@ -1,12 +1,12 @@
 #' @title Heuman Lambda function
 #' @description Evaluates the Heuman Lambda function.
 #'
-#' @param phi Jacobi amplitude, a complex number
-#' @param m parameter, a complex number
+#' @param phi Jacobi amplitude, a complex number/vector
+#' @param m parameter, a complex number/vector
 #' @param minerror the bound on the relative error passed to
 #'   \code{\link{elliptic_F}} and \code{\link{elliptic_Z}}
 #'
-#' @return A complex number.
+#' @return A complex number or vector.
 #' @export
 Lambda0 <- function(phi, m, minerror = 1e-14) {
   if(m == 0) {
diff --git a/R/elliptic_E.R b/R/elliptic_E.R
index 59d7fec..9d8f6d0 100644
--- a/R/elliptic_E.R
+++ b/R/elliptic_E.R
@@ -1,13 +1,13 @@
 #' Incomplete elliptic integral of the second kind
 #' @description Evaluate the incomplete elliptic integral of the second kind.
 #'
-#' @param phi amplitude, real or complex number
-#' @param m parameter, real or complex number
+#' @param phi amplitude, real or complex number/vector
+#' @param m parameter, real or complex number/vector
 #' @param minerror the bound on the relative error passed to
 #'   \code{\link{Carlson_RF}} and \code{\link{Carlson_RD}}
 #'
-#' @return A complex number, the value of the incomplete elliptic integral
-#'   \ifelse{html}{\out{E(&phi;,m)}}{\eqn{E(\phi,m)}{E(phi,m)}}.
+#' @return A complex number or vector, the value(s) of the incomplete elliptic
+#'   integral \ifelse{html}{\out{E(&phi;,m)}}{\eqn{E(\phi,m)}{E(phi,m)}}.
 #' @export
 #'
 #' @examples elliptic_E(1, 0.2)
diff --git a/R/elliptic_PI.R b/R/elliptic_PI.R
index 6420608..9a2ffe8 100644
--- a/R/elliptic_PI.R
+++ b/R/elliptic_PI.R
@@ -1,28 +1,36 @@
 #' Incomplete elliptic integral of the third kind
 #' @description Evaluate the incomplete elliptic integral of the third kind.
 #'
-#' @param phi amplitude, real or complex number
-#' @param n characteristic, real or complex number
-#' @param m parameter, real or complex number
+#' @param phi amplitude, real or complex number/vector
+#' @param n characteristic, real or complex number/vector
+#' @param m parameter, real or complex number/vector
 #' @param minerror the bound on the relative error passed to
 #'   \code{\link{Carlson_RF}} and \code{\link{Carlson_RJ}}
 #'
-#' @return A complex number, the value of the incomplete elliptic integral
+#' @return A complex number or vector, the value(s) of the incomplete elliptic
+#'   integral
 #'   \ifelse{html}{\out{&Pi;(&phi;,n,m)}}{\eqn{\Pi(\phi,n,m)}{PI(phi,n,m)}}.
 #' @export
 #'
 #' @examples elliptic_PI(1, 0.8, 0.2)
 #' gsl::ellint_P(1, sqrt(0.2), -0.8)
-elliptic_PI <- function(phi, n, m, minerror = 1e-15){
+elliptic_PI <- function(phi, n, m, minerror = 1e-15) {
   phi <- as.complex(phi)
   n   <- as.complex(n)
   m   <- as.complex(m)
-  if(length(phi) == 1L) {
-    phi <- rep(phi, length(m))
-  } else if(length(m) == 1L) {
-    m <- rep(m, length(phi))
-  } else if(length(phi) != length(m)) {
-    stop("Incompatible lengths of `phi` and `m`.")
+  lgths <- c(length(phi), length(n), length(m))
+  L <- max(lgths)
+  stopifnot(all(lgths %in% c(1L, L)))
+  if(L != 1L && any(lgths == 1L)) {
+    if(length(phi) == 1L) {
+      phi <- rep(phi, L)
+    }
+    if(length(n) == 1L) {
+      n <- rep(n, L)
+    }
+    if(length(m) == 1L) {
+      m <- rep(m, L)
+    }
   }
   ellPIcpp(phi, n, m, minerror)
   # if(phi == 0 || (is.infinite(Re(m)) && Im(m) == 0) ||
diff --git a/R/elliptic_Z.R b/R/elliptic_Z.R
index 2897b4f..4243989 100644
--- a/R/elliptic_Z.R
+++ b/R/elliptic_Z.R
@@ -1,12 +1,12 @@
 #' Jacobi zeta function
 #' @description Evaluate the Jacobi zeta function.
 #'
-#' @param phi amplitude, real or complex number
-#' @param m parameter, real or complex number
+#' @param phi amplitude, real or complex number/vector
+#' @param m parameter, real or complex number/vector
 #' @param minerror bound on relative error passed to \code{\link{elliptic_E}}
 #'   and \code{\link{elliptic_F}}
 #'
-#' @return A complex number, the value of the Jacobi zeta function
+#' @return A complex number or vector, the value(s) of the Jacobi zeta function
 #'   \ifelse{html}{\out{Z(&phi;,m)}}{\eqn{Z(\phi,m)}{Z(phi,m)}}.
 #' @export
 elliptic_Z <- function(phi, m, minerror = 1e-15) {
diff --git a/man/Lambda0.Rd b/man/Lambda0.Rd
index 0b8f97b..0970db8 100644
--- a/man/Lambda0.Rd
+++ b/man/Lambda0.Rd
@@ -7,15 +7,15 @@
 Lambda0(phi, m, minerror = 1e-14)
 }
 \arguments{
-\item{phi}{Jacobi amplitude, a complex number}
+\item{phi}{Jacobi amplitude, a complex number/vector}
 
-\item{m}{parameter, a complex number}
+\item{m}{parameter, a complex number/vector}
 
 \item{minerror}{the bound on the relative error passed to
 \code{\link{elliptic_F}} and \code{\link{elliptic_Z}}}
 }
 \value{
-A complex number.
+A complex number or vector.
 }
 \description{
 Evaluates the Heuman Lambda function.
diff --git a/man/elliptic_E.Rd b/man/elliptic_E.Rd
index 165a12b..8e4ece9 100644
--- a/man/elliptic_E.Rd
+++ b/man/elliptic_E.Rd
@@ -7,16 +7,16 @@
 elliptic_E(phi, m, minerror = 1e-15)
 }
 \arguments{
-\item{phi}{amplitude, real or complex number}
+\item{phi}{amplitude, real or complex number/vector}
 
-\item{m}{parameter, real or complex number}
+\item{m}{parameter, real or complex number/vector}
 
 \item{minerror}{the bound on the relative error passed to
 \code{\link{Carlson_RF}} and \code{\link{Carlson_RD}}}
 }
 \value{
-A complex number, the value of the incomplete elliptic integral
-  \ifelse{html}{\out{E(&phi;,m)}}{\eqn{E(\phi,m)}{E(phi,m)}}.
+A complex number or vector, the value(s) of the incomplete elliptic
+  integral \ifelse{html}{\out{E(&phi;,m)}}{\eqn{E(\phi,m)}{E(phi,m)}}.
 }
 \description{
 Evaluate the incomplete elliptic integral of the second kind.
diff --git a/man/elliptic_F.Rd b/man/elliptic_F.Rd
index 7937f55..3469a9f 100644
--- a/man/elliptic_F.Rd
+++ b/man/elliptic_F.Rd
@@ -7,16 +7,16 @@
 elliptic_F(phi, m, minerror = 1e-15)
 }
 \arguments{
-\item{phi}{amplitude, real or complex number}
+\item{phi}{amplitude, real or complex number/vector}
 
-\item{m}{parameter, real or complex number}
+\item{m}{parameter, real or complex number/vectot}
 
 \item{minerror}{the bound on the relative error passed to
 \code{\link{Carlson_RF}}}
 }
 \value{
-A complex number, the value of the incomplete elliptic integral
-  \ifelse{html}{\out{F(&phi;,m)}}{\eqn{F(\phi,m)}{F(phi,m)}}.
+A complex number or vector, the value(s) of the incomplete elliptic
+  integral \ifelse{html}{\out{F(&phi;,m)}}{\eqn{F(\phi,m)}{F(phi,m)}}.
 }
 \description{
 Evaluate the incomplete elliptic integral of the first kind.
diff --git a/man/elliptic_PI.Rd b/man/elliptic_PI.Rd
index 7c9fd42..d7ca268 100644
--- a/man/elliptic_PI.Rd
+++ b/man/elliptic_PI.Rd
@@ -7,17 +7,18 @@
 elliptic_PI(phi, n, m, minerror = 1e-15)
 }
 \arguments{
-\item{phi}{amplitude, real or complex number}
+\item{phi}{amplitude, real or complex number/vector}
 
-\item{n}{characteristic, real or complex number}
+\item{n}{characteristic, real or complex number/vector}
 
-\item{m}{parameter, real or complex number}
+\item{m}{parameter, real or complex number/vector}
 
 \item{minerror}{the bound on the relative error passed to
 \code{\link{Carlson_RF}} and \code{\link{Carlson_RJ}}}
 }
 \value{
-A complex number, the value of the incomplete elliptic integral
+A complex number or vector, the value(s) of the incomplete elliptic
+  integral
   \ifelse{html}{\out{&Pi;(&phi;,n,m)}}{\eqn{\Pi(\phi,n,m)}{PI(phi,n,m)}}.
 }
 \description{
diff --git a/man/elliptic_Z.Rd b/man/elliptic_Z.Rd
index c5095ea..14c2e35 100644
--- a/man/elliptic_Z.Rd
+++ b/man/elliptic_Z.Rd
@@ -7,15 +7,15 @@
 elliptic_Z(phi, m, minerror = 1e-15)
 }
 \arguments{
-\item{phi}{amplitude, real or complex number}
+\item{phi}{amplitude, real or complex number/vector}
 
-\item{m}{parameter, real or complex number}
+\item{m}{parameter, real or complex number/vector}
 
 \item{minerror}{bound on relative error passed to \code{\link{elliptic_E}}
 and \code{\link{elliptic_F}}}
 }
 \value{
-A complex number, the value of the Jacobi zeta function
+A complex number or vector, the value(s) of the Jacobi zeta function
   \ifelse{html}{\out{Z(&phi;,m)}}{\eqn{Z(\phi,m)}{Z(phi,m)}}.
 }
 \description{