docs: public float[] getQuantiles(float[] probabilities)
    {
      for (float p : probabilities) {
        Preconditions.checkArgument(0 < p && p < 1, "quantile probabilities must be strictly between 0 and 1");
      }
  
      float[] quantiles = new float[probabilities.length];
      Arrays.fill(quantiles, Float.NaN);
  
      if (this.count() == 0) {
        return quantiles;
      }
  
      final long[] bins = this.bins();
  
      for (int j = 0; j < probabilities.length; ++j) {
        // Adapted from Ben-Haiv/Tom-Tov (Algorithm 4: Uniform Procedure)
        // Our histogram has a set of bins {(p1, m1), (p2, m2), ... (pB, mB)}
        // and properties of min and max saved during construction, such
        // that min <= p1 and pB <= max.
        //
        // When min < p1 or pB < max, these are special cases of using the
        // dummy bins (p0, 0) or (pB+1, 0) respectively, where p0 == min
        // and pB+1 == max.
        //
        // This histogram gives an ordered set of numbers {pi, pi+1, ..., pn},
        // such that min <= pi < pn <= max, and n is the sum({m1, ..., mB}).
        // We use s to determine which pairs of (pi, mi), (pi+1, mi+1) to
        // calculate our quantile from by computing sum([pi,mi]) < s < sum
        // ([pi+1, mi+1])
        final double s = probabilities[j] * this.count();
  
        int i = 0;
        int sum = 0;
        int k = 1;
        long count;
        while (k <= this.binCount()) {
          count = bins[k - 1];
          if (sum + count > s) {
            i = k - 1;
            break;
          } else {
            sum += count;
          }
          ++k;
        }
  
        if (i == 0) {
          // At the first bin, there are no points to the left of p (min)
          quantiles[j] = this.min();
        } else {
          final double d = s - sum;
          final double c = -2 * d;
          final long a = bins[i] - bins[i - 1];
          final long b = 2 * bins[i - 1];
          double z;
          if (a == 0) {
            z = -c / b;
          } else {
            z = (-b + Math.sqrt(b * b - 4 * a * c)) / (2 * a);
          }
          final double uj = this.positions[i - 1] + (this.positions[i] - this.positions[i - 1]) * z;
          // A given bin (p, m) has m/2 points to the left and right of p, and
          // uj could be one of those points. However, uj is still subject to:
          // [min] (p0, 0) < uj < (p, m) or
          //                      (p, m) < uj < (pB+1, 0) [max]
          quantiles[j] = ((float) uj < this.max()) ? (float) uj : this.max();
        }
      }
  
      return quantiles;
    }