functions: Factorial, FactorialLn

Signed-off-by: Christoph Ruegg <git@cdrnet.ch>

src/Numerics/Numerics.csproj

@@ -108,6 +108,7 @@
     <Compile Include="Sorting.cs" />
     <Compile Include="SpecialFunctions.cs" />
     <Compile Include="SpecialFunctions\Erf.cs" />
+    <Compile Include="SpecialFunctions\Factorial.cs" />
     <Compile Include="Statistics\DescriptiveStatistics.cs" />
     <Compile Include="Statistics\Statistics.cs" />
     <Compile Include="Threading\AggregateException.cs" />

src/Numerics/SpecialFunctions/Factorial.cs

@@ -0,0 +1,111 @@
+// <copyright file="Factorial.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+namespace MathNet.Numerics
+{
+    using System;
+
+    public partial class SpecialFunctions
+    {
+        /// <summary>
+        /// 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.
+        /// </summary>
+        /// <returns>A value value! for value > 0</returns>
+        /// <remarks>
+        /// If you need to multiply or divide various such factorials, consider using the logarithmic version 
+        /// <see cref="FactorialLn"/> instead so you can add instead of multiply and subtract instead of divide, and
+        /// then exponentiate the result using <see cref="System.Math.Exp"/>. This will also circumvent the problem that
+        /// factorials become very large even for small parameters.
+        /// </remarks>
+        /// <exception cref="ArgumentOutOfRangeException" />
+        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;
+        }
+
+        /// <summary>
+        /// Computes the logarithmic factorial function x -> ln(x!) of an integer number > 0.
+        /// </summary>
+        /// <returns>A value value! for value > 0</returns>
+        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;
+    }
+}

src/UnitTests/SpecialFunctionsTests/FactorialTest.cs

@@ -0,0 +1,82 @@
+// <copyright file="FactorialTest.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+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<ArgumentOutOfRangeException>(() => SpecialFunctions.Factorial(Int32.MinValue));
+            Assert.Throws<ArgumentOutOfRangeException>(() => SpecialFunctions.Factorial(-1));
+            Assert.Throws<ArgumentOutOfRangeException>(() => 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);
+        }
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -93,8 +93,9 @@
     <Compile Include="PrecisionTest.cs" />
     <Compile Include="Properties\AssemblyInfo.cs" />
     <Compile Include="SortingTests.cs" />
-    <Compile Include="SpecialFunctionsTest\ErfTests.cs" />
+    <Compile Include="SpecialFunctionsTests\ErfTests.cs" />
-    <Compile Include="SpecialFunctionsTest\SpecialFunctionsTests.cs" />
+    <Compile Include="SpecialFunctionsTests\FactorialTest.cs" />
+    <Compile Include="SpecialFunctionsTests\SpecialFunctionsTests.cs" />
     <Compile Include="StatisticsTests\DescriptiveStatisticsTests.cs" />
     <Compile Include="StatisticsTests\StatisticsTests.cs" />
     <Compile Include="StatisticsTests\StatTestData.cs" />