docs: def ks_unif_pelz_good(samples, statistic):
      """
      Approximates the statistic distribution by a transformed Li-Chien formula.
  
      This ought to be a bit more accurate than using the Kolmogorov limit, but
      should only be used with large squared sample count times statistic.
      See: doi:10.18637/jss.v039.i11 and http://www.jstor.org/stable/2985019.
      """
      x = 1 / statistic
      r2 = 1 / samples
      rx = sqrt(r2) * x
      r2x = r2 * x
      r2x2 = r2x * x
      r4x = r2x * r2
      r4x2 = r2x2 * r2
      r4x3 = r2x2 * r2x
      r5x3 = r4x2 * rx
      r5x4 = r4x3 * rx
      r6x3 = r4x2 * r2x
      r7x5 = r5x4 * r2x
      r9x6 = r7x5 * r2x
      r11x8 = r9x6 * r2x2
      a1 = rx * (-r6x3 / 108 + r4x2 / 18 - r4x / 36 - r2x / 3 + r2 / 6 + 2)
      a2 = pi2 / 3 * r5x3 * (r4x3 / 8 - r2x2 * 5 / 12 - r2x * 4 / 45 + x + 1 / 6)
      a3 = pi4 / 9 * r7x5 * (-r4x3 / 6 + r2x2 / 4 + r2x * 53 / 90 - 1 / 2)
      a4 = pi6 / 108 * r11x8 * (r2x2 / 6 - 1)
      a5 = pi2 / 18 * r5x3 * (r2x / 2 - 1)
      a6 = -pi4 * r9x6 / 108
      w = -pi2 / 2 * r2x2
      return hpi1d2 * ((a1 + (a2 + (a3 + a4 * hs2) * hs2) * hs2) * exp(w * hs2) +
                       (a5 + a6 * is2) * is2 * exp(w * is2)).sum()