///|
/// Erlang loss probability by bounded recursion; no factorials or a^c.
/// Offered traffic is dimensionless lambda/mu, not an arrival rate.
pub fn erlang_b(offered : Double, channels : Int) -> Double raise QueueError {
  nonnegative(offered, "offered traffic")
  if offered > 1.0e12 || channels < 0 || channels > 10000 {
    raise Limit("Erlang B: traffic <=1e12, channels 0..10000")
  }
  let mut loss = 1.0
  for k in 1..<=channels {
    let numerator = offered * loss
    loss = numerator / (k.to_double() + numerator)
  }
  loss
}

///|
/// Delay probability in an infinite-buffer M/M/c model. Stability is strict.
pub fn erlang_c(offered : Double, channels : Int) -> Double raise QueueError {
  servers(channels)
  nonnegative(offered, "offered traffic")
  let rho = offered / channels.to_double()
  if rho >= 1.0 {
    raise Unstable("M/M/c requires offered traffic < server count")
  }
  let b = erlang_b(offered, channels)
  let c = b / (1.0 - rho + rho * b)
  finite(c, "Erlang C")
  c
}

///|
/// Minimum channels for a loss target. Zero channels may serve zero traffic
/// but blocks every attempted arrival; use at least one for actual service.
pub fn loss_staffing(
  offered : Double,
  maximum_loss : Double,
  maximum_channels : Int,
) -> Int raise QueueError {
  nonnegative(offered, "offered traffic")
  probability(maximum_loss, "maximum loss")
  if offered > 1.0e12 || maximum_channels < 0 || maximum_channels > 10000 {
    raise Limit("loss staffing bounds")
  }
  let mut loss = 1.0
  if loss <= maximum_loss {
    return 0
  }
  for c in 1..<=maximum_channels {
    let x = offered * loss
    loss = x / (c.to_double() + x)
    if loss <= maximum_loss {
      return c
    }
  }
  raise Infeasible("loss target unattainable within channel bound")
}

///|
/// Invert Erlang B in offered traffic on a caller-supplied bounded interval.
/// Returns the greatest feasible traffic to relative interval tolerance.
pub fn loss_traffic_limit(
  channels : Int,
  maximum_loss : Double,
  upper : Double,
) -> Double raise QueueError {
  servers(channels)
  probability(maximum_loss, "maximum loss")
  nonnegative(upper, "traffic upper bound")
  if upper > 1.0e12 {
    raise Limit("traffic upper bound")
  }
  if erlang_b(upper, channels) <= maximum_loss {
    return upper
  }
  let mut lo = 0.0
  let mut hi = upper
  for _ in 0..<100 {
    let mid = lo + (hi - lo) * 0.5
    if erlang_b(mid, channels) <= maximum_loss {
      lo = mid
    } else {
      hi = mid
    }
  }
  lo
}