diff --git a/R/dist_in_torus.R b/R/dist_in_torus.R
new file mode 100644
index 0000000..db62966
--- /dev/null
+++ b/R/dist_in_torus.R
@@ -0,0 +1,82 @@
+#' A general function to calculate distances in a n-dimensional toroid.
+#'
+#' @description
+#' By default (i.e. if the arguments `lower` and `upper` are not provided), this
+#' function returns the distance in the Euclidean space by assuming borders
+#' infinitely apart (i.e. points in a small portion of an infinitely large
+#' toroid).
+#'
+#' The shortest distance in the toroid is the hypotenuse of the smallest
+#' hyper-triangle. The 'internal' distance is the typical distance based on the
+#' coordinates, as in the Euclidean space. The 'external' distance is crossing
+#' borders, going around. There are only two ways of measuring distance along
+#' each dimension.
+#'
+#' @param x A numeric matrix giving the coordinates (positions) of the points.
+#' @param lower,upper Numeric vectors of length `ncol(x)`. The minimum and
+#'   maximum possible values of the coordinates along each dimension.
+#'
+#' @author Gabriel Arellano
+#'
+#' @return A numeric matrix.
+#' @export
+#'
+#' @examples
+#' numeric_vec <- c(runif(10, min = 3, max = 5), runif(10, min = 13, max = 15))
+#' x <- matrix(numeric_vec, ncol = 2)
+#'
+#' # Euclidean distances
+#' d0 <- dist(x)
+#' # default behaviour
+#' d1 <- dist_in_torus(x)
+#' # distances in the toroid
+#' d2 <- dist_in_torus(x, lower = c(3, 13), upper = c(5, 15))
+#'
+#' par(mfrow = c(1, 3))
+#' plot(x, xlim = c(3, 5), ylim = c(13, 15), xlab = "x", ylab = "y")
+#' plot(c(d0), c(as.dist(d1)), main = "default = Euclidean = infinite toroid")
+#' abline(0, 1)
+#' plot(c(d0), c(as.dist(d2)), main = "finite toroid")
+#' abline(0, 1)
+#'
+#' # `upper` and `lower` must be as long as `ncol(x)`
+#' x <- matrix(runif(9), ncol = 3)
+#' dist_in_torus(x, lower = c(0, 0, 0), upper = c(1, 1, 1))
+dist_in_torus <- function(x,
+                          lower = rep(-Inf, ncol(x)),
+                          upper = rep(Inf, ncol(x))) {
+  if (!is.matrix(x)) {
+    warn(paste0("Coercing `x` to matrix.\n* `x` was of class ", class(x)))
+    x <- as.matrix(x)
+  }
+  check_dist_in_torus(x = x, lower = lower, upper = upper)
+
+  # Number of dimensions
+  n <- ncol(x)
+  # Size of the n-dimensional space considered
+  ranges <- upper - lower
+
+  # Internal and external cathetuses along each dimension:
+  internal_cats <- sapply(
+    1:n, function(i) abs(outer(x[, i], x[, i], "-")),
+    simplify = "array"
+  )
+  external_cats <- sapply(1:n, function(i) ranges[i] - internal_cats[, , i])
+
+  # The shortest cathetuses along each dimension define the smallest
+  # hyper-triangle:
+  shortest_cats <- pmin(internal_cats, external_cats)
+
+  # Application of the Pythagorean theorem across layers:
+  hypo <- sqrt(rowSums(shortest_cats^2, dims = 2))
+  hypo
+}
+
+check_dist_in_torus <- function(x, lower, upper) {
+  if (!is.numeric(x)) {
+    msg <- paste0("`x` must be numeric.\n", "* It has type ", typeof(x))
+    abort(msg)
+  }
+
+  stopifnot(length(lower) == ncol(x), length(upper) == ncol(x))
+}
diff --git a/tests/testthat/test-dist_in_torus.R b/tests/testthat/test-dist_in_torus.R
new file mode 100644
index 0000000..342bde3
--- /dev/null
+++ b/tests/testthat/test-dist_in_torus.R
@@ -0,0 +1,51 @@
+context("dist_in_torus")
+
+test_that("outputs is a matrix of doubles", {
+  x <- matrix(0:3, nrow = 2)
+  d1 <- dist_in_torus(x)
+  expect_type(d1, "double")
+  expect_true("matrix" %in% class(d1))
+
+  d2 <- dist_in_torus(x, lower = c(3, 13), upper = c(5, 15))
+  expect_type(d2, "double")
+  expect_true("matrix" %in% class(d2))
+})
+
+test_that("fails with wrong input", {
+  chr <- matrix(letters[1:4], nrow = 2)
+  expect_error(
+    dist_in_torus(chr), "`x` must be numeric"
+  )
+
+  x <- matrix(1:9, nrow = 3)
+  expect_error(dist_in_torus(x, lower = c(0, 0), upper = c(1, 1, 1)))
+  expect_error(dist_in_torus(x, lower = c(0, 0, 0), upper = c(1, 1)))
+})
+
+test_that("warns if input is not a matrix", {
+  x <- matrix(1:4, nrow = 2)
+  expect_warning(
+    dist_in_torus(as.data.frame(x)), "Coercing `x` to matrix."
+  )
+
+  x <- 1:4
+  expect_warning(
+    dist_in_torus(x)
+  )
+})
+
+test_that("behaves in particular ways with extreeme conditions", {
+  # This is mainly to document the behaviour. Not to say it is OK. I don't know.
+  # Passes
+  expect_silent(dist_in_torus(matrix(c(NaN, NaN))))
+  expect_silent(dist_in_torus(matrix(c(Inf, Inf))))
+  expect_silent(dist_in_torus(matrix(c(1, 1))))
+  expect_silent(dist_in_torus(matrix(c(-1, -1))))
+  expect_silent(dist_in_torus(matrix(c(1, 1))))
+  # Warns
+  expect_warning(dist_in_torus(data.frame(a = c(1, 1))))
+  # Fails
+  expect_error(dist_in_torus(matrix(c(NA, NA))))
+  expect_error(dist_in_torus(matrix(1)))
+  expect_error(dist_in_torus(matrix(c(NULL, NULL))))
+})