minQuadratic <- function (y, f, g, k, lw, up) {
  func <- function (x, f, g, k) {
    return (f + g * (x - y) + 0.5 * k * (x - y) ^ 2)
  }
  fup <- func (up, f, g, k)
  flw <- func (lw, f, g, k)
  if (k <= 0) {
    if (fup <= flw)
      return (list (x = up, f = fup))
    if (fup >= flw)
      return (list (x = lw, f = flw))
  }
  xmn <- y - g / k
  fmn <- func (xmn, f, g, k)
  if (xmn >= up)
    return (list (x = up, f = fup))
  if (xmn <= lw)
    return (list (x = lw, f = flw))
  return (list (x =  xmn, f = fmn))
}

minimaxQ <- function (y, f, g, k, lw, up) {
  n <- length (f)
  res <- Inf
  for (i in 1:n) {
    xlw <- lw
    xup <- up
    fail <- FALSE
    for (j in 1:n) {
      if (j == i)
        next
      df <- f[i] - f[j]
      dg <- g[i] - g[j]
      if (dg > 0)
        xlw <- max (xlw, y - df / dg)
      if (dg < 0)
        xup <- min (xup, y - df / dg)
      if ((dg == 0) && (f[i] < f[j]))
        fail <- TRUE
    }
    if (xlw >= xup)
      fail <- TRUE
    if (!fail)
      h <- minQuadratic (y, f[i], g[i], k, xlw, xup)
    if (h$f < res) {
      res <- h$f
      sol <- h$x
    }
  }
  return (sol)
}

minimaxB <- function (y, f, g, k, lw, up) {
  n <- length (f)
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      df <- f[i] - f[j]
      dg <- g[i] - g[j]
      dk <- k[i] - k[j]
      rt <- y + solve (polynomial (c (df, dg, dk / 2)))
      rt <-
        rt[which ((rt <= up) & (rt >= lw) & (!is.complex (rt)))]
    }
  }
  rt <- sort (unique (c(lw, up, rt)))
  m <- length (rt)
  fmin <- Inf
  for (j in 1:(m - 1)) {
    x <- (rt[j] + rt[j + 1]) / 2
    smax <- -Inf
    for (i in 1:n) {
      si <- f[i] + g[i] * (x - y) + 0.5 * k[i] * (x - y) ^ 2
      if (si > smax) {
        smax <- si
        l <- i
      }
    }
    h <- minQuadratic (y, f[l], g[l], k[l], rt[j], rt[j + 1])
    if (h$f < fmin) {
      fmin <- h$f
      xmin <- h$x
    }
  }
  return (xmin)
}