From 3db0e003af6d407a02bda28fe2e8ee10e15f9b54 Mon Sep 17 00:00:00 2001
From: Thorsten Sommer <mail@tsommer.org>
Date: Sat, 31 Oct 2020 12:19:02 +0100
Subject: [PATCH] Added math functions: gamma and factorial

---
 FastRng/Double/MathTools.cs           |  61 +++++
 FastRngTests/Double/MathToolsTests.cs | 341 ++++++++++++++++++++++++++
 2 files changed, 402 insertions(+)
 create mode 100644 FastRng/Double/MathTools.cs
 create mode 100644 FastRngTests/Double/MathToolsTests.cs

diff --git a/FastRng/Double/MathTools.cs b/FastRng/Double/MathTools.cs
new file mode 100644
index 0000000..3de59d6
--- /dev/null
+++ b/FastRng/Double/MathTools.cs
@@ -0,0 +1,61 @@
+using System;
+
+namespace FastRng.Double
+{
+    public static class MathTools
+    {
+        private static readonly double SQRT_2 = Math.Sqrt(2.0);
+        private static readonly double SQRT_PI = Math.Sqrt(Math.PI);
+        
+        public static double Gamma(double z)
+        {
+            // Source: http://rosettacode.org/wiki/Gamma_function#Go
+            
+            const double F1 = 6.5;
+            const double A1 = .99999999999980993;
+            const double A2 = 676.5203681218851;
+            const double A3 = 1259.1392167224028;
+            const double A4 = 771.32342877765313;
+            const double A5 = 176.61502916214059;
+            const double A6 = 12.507343278686905;
+            const double A7 = .13857109526572012;
+            const double A8 = 9.9843695780195716e-6;
+            const double A9 = 1.5056327351493116e-7;
+
+            var t = z + F1;
+            var x =  A1 +
+                     A2 / z -
+                     A3 / (z + 1) +
+                     A4 / (z + 2) -
+                     A5 / (z + 3) +
+                     A6 / (z + 4) -
+                     A7 / (z + 5) +
+                     A8 / (z + 6) +
+                     A9 / (z + 7);
+
+            return MathTools.SQRT_2 * MathTools.SQRT_PI * Math.Pow(t, z - 0.5) * Math.Exp(-t) * x;
+        }
+        
+        public static double Factorial(double x) => MathTools.Gamma(x + 1.0);
+
+        public static ulong Factorial(uint x)
+        {
+            if (x > 20)
+                throw new ArgumentOutOfRangeException(nameof(x), $"Cannot compute {x}!, since ulong.max is 18_446_744_073_709_551_615.");
+            
+            ulong accumulator = 1;
+            for (uint factor = 1; factor <= x; factor++)
+                accumulator *= factor;
+
+            return accumulator;
+        }
+
+        public static ulong Factorial(int x)
+        {
+            if(x < 0)
+                throw new ArgumentOutOfRangeException(nameof(x), "Given value must be greater as zero.");
+
+            return MathTools.Factorial((uint) x);
+        }
+    }
+}
\ No newline at end of file
diff --git a/FastRngTests/Double/MathToolsTests.cs b/FastRngTests/Double/MathToolsTests.cs
new file mode 100644
index 0000000..85f7bff
--- /dev/null
+++ b/FastRngTests/Double/MathToolsTests.cs
@@ -0,0 +1,341 @@
+using System;
+using System.Diagnostics.CodeAnalysis;
+using FastRng.Double;
+using NUnit.Framework;
+
+namespace FastRngTests.Double
+{
+    [ExcludeFromCodeCoverage]
+    public class MathToolsTests
+    {
+        #region Gamma
+
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void GammaTest01()
+        {
+            Assert.That(MathTools.Gamma(-0.5), Is.EqualTo(-3.544907701811087).Within(1e-6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void GammaTest02()
+        {
+            Assert.That(MathTools.Gamma(0.1), Is.EqualTo(9.513507698668752).Within(1e-6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void GammaTest03()
+        {
+            Assert.That(MathTools.Gamma(0.5), Is.EqualTo(1.772453850905517).Within(1e-6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void GammaTest04()
+        {
+            Assert.That(MathTools.Gamma(1.0), Is.EqualTo(1.0).Within(1e-6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void GammaTest05()
+        {
+            Assert.That(MathTools.Gamma(1.5), Is.EqualTo(0.8862269254527587).Within(1e-6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void GammaTest06()
+        {
+            Assert.That(MathTools.Gamma(2.0), Is.EqualTo(1.0).Within(1e-6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void GammaTest07()
+        {
+            Assert.That(MathTools.Gamma(3.0), Is.EqualTo(2.0).Within(1e-6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void GammaTest08()
+        {
+            Assert.That(MathTools.Gamma(10.0), Is.EqualTo(362_880.0000000015).Within(1e-6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void GammaTest09()
+        {
+            Assert.That(MathTools.Gamma(140.0), Is.EqualTo(9.6157231969402357e+238).Within(1e-6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void GammaTest10()
+        {
+            Assert.That(MathTools.Gamma(170.0), Is.EqualTo(double.PositiveInfinity));
+        }
+
+        #endregion
+
+        #region Factorial (integer)
+
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger01()
+        {
+            Assert.That(MathTools.Factorial(0), Is.EqualTo(1));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger02()
+        {
+            Assert.That(MathTools.Factorial(1), Is.EqualTo(1));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger03()
+        {
+            Assert.That(MathTools.Factorial(2), Is.EqualTo(2));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger04()
+        {
+            Assert.That(MathTools.Factorial(3), Is.EqualTo(6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger05()
+        {
+            Assert.That(MathTools.Factorial(4), Is.EqualTo(24));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger06()
+        {
+            Assert.That(MathTools.Factorial(5), Is.EqualTo(120));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger07()
+        {
+            Assert.That(MathTools.Factorial(6), Is.EqualTo(720));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger08()
+        {
+            Assert.That(MathTools.Factorial(7), Is.EqualTo(5_040));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger09()
+        {
+            Assert.That(MathTools.Factorial(8), Is.EqualTo(40_320));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger10()
+        {
+            Assert.That(MathTools.Factorial(9), Is.EqualTo(362_880));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger11()
+        {
+            Assert.That(MathTools.Factorial(10), Is.EqualTo(3_628_800));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger12()
+        {
+            Assert.That(MathTools.Factorial(11), Is.EqualTo(39_916_800));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger13()
+        {
+            Assert.That(MathTools.Factorial(12), Is.EqualTo(479_001_600));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger14()
+        {
+            Assert.That(MathTools.Factorial(13), Is.EqualTo(6_227_020_800));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger15()
+        {
+            Assert.That(MathTools.Factorial(14), Is.EqualTo(87_178_291_200));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger16()
+        {
+            Assert.That(MathTools.Factorial(15), Is.EqualTo(1_307_674_368_000));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger17()
+        {
+            Assert.That(MathTools.Factorial(16), Is.EqualTo(20_922_789_888_000));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger18()
+        {
+            Assert.That(MathTools.Factorial(17), Is.EqualTo(355_687_428_096_000));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger19()
+        {
+            Assert.That(MathTools.Factorial(18), Is.EqualTo(6_402_373_705_728_000));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger20()
+        {
+            Assert.That(MathTools.Factorial(19), Is.EqualTo(121_645_100_408_832_000));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger21()
+        {
+            Assert.That(MathTools.Factorial(20), Is.EqualTo(2_432_902_008_176_640_000));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger22()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(() => MathTools.Factorial(21));
+
+            // Note: 21! is not possible in C# until we got 128 bit integers, since:
+            //       ulong.max == 18_446_744_073_709_551_615 < 51_090_942_171_709_400_000
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger23()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(() => MathTools.Factorial(45_646));
+
+            // Note: 45_646! is not possible in C# since:
+            //       ulong.max == 18_446_744_073_709_551_615
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger24()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(() => MathTools.Factorial(-1));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialInteger25()
+        {
+            Assert.Throws<ArgumentOutOfRangeException>(() => MathTools.Factorial(-6_565));
+        }
+
+        #endregion
+
+        #region Factorial (floating point)
+
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialFloatingPoint01()
+        {
+            Assert.That(MathTools.Factorial(0.5), Is.EqualTo(0.886226925).Within(1e6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialFloatingPoint02()
+        {
+            Assert.That(MathTools.Factorial(1.5), Is.EqualTo(1.329340388).Within(1e6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialFloatingPoint03()
+        {
+            Assert.That(MathTools.Factorial(-1.5), Is.EqualTo(-1.329340388).Within(1e6));
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public void FactorialFloatingPoint04()
+        {
+            Assert.That(MathTools.Factorial(7.5), Is.EqualTo(14_034.407293483).Within(1e6));
+        }
+
+        #endregion
+    }
+}
\ No newline at end of file