public static partial class DoubleExtensions
{
    public const double Pi     = 3.14159265358979323846d; // System.Math.PI;
    public const double SqrtPi = 1.77245385090551602729d; // System.Math.Sqrt(System.Math.PI);

    public static double Pow(this double value, int power) { double r = 1d; while (--power >= 0) r *= value; while (++power < 0) r /= value; return r; }
    public static double Pochhammer(this double value, int n) { double r = 1d; while (--n >= 0) r *= value + n; while (++n < 0) r /= value + n; return r; }
    public static double Pochhammer(this double value, double n) { return (value + n).Gamma() / value.Gamma(); }
    public static double Gamma(this double value)
    {
        if (double.IsNaN(value)) return double.NaN;
        if (double.IsPositiveInfinity(value)) return double.PositiveInfinity;
        if (double.IsNegativeInfinity(value)) return double.NaN;
        double q = 171.62437695630272d;
        if (value > q) return double.PositiveInfinity;
        q = value % 1d;
        if (q < 0d) ++q;
        var z = (int)(value - q);
        value = q + 8d;
        return q.Equals(0d) || q.Equals(1d) ? (q + 1d).Pochhammer(z - 1) :
            q.Equals(0.5d) ? 0.5d * SqrtPi * (q + 1d).Pochhammer(z - 1) :
            System.Math.Exp(
            0.5d * System.Math.Log(Pi + Pi)
            + (value - 0.5d) * System.Math.Log(value) - value
            + value.Pow(-1) / 12d    // There we using "http://oeis.org/A046969/list" like denominator
            - value.Pow(-3) / 360d    // and "http://oeis.org/A001067/list" like numerator.
            + value.Pow(-5) / 1260d    // More info: "http://mathworld.wolfram.com/StirlingsSeries.html"
            - value.Pow(-7) / 1680d
            + value.Pow(-9) / 1188d
            - value.Pow(-11) / 360360d * 691d
            + value.Pow(-13) / 156d
        //  - value.Pow(-15) / 122400 * 3617
        //  + value.Pow(-17) / 244188 * 43867
        //  ...
            ) * value.Pochhammer(z - 8);
    }
    public static double Combination(this double n, uint k) { return (n - k + 1).Pochhammer(k) / k.Factorial(); }
    public static double Combination(this double n, double k) { return (n - k + 1).Pochhammer(k) / (++k).Gamma(); }
}

public static partial class IntExtensions
{
    public static ulong Pow(this uint value, uint power)
    {
        ulong r = 1uL, m = value;
        while (true)
        {
            if ((power & 1) != 0) r *= m;
            if ((power >>= 1) == 0) return r;
            m *= m;
        }
    }
    public static long Pow(this int value, uint power)
    {
        long r = 1L, m = value;
        while (true)
        {
            if ((power & 1) != 0) r *= m;
            if ((power >>= 1) == 0) return r;
            m *= m;
        }
    }
    public static ulong Pochhammer(this uint value, uint n) { var r = 1uL; while (n > 0u) r *= value + --n; return r; }
    public static long Pochhammer(this int value, uint n) { var r = 1L; while (n > 0u) r *= value + --n; return r; }
    public static ulong Factorial(this uint value) { var r = 1uL; while (value > 1u) r *= value--; return r; }
    public static ulong Factorial(this ushort value) { var r = 1uL; while (value > 1u) r *= value--; return r; }
    public static ulong Factorial(this byte value) { var r = 1uL; while (value > 1u) r *= value--; return r; }
    public static ulong Combination(this uint n, uint k) { return (n - k + 1u).Pochhammer(k) / k.Factorial(); }
    public static double Combination(this uint n, double k) { return (n - k + 1d).Pochhammer(k) / (++k).Gamma(); }
}