///|
pub struct Distribution {
  priv probabilities : Array[Double]
  priv tail : Double
  priv terms : Int32
  priv segments : Int32
}

///|
pub fn Distribution::masses(self : Distribution) -> Array[Double] {
  self.probabilities.copy()
}

///|
/// Conservative Poisson truncation bound, not a floating-point rounding bound.
pub fn Distribution::tail_bound(self : Distribution) -> Double {
  self.tail
}

///|
pub fn Distribution::work_terms(self : Distribution) -> Int32 {
  self.terms
}

///|
pub fn Distribution::segment_count(self : Distribution) -> Int32 {
  self.segments
}

///|
fn uniform_step(
  g : Generator,
  p : Array[Double],
  lambda : Double,
) -> Array[Double] raise KineticsError {
  let result = Array::make(g.size(), 0.0)
  for i in 0.. Distribution raise KineticsError {
  validate_probability(initial, self.size())
  nonnegative(time, "CTMC time")
  finite(tolerance, "Poisson tolerance")
  if tolerance < 1.0e-15 || tolerance > 1.0e-3 {
    raise Invalid("tolerance outside 1e-15..1e-3")
  }
  if max_terms < 1 ||
    max_terms > 1000000 ||
    max_work < 1 ||
    max_work > 100000000 {
    raise Invalid("uniformization work budget")
  }
  let mut lambda = 0.0
  let mut per_term = self.size()
  for i in 0.. lambda {
      lambda = self.exits[i]
    }
    per_term += self.edges[i].length()
  }
  if lambda == 0.0 || time == 0.0 {
    return { probabilities: initial.copy(), tail: 0.0, terms: 0, segments: 0, }
  }
  let mean = lambda * time
  finite(mean, "CTMC rate-time product")
  if mean > 160000.0 {
    raise Budget("uniformization at most 10000 time segments")
  }
  let mut segments = (mean / 16.0).to_int()
  if segments.to_double() * 16.0 < mean {
    segments += 1
  }
  if segments == 0 {
    segments = 1
  }
  let mu = mean / segments.to_double()
  let epsilon = tolerance / segments.to_double()
  let mut p = initial.copy()
  let mut tail = 0.0
  let mut terms = 0
  let mut work = 0
  for _ in 0.. x * weight)
    let mut k = 0
    while true {
      if terms >= max_terms || per_term > max_work - work {
        raise Budget("uniformization term/work budget exhausted")
      }
      terms += 1
      work += per_term
      let next_weight = weight * mu / (k + 1).to_double()
      let ratio = mu / (k + 2).to_double()
      if ratio < 1.0 {
        let bound = mass * next_weight / (1.0 - ratio)
        if bound <= epsilon {
          tail += bound
          break
        }
      }
      term = uniform_step(self, term, lambda)
      weight = next_weight
      k += 1
      for i in 0..