public static class NumericalRoutines
{  
  /// <summary>
  /// Точность вычислений гамма-функции по алгоритму Lanczos’а
  /// </summary>
  public const double NUM_LANCZOS_E = 1.0e-15;

  // Размерность вектора, задающего точность аппроксимации Lanczos’а (g)
  private const int NUM_LANCZOS_G = 7;

  // Вектор множителей P для аппроксимации Lanczos’а
  private static readonly double [] NUM_LANCZOS_P = new double []  
  {
  9.9999999999980993e-1,  
  6.7652036812188510e+2,  
  -1.2591392167224028e+3,  
  7.7132342877765313e+2,  
  -1.7661502916214059e+2,  
  1.2507343278686905e+1,  
  -1.3857109526572012e-1,  
  9.9843695780195716e-6,  
  1.5056327351493116e-7
  };

  /// <summary>
  /// Натуральный логарифм гамма-функции от x > 0, рассчитывается через аппроксимацию Lanczos’а (с точностью до  
  /// 15-ого знака)
  /// </summary>
  public static double GammaLogarithmByLanczos (double x)
  {
  if (x > 0.5)
  {
  double arg = x - 1.0;
  double XP = NUM_LANCZOS_P [0];
  double shift = arg + NUM_LANCZOS_G + 0.5;

  for (int i = 1; i < NUM_LANCZOS_P.Length; i++)
  {
  XP += NUM_LANCZOS_P [i] / (arg + i);
  }

  return MathConst.MATH_LN_SQRT_2_PI + Math.Log (XP) + (arg + 0.5) * Math.Log (shift) - shift;
  }

  else if (x < 0.5)
  {
  return MathConst.MATH_LN_PI / ((NumericalRoutines.GammaLogarithmByLanczos (Math.Sin (Math.PI * x)) +
  NumericalRoutines.GammaLogarithmByLanczos (1.0 - x)));
  }

  else return MathConst.MATH_LN_SQRT_PI;
  }

  /// <summary>
  /// Неполная бета-функция
  /// </summary>
  /// <param name="x">Аргумент функции</param>
  /// <param name="a">Параметр функции</param>
  /// <param name="b">Параметр функции</param>
  /// <param name="epsilon">Точность вычислений</param>
  public static double IncompleteBeta (double x, double a, double b, double epsilon)
  {
  if ((x < 0.0) || (x > 1.0)) throw new ArgumentOutOfRangeException ();

  double coeffB;
  if ((x == 0.0) || (x == 1.0)) coeffB = 0.0;

  coeffB = Math.Exp (NumericalRoutines.GammaLogarithmByLanczos (a + b) -  
  NumericalRoutines.GammaLogarithmByLanczos (a) -
  NumericalRoutines.GammaLogarithmByLanczos (b) + a * Math.Log (x) + b * Math.Log (1.0 - x));

  if (x < (a + 1.0) / (a + b + 2.0))
  {
  return coeffB * NumericalRoutines.IncompleteBetaContinuedFraction (x, a, b, epsilon) / a;
  }

  else
  {
  return 1.0 - coeffB * NumericalRoutines.IncompleteBetaContinuedFraction (1.0 - x, b, a, epsilon) / b;
  }
  }

  /// <summary>
  /// Расчёт дроби неполной бета-функции с помощью разложения в цепные дроби
  /// </summary>
  /// <param name="x">Аргумент функции</param>
  /// <param name="a">Параметр функции</param>
  /// <param name="b">Параметр функции</param>
  /// <param name="epsilon">Точность вычислений</param>
  private static double IncompleteBetaContinuedFraction (double x, double a, double b, double epsilon)
  {
  double qab = a + b;
  double qap = a + 1.0;
  double qam = a - 1.0;
  double c = 1.0;
  double d = 1.0 / (1.0 - qab * x / qap);
   
  double aa, del;
  double result = d;

  int m = 1, m2 = 2;

  do
  {
  aa = m * (b - m) * x / ((qam + m2) * (a + m2));
  d = 1.0 / (1.0 + aa * d);
  c = 1.0 + aa / c;
  result *= d * c;
   
  aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
  d = 1.0 / (1.0 + aa * d);
  c = 1.0 + aa / c;
  del = d * c;
  result *= del;

  m++;
  m2 += 2;

  } while (Math.Abs (1.0 - del) >= epsilon);

  return result;
  }
}