mm <- function(x) {
  x    <- sort(x)       
  n    <- length(x)
  xbar <- mean(x)
  xvar <- var(x)
  r    <- (2 * x - x[n] - x[1]) / (x[n] - x[1])

  rlo  <- r[1:(n - 1)]     
  rhi  <- r[2:n]
  rmid <- r[2:(n - 1)]

  xlo  <- x[1:(n - 1)]
  xhi  <- x[2:n]
  xmid <- x[2:(n - 1)]

  aa   <- 1 / (3 * (n - 1)) * (sum(rlo * rlo) + sum(rlo * rhi) + sum(rhi * rhi)) - 1 / (n - 1) ^ 2 * (sum(rmid) ^ 2)
  bb   <- 1 / (3 * (n - 1)) * (sum(2 * xlo * rlo) + sum(xlo * rhi + xhi * rlo) + sum(2 * xhi * rhi)) - 1 / (n - 1) ^ 2 * sum(rmid) * (x[1] + x[n] + 2 * sum(xmid))
  cc   <- 1 / (3 * (n - 1)) * (sum(xlo ^ 2) + sum(xlo * xhi) + sum(xhi ^ 2)) - 1 / (4 * (n - 1) ^ 2) * ((x[1] + x[n] + 2 * sum(xmid)) ^ 2) - xvar

  del  <- (-bb + sqrt(bb ^ 2 - 4 * aa * cc)) / (2 * aa)
  xp   <- x + r * del
  xpp  <- xp - ((sum(xp) - xp[1] / 2 - xp[n] / 2) / (n - 1) - xbar)
  xpp
}