Fixed DiGammaInv bug for 0 and infinity arguments.

Added more more special functions. Signed-off-by: jvangael <jurgen.vangael@gmail.com> Signed-off-by: jvangael <jurgen.vangael@gmail.com>
17 years ago · f5f81a11ab
5 changed files with 167 additions and 3 deletions
--- a/src/Numerics/SpecialFunctions/Erf.cs
+++ b/src/Numerics/SpecialFunctions/Erf.cs
@ -112,6 +112,49 @@ namespace MathNet.Numerics

            return ErfImp(x, true);
        }
+        ///<summary>Calculates the inverse error function evaluated at z.</summary>
+        /// <returns>The inverse error function evaluated at given value.</returns>
+        /// <remarks>
+        /// 	<list type="bullet">
+        /// 		<item>returns Double.PositiveInfinity if <c>z &gt;= 1.0</c>.</item>
+        /// 		<item>returns Double.NegativeInfinity if <c>z &lt;= -1.0</c>.</item>
+        /// 	</list>
+        /// </remarks>
+        ///<summary>Calculates the inverse error function evaluated at z.</summary>
+        ///<param name="z">value to evaluate.</param>
+        ///<returns>the inverse error function evaluated at Z.</returns>
+        public static double ErfInv(double z)
+        {
+            if (z == 0.0)
+            {
+                return 0.0;
+            }
+
+            if (z >= 1.0)
+            {
+                return double.PositiveInfinity;
+            }
+            if (z <= -1.0)
+            {
+                return double.NegativeInfinity;
+            }
+
+            double p, q, s;
+            if (z < 0)
+            {
+                p = -z;
+                q = 1 - p;
+                s = -1;
+            }
+            else
+            {
+                p = z;
+                q = 1 - z;
+                s = 1;
+            }
+
+            return ErfInvImpl(p, q, s);
+        }

        /// <summary>
        /// Implementation of the error function.
--- a/src/Numerics/SpecialFunctions/Factorial.cs
+++ b/src/Numerics/SpecialFunctions/Factorial.cs
@ -29,6 +29,7 @@
 namespace MathNet.Numerics
 {
    using System;
+    using Properties;

    public partial class SpecialFunctions
    {
@ -127,5 +128,47 @@ namespace MathNet.Numerics

            return FactorialLn(n) - FactorialLn(k) - FactorialLn(n - k);
        }
+
+        /// <summary>
+        /// Computes the multinomial coefficient: n choose n1, n2, n3, ...
+        /// </summary>
+        /// <param name="n">A nonnegative value n.</param>
+        /// <param name="ni">An array of nonnegative values that sum to <paramref name="n"/>.</param>
+        /// <returns>The multinomial coefficient.</returns>
+        /// <exception cref="ArgumentNullException">if <paramref name="ni"/> is <see langword="null" />.</exception>   
+        /// <exception cref="ArgumentException">If <paramref name="n"/> or any of the <paramref name="ni"/> are negative.</exception>
+        /// <exception cref="ArgumentException">If the sum of all <paramref name="ni"/> is not equal to <paramref name="n"/>.</exception>
+        public static double Multinomial(int n, int[] ni)
+        {
+            if (n < 0)
+            {
+                throw new ArgumentException(Resources.ArgumentMustBePositive, "n");
+            }
+            if (ni == null)
+            {
+                throw new ArgumentNullException("ni");
+            }
+
+            int sum = 0;
+            double ret = FactorialLn(n);
+            for (int i = 0; i < ni.Length; i++)
+            {
+                if (ni[i] < 0)
+                {
+                    throw new ArgumentException(Resources.ArgumentMustBePositive, "ni[" + i + "]");
+                }
+
+                ret -= FactorialLn(ni[i]);
+                sum += ni[i];
+            }
+
+            // Before returning, check that the sum of all elements was equal to n.
+            if (sum != n)
+            {
+                throw new ArgumentException(Resources.ArgumentParameterSetInvalid , "ni");
+            }
+
+            return System.Math.Floor(0.5 + System.Math.Exp(ret));
+        }
    }
 }
--- a/src/UnitTests/SpecialFunctionsTests/ErfTests.cs
+++ b/src/UnitTests/SpecialFunctionsTests/ErfTests.cs
@ -26,7 +26,7 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>

