diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 145b3afd..f240ced1 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -108,6 +108,7 @@ + diff --git a/src/Numerics/SpecialFunctions/Factorial.cs b/src/Numerics/SpecialFunctions/Factorial.cs new file mode 100644 index 00000000..6807062a --- /dev/null +++ b/src/Numerics/SpecialFunctions/Factorial.cs @@ -0,0 +1,111 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics +{ + using System; + + public partial class SpecialFunctions + { + /// + /// Computes the factorial function x -> x! of an integer number > 0. The function can represent all number up + /// to 22! exactly, all numbers up to 170! using a double representation. All larger values will overflow. + /// + /// A value value! for value > 0 + /// + /// If you need to multiply or divide various such factorials, consider using the logarithmic version + /// instead so you can add instead of multiply and subtract instead of divide, and + /// then exponentiate the result using . This will also circumvent the problem that + /// factorials become very large even for small parameters. + /// + /// + public static double Factorial(int x) + { + if (x < 0) + { + throw new ArgumentOutOfRangeException("x", Properties.Resources.ArgumentPositive); + } + + if (x <= FactorialMaxArgument) + { + if (factorialCache == null) + { + factorialCache = GenerateFactorials(FactorialMaxArgument); + } + + return factorialCache[x]; + } + + return Double.PositiveInfinity; + } + + /// + /// Computes the logarithmic factorial function x -> ln(x!) of an integer number > 0. + /// + /// A value value! for value > 0 + public static double FactorialLn(int x) + { + if (x < 0) + { + throw new ArgumentOutOfRangeException("x", Properties.Resources.ArgumentPositive); + } + + if (x <= 1) + { + return 0d; + } + + if (x <= FactorialMaxArgument) + { + if (factorialCache == null) + { + factorialCache = GenerateFactorials(FactorialMaxArgument); + } + + return Math.Log(factorialCache[x]); + } + + return GammaLn(x + 1.0); + } + + private static double[] GenerateFactorials(int max) + { + var cache = new double[max + 1]; + cache[0] = 1.0; + for (int i = 1; i < cache.Length; i++) + { + cache[i] = cache[i - 1] * i; + } + + return cache; + } + + private const int FactorialMaxArgument = 170; + private static double[] factorialCache; + } +} \ No newline at end of file diff --git a/src/UnitTests/SpecialFunctionsTest/ErfTests.cs b/src/UnitTests/SpecialFunctionsTests/ErfTests.cs similarity index 100% rename from src/UnitTests/SpecialFunctionsTest/ErfTests.cs rename to src/UnitTests/SpecialFunctionsTests/ErfTests.cs diff --git a/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs b/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs new file mode 100644 index 00000000..38dc6448 --- /dev/null +++ b/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs @@ -0,0 +1,82 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests +{ + using System; + using MbUnit.Framework; + + [TestFixture] + public class FactorialTest + { + [Test] + public void CanComputeFactorial() + { + // exact + double factorial = 1.0; + for (int i = 1; i < 23; i++) + { + factorial *= i; + AssertHelpers.AlmostEqual(factorial, SpecialFunctions.Factorial(i), 14); + AssertHelpers.AlmostEqual(Math.Log(factorial), SpecialFunctions.FactorialLn(i), 14); + } + + // approximation + for (int i = 23; i < 171; i++) + { + factorial *= i; + AssertHelpers.AlmostEqual(factorial, SpecialFunctions.Factorial(i), 14); + AssertHelpers.AlmostEqual(Math.Log(factorial), SpecialFunctions.FactorialLn(i), 14); + } + } + + [Test] + public void ThrowsOnNegativeArgument() + { + Assert.Throws(() => SpecialFunctions.Factorial(Int32.MinValue)); + Assert.Throws(() => SpecialFunctions.Factorial(-1)); + Assert.Throws(() => SpecialFunctions.FactorialLn(-1)); + } + + [Test] + public void FactorialOverflowsToInfinity() + { + Assert.AreEqual(Double.PositiveInfinity, SpecialFunctions.Factorial(172)); + Assert.AreEqual(Double.PositiveInfinity, SpecialFunctions.Factorial(Int32.MaxValue)); + } + + [Test] + public void FactorialLnDoesNotOverflow() + { + AssertHelpers.AlmostEqual(6078.2118847500501140, SpecialFunctions.FactorialLn(1 << 10), 14); + AssertHelpers.AlmostEqual(29978.648060844048236, SpecialFunctions.FactorialLn(1 << 12), 14); + AssertHelpers.AlmostEqual(307933.81973375485425, SpecialFunctions.FactorialLn(1 << 15), 14); + AssertHelpers.AlmostEqual(1413421.9939462073242, SpecialFunctions.FactorialLn(1 << 17), 14); + } + } +} \ No newline at end of file diff --git a/src/UnitTests/SpecialFunctionsTest/SpecialFunctionsTests.cs b/src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs similarity index 100% rename from src/UnitTests/SpecialFunctionsTest/SpecialFunctionsTests.cs rename to src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 42972add..c3364370 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -93,8 +93,9 @@ - - + + +