diff --git a/src/Mathematics.NET/Mathematics.cs b/src/Mathematics.NET/Mathematics.cs index 52404f69..e9d8d688 100644 --- a/src/Mathematics.NET/Mathematics.cs +++ b/src/Mathematics.NET/Mathematics.cs @@ -25,6 +25,10 @@ // SOFTWARE. // +using System.Numerics; +using System.Runtime.CompilerServices; +using Complex = Mathematics.NET.Core.Complex; + namespace Mathematics.NET; /// Provides Mathematics.NET functionality @@ -84,4 +88,89 @@ public static class Mathematics /// public static Real ZetaOf4 => Real.ZetaOf4; + + // + // Methods + // + + // Integer-related + + /// The binomial coefficient + /// This method works with BigInteger. Be aware that this method may take a long time. + /// A type that implements + /// A total of elements + /// An unordered subset of elements + /// A cancellation token + /// choose + /// + public static T Binomial(T n, T k, CancellationToken cancellationToken) + where T : IBinaryInteger + { + var result = T.One; + checked + { + T u = T.One; + for (T i = T.Zero; i < k; i++) + { + cancellationToken.ThrowIfCancellationRequested(); + result *= n - i; + u *= i + T.One; + + var gcd = GCD(result, u); + if (gcd != T.One) + { + result /= gcd; + u /= gcd; + } + } + } + return result; + } + + /// Compute the greatest common divisor of two integers. + /// A type that implements + /// An integer + /// An integer + /// The GCD of the two values + [MethodImpl(MethodImplOptions.AggressiveInlining)] + public static T GCD(T p, T q) + where T : IBinaryInteger + { + p = T.Abs(p); + q = T.Abs(q); + while (p != T.Zero && q != T.Zero) + { + if (p > q) + { + p %= q; + } + else + { + q %= p; + } + } + return p | q; + } + + /// The multinomial coefficient + /// An array of positive integers + /// A cancellation token + /// The coefficient associated with a term with powers given by + /// + public static T Multinomial(ReadOnlySpan α, CancellationToken cancellationToken) + where T : IBinaryInteger + { + var result = T.One; + checked + { + var partialSum = α[0]; + for (int i = 1; i < α.Length; i++) + { + var element = α[i]; + partialSum += element; + result *= Binomial(partialSum, element, cancellationToken); + } + } + return result; + } }