-namespace MathNet.Numerics.UnitTests.SpecialFunctionTests
+namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
 {
    using MbUnit.Framework;
    using MathNet.Numerics;
@ -106,5 +106,29 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionTests
        {
            AssertHelpers.AlmostEqual(f, SpecialFunctions.ErfcInv(x), 8);
        }
+
+        [Test]
+        [Row(double.NaN, double.NaN)]
+        [Row(-1.0, -0.84270079294971486934122063508260925929606699796630291)]
+        [Row(0.0, 0.0)]
+        [Row(1e-15, 0.0000000000000011283791670955126615773132947717431253912942469337536)]
+        [Row(0.1, 0.1124629160182848984047122510143040617233925185058162)]
+        [Row(0.2, 0.22270258921047846617645303120925671669511570710081967)]
+        [Row(0.3, 0.32862675945912741618961798531820303325847175931290341)]
+        [Row(0.4, 0.42839235504666847645410962730772853743532927705981257)]
+        [Row(0.5, 0.5204998778130465376827466538919645287364515757579637)]
+        [Row(1.0, 0.84270079294971486934122063508260925929606699796630291)]
+        [Row(1.5, 0.96610514647531072706697626164594785868141047925763678)]
+        [Row(2.0, 0.99532226501895273416206925636725292861089179704006008)]
+        [Row(2.5, 0.99959304798255504106043578426002508727965132259628658)]
+        [Row(3.0, 0.99997790950300141455862722387041767962015229291260075)]
+        [Row(4.0, 0.99999998458274209971998114784032651311595142785474641)]
+        [Row(5.0, 0.99999999999846254020557196514981165651461662110988195)]
+        [Row(double.PositiveInfinity, 1.0)]
+        [Row(double.NegativeInfinity, -1.0)]
+        public void ErfInv(double x, double f)
+        {
+            AssertHelpers.AlmostEqual(x, SpecialFunctions.ErfInv(f), 6);
+        }
    }
 }
--- a/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs
+++ b/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs
@ -102,5 +102,14 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
            AssertHelpers.AlmostEqual(Math.Log(0), SpecialFunctions.BinomialLn(5, 7), 14);
            AssertHelpers.AlmostEqual(Math.Log(0), SpecialFunctions.BinomialLn(5, -7), 14);
        }
+
+        [Test]
+        public void CanComputeMultinomial()
+        {
+            AssertHelpers.AlmostEqual(1, SpecialFunctions.Multinomial(1, new int[] { 1, 0 }), 14);
+            AssertHelpers.AlmostEqual(10, SpecialFunctions.Multinomial(5, new int[] { 3, 2 }), 14);
+            AssertHelpers.AlmostEqual(10, SpecialFunctions.Multinomial(5, new int[] { 2, 3 }), 14);
+            AssertHelpers.AlmostEqual(35, SpecialFunctions.Multinomial(7, new int[] { 3, 4 }), 14);
+        }
    }
 }
--- a/src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs
+++ b/src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs
@ -1,4 +1,4 @@
-// <copyright file="SpecialFunctionsTests.cs" company="Math.NET">
+// <copyright file="SpecialFunctionsTests.cs" company="Math.NET">
 // Math.NET Numerics, part of the Math.NET Project
 // http://mathnet.opensourcedotnet.info
 //
@ -26,7 +26,7 @@
 // OTHER DEALINGS IN THE SOFTWARE.
 // </copyright>

-namespace MathNet.Numerics.UnitTests.SpecialFunctionTests
+namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
 {
    using System;
    using MbUnit.Framework;
@ -103,6 +103,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionTests

        [Test]
        [Row(Double.NaN, Double.NaN)]
+        [Row(0.0, Double.NegativeInfinity)]
        [Row(0.1, -10.423754940411076232100295314502760886768558023951363)]
        [Row(1.0, -0.57721566490153286060651209008240243104215933593992359)]
        [Row(1.5, 0.036489973978576520559023667001244432806840395339565888)]
@ -117,9 +118,53 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionTests
        [Row(5.0, 1.5061176684318004727268212432509309022911739973934097)]
        [Row(5.5, 1.6110931485817511237336268416044190359814435699427405)]
        [Row(10.1, 2.2622143570941481235561593642219403924532310597356171)]
+        [Row(Double.PositiveInfinity, Double.PositiveInfinity)]
        public void DiGammaInv(double x, double f)
        {
            AssertHelpers.AlmostEqual(x, SpecialFunctions.DiGammaInv(f), 13);
        }
+
+        /// <summary>
+        /// Compute the t'th harmonic number using a loop.
+        /// </summary>
+        private double ExactHarmonic(int t)
+        {
+            double r = 0.0;
+            for (int i = 1; i <= t; i++)
+            {
+                r += 1.0 / i;
+            }
+            return r;
+        }
+
+        [Test]
+        [Row(1)]
+        [Row(2)]
+        [Row(4)]
+        [Row(8)]
+        [Row(16)]
+        [Row(100)]
+        [Row(1000)]
+        [Row(10000)]
+        [Row(100000)]
+        [Row(1000000)]
+        public void Harmonic(int i)
+        {
+            AssertHelpers.AlmostEqual(ExactHarmonic(i), SpecialFunctions.Harmonic(i), 13);
+        }
+
+        [Test]
+        public void BetaLn()
+        {
+            AssertHelpers.AlmostEqual(System.Math.Log(0.5), SpecialFunctions.BetaLn(1.0, 2.0), 14);
+            AssertHelpers.AlmostEqual(System.Math.Log(1.0), SpecialFunctions.BetaLn(1.0, 1.0), 14);
+        }
+
+        [Test]
+        public void Beta()
+        {
+            AssertHelpers.AlmostEqual(0.5, SpecialFunctions.Beta(1.0, 2.0), 14);
+            AssertHelpers.AlmostEqual(1.0, SpecialFunctions.Beta(1.0, 1.0), 14);
+        }
    }
 }