diff --git a/R/beeswarm.R b/R/beeswarm.R
index 39ff24f..1cca008 100644
--- a/R/beeswarm.R
+++ b/R/beeswarm.R
@@ -20,6 +20,7 @@ beeswarm.default <- function(x,
     corralWidth, side = 0L, 
     priority = c("ascending", "descending", "density", "random", "none"),
     fast = TRUE,
+    epsilon = 0,
     pch = par("pch"), col = par("col"), bg = NA, 
     pwpch = NULL, pwcol = NULL, pwbg = NULL, pwcex = NULL,
     do.plot = TRUE, add = FALSE, axes = TRUE, log = FALSE, 
@@ -196,11 +197,11 @@ beeswarm.default <- function(x,
     if(horizontal) {
       g.offset <- lapply(x, function(a) swarmy(x = a, y = rep(0, length(a)), 
           cex = sizeMultiplier, side = side, priority = priority,
-          fast = fast, compact = compact)$y)
+          fast = fast, compact = compact, epsilon = epsilon)$y)
     } else {
       g.offset <- lapply(x, function(a) swarmx(x = rep(0, length(a)), y = a, 
           cex = sizeMultiplier, side = side, priority = priority,
-          fast = fast, compact = compact)$x)
+          fast = fast, compact = compact, epsilon = epsilon)$x)
     }
     d.pos <- x
   } else {          ####   non-swarm methods
@@ -402,7 +403,7 @@ beeswarm.formula <- function (formula, data = NULL, subset, na.action = NULL,
 }
 
 
-.calculateCompactSwarm <- function(x, dsize, gsize, side = 0L, priority = "ascending") {
+.calculateCompactSwarm <- function(x, dsize, gsize, side = 0L, priority = "ascending", epsilon = 0) {
   if(length(x) == 0) return(numeric(0))
   stopifnot(side %in% -1:1)
 
@@ -428,8 +429,9 @@ beeswarm.formula <- function (formula, data = NULL, subset, na.action = NULL,
   #### place the points
   if(nrow(out) > 1) {
     for (iter in 1:nrow(out)) {          ## we will place one point at a time
-      i <- which.min(abs(out$y.best))    ## Choose a point that can be placed
-                                         ## close to non-data axis
+      min_abs = min(abs(out$y.best))
+      ## Choose a point that can be placed close to non-data axis
+      i <- min(which(abs(out$y.best) <= min_abs + epsilon))
       xi <- out$x[i]
       yi <- out$y[i] <- out$y.best[i]
       out$placed[i] <- TRUE
@@ -461,7 +463,8 @@ beeswarm.formula <- function (formula, data = NULL, subset, na.action = NULL,
 
 
 .calculateSwarmUsingC <- function(x, dsize, gsize, side = 0L, priority = "ascending",
-    compact = FALSE) {
+    compact = FALSE, epsilon = 0) {
+  print(x)
   if(length(x) == 0) return(numeric(0))
   stopifnot(side %in% -1:1)
 
@@ -489,11 +492,12 @@ beeswarm.formula <- function (formula, data = NULL, subset, na.action = NULL,
                n = n,
                compact = as.logical(compact),
                side = as.integer(side),
+               epsilon = epsilon,
                placed = integer(n),         # used internally by C implementations
                workspace = numeric(n * 4),  # used internally by C implementations
                y = numeric(n))
   y <- rep(NA, length(x))
-  y[out$index] <- result[[7]] * gsize
+  y[out$index] <- result[[8]] * gsize
   y
 }
 
@@ -505,7 +509,10 @@ swarmx <- function(x, y,
     log = NULL, cex = par("cex"), side = 0L, 
     priority = c("ascending", "descending", "density", "random", "none"),
     fast = TRUE,
-    compact = FALSE) { 
+    compact = FALSE,
+    epsilon = 0) { 
+  stopifnot(is.numeric(epsilon))
+  stopifnot(!is.na(epsilon))
   priority <- match.arg(priority)
   if(is.null(log)) 
     log <- paste(ifelse(par('xlog'), 'x', ''), ifelse(par('ylog'), 'y', ''), sep = '')
@@ -517,11 +524,16 @@ swarmx <- function(x, y,
   if(ylog) xy$y <- log10(xy$y)
   if (fast) {
     x.new <- xy$x + .calculateSwarmUsingC(xy$y, dsize = ysize * cex,
-      gsize = xsize * cex, side = side, priority = priority, compact = compact)
+      gsize = xsize * cex, side = side, priority = priority, compact = compact,
+      epsilon = epsilon)
   } else {
-    swarmFn <- ifelse(compact, .calculateCompactSwarm, .calculateSwarm)
-    x.new <- xy$x + swarmFn(xy$y, dsize = ysize * cex, gsize = xsize * cex,
-      side = side, priority = priority)
+    if(compact){
+      x.new <- xy$x + .calculateCompactSwarm(xy$y, dsize = ysize * cex, gsize = xsize * cex,
+        side = side, priority = priority, epsilon = epsilon)
+    } else {
+      x.new <- xy$x + .calculateSwarm(xy$y, dsize = ysize * cex, gsize = xsize * cex,
+        side = side, priority = priority)
+    }
   }
   out <- data.frame(x = x.new, y = y)
   if(xlog) out$x <- 10 ^ out$x
@@ -535,7 +547,10 @@ swarmy <- function(x, y,
     log = NULL, cex = par("cex"), side = 0L, 
     priority = c("ascending", "descending", "density", "random", "none"),
     fast = TRUE,
-    compact = FALSE) { 
+    compact = FALSE,
+    epsilon = 0) { 
+  stopifnot(is.numeric(epsilon))
+  stopifnot(!is.na(epsilon))
   priority <- match.arg(priority)
   if(is.null(log)) 
     log <- paste(ifelse(par('xlog'), 'x', ''), ifelse(par('ylog'), 'y', ''), sep = '')
@@ -547,11 +562,16 @@ swarmy <- function(x, y,
   if(ylog) xy$y <- log10(xy$y)
   if (fast) {
     y.new <- xy$y + .calculateSwarmUsingC(xy$x, dsize = xsize * cex,
-      gsize = ysize * cex, side = side, priority = priority, compact = compact)
+      gsize = ysize * cex, side = side, priority = priority, compact = compact,
+      epsilon = epsilon)
   } else {
-    swarmFn <- ifelse(compact, .calculateCompactSwarm, .calculateSwarm)
-    y.new <- xy$y + swarmFn(xy$x, dsize = xsize * cex, gsize = ysize * cex,
-      side = side, priority = priority)
+    if(compact){
+      y.new <- xy$y + .calculateCompactSwarm(xy$x, dsize = xsize * cex, gsize = ysize * cex,
+        side = side, priority = priority, epsilon = epsilon)
+    } else {
+      y.new <- xy$y + .calculateSwarm(xy$x, dsize = xsize * cex, gsize = ysize * cex,
+        side = side, priority = priority)
+    }
   }
   out <- data.frame(x = x, y = y.new)
   if(ylog) out$y <- 10 ^ out$y
diff --git a/man/beeswarm.Rd b/man/beeswarm.Rd
index c7ceecc..dd9e8af 100644
--- a/man/beeswarm.Rd
+++ b/man/beeswarm.Rd
@@ -24,6 +24,7 @@ beeswarm(x, \dots)
     corralWidth, side = 0L, 
     priority = c("ascending", "descending", "density", "random", "none"),
     fast = TRUE,
+    epsilon = 0,
     pch = par("pch"), col = par("col"), bg = NA, 
     pwpch = NULL, pwcol = NULL, pwbg = NULL, pwcex = NULL,
     do.plot = TRUE, add = FALSE, axes = TRUE, log = FALSE,
@@ -56,6 +57,7 @@ beeswarm(x, \dots)
   \item{side}{ Direction to perform jittering: 0: both directions; 1: to the right or upwards; -1: to the left or downwards.}
   \item{priority}{ Order used to perform point layout when method is \code{"swarm"} or \code{"compactswarm"}; ignored otherwise (see Details).}
   \item{fast}{ Use compiled version of algorithm?  This option is ignored for all methods except \code{"swarm"} and \code{"compactswarm"}.}
+  \item{epsilon}{Tolerance for floating-point comparisons.  The default value of zero usually works well, but in some cases a tiny value of \code{epsilon} such as 0.0000001 gives a more compact arrangement of points.  This option is ignored for all methods except \code{"compactswarm"}.}
   \item{pch, col, bg}{ Plotting characters and colors, specified by group. Recycled if necessary (see Details). }
   \item{pwpch, pwcol, pwbg, pwcex}{ \dQuote{Point-wise} plotting characteristics, specified for each data point (see Details). }
   \item{do.plot}{ Draw a plot? }
diff --git a/man/swarmx.Rd b/man/swarmx.Rd
index 0583da2..ffc4f6f 100644
--- a/man/swarmx.Rd
+++ b/man/swarmx.Rd
@@ -9,13 +9,13 @@ swarmx(x, y,
     ysize = yinch(0.08, warn.log = FALSE),
     log = NULL, cex = par("cex"), side = 0L, 
     priority = c("ascending", "descending", "density", "random", "none"),
-    fast = TRUE, compact = FALSE)
+    fast = TRUE, compact = FALSE, epsilon = 0)
 swarmy(x, y, 
     xsize = xinch(0.08, warn.log = FALSE), 
     ysize = yinch(0.08, warn.log = FALSE),
     log = NULL, cex = par("cex"), side = 0L, 
     priority = c("ascending", "descending", "density", "random", "none"),
-    fast = TRUE, compact = FALSE)
+    fast = TRUE, compact = FALSE, epsilon = 0)
 }
 \arguments{
   \item{x, y}{ Coordinate vectors in any format supported by \code{\link{xy.coords}}. }
@@ -26,6 +26,7 @@ swarmy(x, y,
   \item{priority}{ Method used to perform point layout (see below).}
   \item{fast}{ Use compiled version of algorithm?  This option is ignored for all methods except \code{"swarm"} and \code{"compactswarm"}.}
   \item{compact}{ Use compact layout? (see below)}
+  \item{epsilon}{Tolerance for floating-point comparisons.  The default value of zero usually works well, but in some cases a tiny value of \code{epsilon} such as 0.0000001 gives a more compact arrangement of points.  This option is ignored for all methods except \code{"compactswarm"}.}
 }
 \details{ 
 For \code{swarmx}, the input coordinates must lie in a vertical line.  For \code{swarmy}, the input coordinates must lie in a horizontal line.
diff --git a/src/beeswarm.c b/src/beeswarm.c
index d2ac95b..fa89e81 100644
--- a/src/beeswarm.c
+++ b/src/beeswarm.c
@@ -90,7 +90,7 @@ static int which_min_abs(double *y_best, int *placed, int n)
 }
 
 static void compactSwarm(double *x, int n, int side,
-    int *placed, double *workspace, double *y)
+    double epsilon, int *placed, double *workspace, double *y)
 {
   double *y_low = workspace;       // largest permitted negative y value
   double *y_high = workspace + n;  // smallest permitted positive y value
@@ -98,7 +98,16 @@ static void compactSwarm(double *x, int n, int side,
 
   for (int iter=0; iter<n; iter++) {
     R_CheckUserInterrupt();
+
     int i = which_min_abs(y_best, placed, n);
+    if (epsilon > 0) {
+      double min_abs = fabs(y_best[i]);
+      i = 0;
+      while (placed[i] || fabs(y_best[i]) > min_abs + epsilon) ++i;
+      // i is now the index of the first unplaced point whose y_best value is no
+      // greater than min_abs
+    }
+
     double xi = x[i];
     double yi = y_best[i];
     y[i] = yi;
@@ -132,15 +141,18 @@ static void compactSwarm(double *x, int n, int side,
  * n         length of x
  * compact   use compact layout?
  * side      -1, 0, or 1
+ * epsilon   for the compact beeswarm only: allow the placement of circles
+             whose distance from the non-data axis is up to epsilon more than
+             that of the closest point to the axis
  * placed    which circles have been placed (logical type)
  * workspace an array of doubles for internal use
  * y         (output) circle positions on non-data axis
  */
 void attribute_hidden calculateSwarm(double *x, int *n, int *compact, int *side,
-    int *placed, double *workspace, double *y)
+    double *epsilon, int *placed, double *workspace, double *y)
 {
   if (*compact) {
-    compactSwarm(x, *n, *side, placed, workspace, y);
+    compactSwarm(x, *n, *side, *epsilon, placed, workspace, y);
   } else {
     swarm(x, *n, *side, placed, workspace, y);
   }
diff --git a/src/beeswarm.h b/src/beeswarm.h
index b4536fd..599f47e 100644
--- a/src/beeswarm.h
+++ b/src/beeswarm.h
@@ -2,6 +2,6 @@
 #define BEESWARM_BEESWARM_H
 
 void calculateSwarm(double *x, int *n, int *compact, int *side,
-    int *placed, double *workspace, double *y);
+    double *epsilon, int *placed, double *workspace, double *y);
 
 #endif /* BEESWARM_BEESWARM_H */
diff --git a/src/init.c b/src/init.c
index 0d0004b..fd1e2e2 100644
--- a/src/init.c
+++ b/src/init.c
@@ -5,10 +5,10 @@
 #include "beeswarm.h"
 
 static R_NativePrimitiveArgType calculateSwarm_types[] = {
-    REALSXP, INTSXP, LGLSXP, INTSXP, INTSXP, REALSXP, REALSXP};
+    REALSXP, INTSXP, LGLSXP, INTSXP, REALSXP, INTSXP, REALSXP, REALSXP};
 
 static const R_CMethodDef CEntries[] = {
-    {"calculateSwarm", (DL_FUNC) &calculateSwarm, 7, calculateSwarm_types},
+    {"calculateSwarm", (DL_FUNC) &calculateSwarm, 8, calculateSwarm_types},
     {NULL, NULL, 0, NULL}
 };
 
diff --git a/tests/beeswarm-test.R b/tests/beeswarm-test.R
index 09a3d2f..ccb3158 100644
--- a/tests/beeswarm-test.R
+++ b/tests/beeswarm-test.R
@@ -6,24 +6,30 @@ test_swarms <- function(x) {
     for (compact in c(FALSE, TRUE)) {
       for (priority in c("ascending", "descending", "density", "random", "none")) {
         for (side in -1:1) {
-          print(na_positions)
-          print(compact)
-          print(priority)
-          print(side)
-          # compare R and C versions of swarmy and swarmx
-          set.seed(1)
-          y1 <- swarmy(x, numeric(length(x)), compact=compact, side=side, priority=priority, fast=TRUE)$y
-          set.seed(1)
-          y2 <- swarmy(x, numeric(length(x)), compact=compact, side=side, priority=priority, fast=FALSE)$y
-          stopifnot(all.equal(y1, y2))
-          stopifnot(identical(which(is.na(y1)), na_positions))
+          for (epsilon in c(0, 0.00001, 0.1)) {
+            if (epsilon > 0 & !compact) {
+              next;
+            }
+            print(na_positions)
+            print(compact)
+            print(priority)
+            print(side)
+            print(epsilon)
+            # compare R and C versions of swarmy and swarmx
+            set.seed(1)
+            y1 <- swarmy(x, numeric(length(x)), compact=compact, side=side, priority=priority, fast=TRUE, epsilon=epsilon)$y
+            set.seed(1)
+            y2 <- swarmy(x, numeric(length(x)), compact=compact, side=side, priority=priority, fast=FALSE, epsilon=epsilon)$y
+            stopifnot(all.equal(y1, y2))
+            stopifnot(identical(which(is.na(y1)), na_positions))
 
-          set.seed(1)
-          x1 <- swarmx(numeric(length(x)), x, compact=compact, side=side, priority=priority, fast=TRUE)$x
-          set.seed(1)
-          x2 <- swarmx(numeric(length(x)), x, compact=compact, side=side, priority=priority, fast=FALSE)$x
-          stopifnot(all.equal(x1, x2))
-          stopifnot(identical(which(is.na(x1)), na_positions))
+            set.seed(1)
+            x1 <- swarmx(numeric(length(x)), x, compact=compact, side=side, priority=priority, fast=TRUE, epsilon=epsilon)$x
+            set.seed(1)
+            x2 <- swarmx(numeric(length(x)), x, compact=compact, side=side, priority=priority, fast=FALSE, epsilon=epsilon)$x
+            stopifnot(all.equal(x1, x2))
+            stopifnot(identical(which(is.na(x1)), na_positions))
+          }
         }
       }
     }