Special functions use auto-generated truncated series (#207)

* Python tool to generate expressions for truncated series.
* Using the generated series in special functions
* Test projects set Optimize=true in Release configuration

Co-authored-by: Tom Minka <8955276+tminka@users.noreply.github.com>

Infer.sln

@@ -34,17 +34,17 @@ Project("{6EC3EE1D-3C4E-46DD-8F32-0CC8E7565705}") = "TestFSharp", "test\TestFSha
 EndProject
 Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Examples", "Examples", "{DC5F5BC4-CDB0-41F7-8B03-CD4C38C8DEB2}"
 EndProject
-Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "ClickThroughModel", "src\Examples\ClickThroughModel\ClickThroughModel.csproj", "{33D86EA2-2161-4EF0-8F17-59602296273C}"
+Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "ClickThroughModel", "src\Examples\ClickThroughModel\ClickThroughModel.csproj", "{33D86EA2-2161-4EF0-8F17-59602296273C}"
 EndProject
-Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "ClinicalTrial", "src\Examples\ClinicalTrial\ClinicalTrial.csproj", "{B517BBF2-60E6-4C69-885A-AE5C014D877B}"
+Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "ClinicalTrial", "src\Examples\ClinicalTrial\ClinicalTrial.csproj", "{B517BBF2-60E6-4C69-885A-AE5C014D877B}"
 EndProject
 Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "InferNET101", "src\Examples\InferNET101\InferNET101.csproj", "{52D174E7-2407-4FC1-9DDA-4D9D14F18618}"
 EndProject
-Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "MontyHall", "src\Examples\MontyHall\MontyHall.csproj", "{6139CF19-0190-4ED5-AEE3-D3CE7458E517}"
+Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "MontyHall", "src\Examples\MontyHall\MontyHall.csproj", "{6139CF19-0190-4ED5-AEE3-D3CE7458E517}"
 EndProject
 Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "MotifFinder", "src\Examples\MotifFinder\MotifFinder.csproj", "{D2A7B5F5-8D33-45AC-9776-07C23F5859BB}"
 EndProject
-Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Image_Classifier", "src\Examples\ImageClassifier\Image_Classifier.csproj", "{87D09BD4-119E-49C1-B0B4-86DF962A00EE}"
+Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Image_Classifier", "src\Examples\ImageClassifier\Image_Classifier.csproj", "{87D09BD4-119E-49C1-B0B4-86DF962A00EE}"
 EndProject
 Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "LDA", "src\Examples\LDA\LDA.csproj", "{6FF3E672-378C-4D61-B4CA-A5A5E01C2563}"
 EndProject
@@ -90,6 +90,8 @@ Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "RobustGaussianProcess", "sr
 EndProject
 Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "LearnersNuGet", "src\Learners\LearnersNuGet\LearnersNuGet.csproj", "{6BC792EE-0436-494E-9F0C-63A3822320D0}"
 EndProject
+Project("{888888A0-9F3D-457C-B088-3A5042F75D52}") = "Tools.GenerateSeries", "src\Tools\GenerateSeries\Tools.GenerateSeries.pyproj", "{D7562CC3-D48A-481D-B40C-07EE51EBBF50}"
+EndProject
 Global
 	GlobalSection(SolutionConfigurationPlatforms) = preSolution
 		Debug|Any CPU = Debug|Any CPU
@@ -498,6 +500,12 @@ Global
 		{6BC792EE-0436-494E-9F0C-63A3822320D0}.ReleaseCore|Any CPU.Build.0 = Release|Any CPU
 		{6BC792EE-0436-494E-9F0C-63A3822320D0}.ReleaseFull|Any CPU.ActiveCfg = Release|Any CPU
 		{6BC792EE-0436-494E-9F0C-63A3822320D0}.ReleaseFull|Any CPU.Build.0 = Release|Any CPU
+		{D7562CC3-D48A-481D-B40C-07EE51EBBF50}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{D7562CC3-D48A-481D-B40C-07EE51EBBF50}.DebugCore|Any CPU.ActiveCfg = Debug|Any CPU
+		{D7562CC3-D48A-481D-B40C-07EE51EBBF50}.DebugFull|Any CPU.ActiveCfg = Debug|Any CPU
+		{D7562CC3-D48A-481D-B40C-07EE51EBBF50}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{D7562CC3-D48A-481D-B40C-07EE51EBBF50}.ReleaseCore|Any CPU.ActiveCfg = Release|Any CPU
+		{D7562CC3-D48A-481D-B40C-07EE51EBBF50}.ReleaseFull|Any CPU.ActiveCfg = Release|Any CPU
 	EndGlobalSection
 	GlobalSection(SolutionProperties) = preSolution
 		HideSolutionNode = FALSE

src/Runtime/Core/Maths/Series.cs

@@ -1,503 +0,0 @@
-using System;
-using System.Collections;
-using System.Collections.Generic;
-using System.Text;
-using System.Linq;
-
-namespace Microsoft.ML.Probabilistic.Core.Maths
-{
-    /// <summary>
-    /// Represents a truncated power series (effectively, a polynomial) of the form
-    /// c0 + c1*x + c2*x^2 + c3*x^3 + ... + cN*x^N
-    /// </summary>
-    struct TruncatedPowerSeries
-    {
-        /// <summary>
-        /// Coefficients of the truncated power series.
-        /// Coefficients[k] is the coefficient near x^k, so Coefficients[0] is the constant term.
-        /// </summary>
-        public double[] Coefficients { get; }
-
-        /// <summary>
-        /// Constructs a truncated power series.
-        /// </summary>
-        /// <param name="coefficients">
-        /// Coefficients of the truncated power series.
-        /// <paramref name="coefficients"/>[k] is the coefficient near x^k, so <paramref name="coefficients"/>[0] is the constant term.
-        /// </param>
-        public TruncatedPowerSeries(double[] coefficients)
-        {
-            this.Coefficients = coefficients;
-        }
-
-        /// <summary>
-        /// Factory function to generate a truncated power series.
-        /// All coefficients of the generated series are computed, when this function is called.
-        /// </summary>
-        /// <param name="length">
-        /// Number of coefficients including the constant term to generate.
-        /// The highest power of x in the resulting series will be (<paramref name="length"/> - 1)
-        /// </param>
-        /// <param name="coefGenerator">
-        /// Provides the coefficients in the series, where the constant term is index zero.
-        /// Coefficients[k] = <paramref name="coefGenerator"/>(k)
-        /// </param>
-        /// <returns>Generated truncated power series.</returns>
-        public static TruncatedPowerSeries Generate(int length, Func<int, double> coefGenerator)
-        {
-            double[] coefficients = new double[length];
-            for (int i = 0; i < length; ++i)
-                coefficients[i] = coefGenerator(i);
-
-            return new TruncatedPowerSeries(coefficients);
-        }
-
-        /// <summary>
-        /// Factory function to generate a truncated power series.
-        /// All coefficients of the generated series are computed, when this function is called.
-        /// Coefficients with smaller indexes are used to compute coefficients with bigger indexes.
-        /// </summary>
-        /// <param name="length">
-        /// Number of coefficients including the constant term to generate.
-        /// The highest power of x in the resulting series will be (<paramref name="length"/> - 1)
-        /// </param>
-        /// <param name="coefGenerator">
-        /// Provides the coefficients in the series, where the constant term is index zero.
-        /// The second argument is the coefficients with indexes smaller than that of the coefficient being computed.
-        /// Coefficients[k] = <paramref name="coefGenerator"/>(k, Coefficients[0..k-1])
-        /// </param>
-        /// <returns>Generated truncated power series.</returns>
-        public static TruncatedPowerSeries Generate(int length, Func<int, IList<double>, double> coefGenerator)
-        {
-            List<double> coefficients = new List<double>(length);
-            for (int i = 0; i < length; ++i)
-                coefficients.Add(coefGenerator(i, coefficients));
-
-            return new TruncatedPowerSeries(coefficients.ToArray());
-        }
-
-        public double Evaluate(double x)
-        {
-            double sum = 0.0;
-            for (int k = Coefficients.Length - 1; k > 0; k--)
-            {
-                sum = x * (Coefficients[k] + sum);
-            }
-            return sum + Coefficients[0];
-        }
-    }
-
-    /// <summary>
-    /// Represents a mathematical power series of the form c0 + c1*x + c2*x^2 + c3*x^3 + ...
-    /// </summary>
-    struct PowerSeries
-    {
-        readonly Func<int, double> coefGenerator;
-
-        /// <summary>
-        /// Constructs a power series.
-        /// </summary>
-        /// <param name="coefGenerator">Provides the coefficients in the series, where the constant term is index zero. 
-        /// Every zero coefficient must eventually be followed by a nonzero coefficient.</param>
-        public PowerSeries(Func<int, double> coefGenerator)
-        {
-            this.coefGenerator = coefGenerator;
-        }
-
-        public double Evaluate(double x)
-        {
-            double sum = coefGenerator(0);
-            double term = 1;
-            for (int i = 1; ; i++)
-            {
-                double oldSum = sum;
-                term *= x;
-                double coefficient = coefGenerator(i);
-                // To ensure termination, every zero coefficient must eventually be followed by a nonzero coefficient.
-                if (coefficient != 0)
-                {
-                    sum += term * coefficient;
-                    if (sum == oldSum) break;
-                }
-            }
-            return sum;
-        }
-    }
-
-    class Series
-    {
-        public TruncatedPowerSeries GammaAt2 { get; }
-        public TruncatedPowerSeries DigammaAt2 { get; }
-        /// <summary>
-        /// Coefficients of de Moivre's expansion for the digamma function.
-        /// </summary>
-        public TruncatedPowerSeries DigammaAsymptotic { get; }
-        public TruncatedPowerSeries TrigammaAt1 { get; }
-        /// <summary>
-        /// Coefficients of de Moivre's expansion for the trigamma function.
-        /// </summary>
-        public TruncatedPowerSeries TrigammaAsymptotic { get; }
-        public TruncatedPowerSeries TetragammaAt1 { get; }
-        /// <summary>
-        /// Coefficients of de Moivre's expansion for the tetragamma function.
-        /// </summary>
-        public TruncatedPowerSeries TetragammaAsymptotic { get; }
-        public TruncatedPowerSeries GammalnAsymptotic { get; }
-        public PowerSeries Log1Plus { get; }
-        public PowerSeries Log1Minus { get; }
-
-        public Series(uint precisionBits)
-        {
-            uint doublePrecisionCutOff = 53;
-            if (precisionBits <= doublePrecisionCutOff)
-            {
-                GammaAt2 = new TruncatedPowerSeries(gammaAt2Coefficients.Take(26).ToArray());
-                DigammaAsymptotic = new TruncatedPowerSeries(digammaAsymptoticCoefficients.Take(8).ToArray());
-                TrigammaAt1 = new TruncatedPowerSeries(trigammaAt1Coefficients.Take(2).ToArray());
-                TrigammaAsymptotic = new TruncatedPowerSeries(trigammaAsymptoticCoefficients.Take(9).ToArray());
-                TetragammaAt1 = new TruncatedPowerSeries(tetragammaAt1Coefficients.Take(2).ToArray());
-                TetragammaAsymptotic = new TruncatedPowerSeries(tetragammaAsymptoticCoefficients.Take(9).ToArray());
-                GammalnAsymptotic = new TruncatedPowerSeries(gammalnAsymptoticCoefficients.Take(7).ToArray());
-            }
-            else
-            {
-                GammaAt2 = new TruncatedPowerSeries(gammaAt2Coefficients.ToArray());
-                DigammaAsymptotic = new TruncatedPowerSeries(digammaAsymptoticCoefficients.ToArray());
-                TrigammaAt1 = new TruncatedPowerSeries(trigammaAt1Coefficients.ToArray());
-                TrigammaAsymptotic = new TruncatedPowerSeries(trigammaAsymptoticCoefficients.ToArray());
-                TetragammaAt1 = new TruncatedPowerSeries(tetragammaAt1Coefficients.ToArray());
-                TetragammaAsymptotic = new TruncatedPowerSeries(tetragammaAsymptoticCoefficients.ToArray());
-                GammalnAsymptotic = new TruncatedPowerSeries(gammalnAsymptoticCoefficients.ToArray());
-            }
-            DigammaAt2 = TruncatedPowerSeries.Generate(GammaAt2.Coefficients.Length, i => GammaAt2.Coefficients[i] * (i + 1));
-            Log1Plus = new PowerSeries(Log1PlusCoefficient);
-            Log1Minus = new PowerSeries(Log1MinusCoefficient);
-        }
-
-        #region Coefficient generators
-
-        private static double Log1PlusCoefficient(int n) => n == 0 ? 0.0 : (n % 2 == 0 ? -1.0 : 1.0) / n;
-        private static double Log1MinusCoefficient(int n) => n == 0 ? 0.0 : -1.0 / n;
-
-        #endregion
-
-        #region Precomputed coefficients
-
-        /* Python code to generate this table (must not be indented):
-from __future__ import division
-from sympy import *
-print("            0.0,")
-for k in range(2,110):
-		print("            %0.34g," % ((-1)**k*(zeta(k)-1)/k).evalf(34))
-         */
-        private static readonly double[] gammaAt2Coefficients = new[]
-        {
-            0,
-            0.3224670334241132030328458313306328,
-            -0.06735230105319810200992236559613957,
-            0.02058080842778454641606167285772244,
-            -0.00738555102867398567678680620929299,
-            0.002890510330741523359332489917505882,
-            -0.001192753911703261018861788045342109,
-            0.000509669524743042450139196564151689,
-            -0.0002231547584535793857934971029521876,
-            9.945751278180853098033475934158787e-05,
-            -4.492623673813314204598598489148742e-05,
-            2.050721277567069106674904621634425e-05,
-            -9.439488275268396715186013101739348e-06,
-            4.37486678990748817440030460357292e-06,
-            -2.039215753801366189692763700169742e-06,
-            9.551412130407419353005918014054565e-07,
-            -4.492469198764566185487696111516831e-07,
-            2.120718480555466464491199763200968e-07,
-            -1.004322482396809908404802422730862e-07,
-            4.769810169363980398286696825264908e-08,
-            -2.271109460894316350425416834932224e-08,
-            1.083865921489695459289913246009665e-08,
-            -5.183475041970046644229946111422438e-09,
-            2.483674543802478475235497043321041e-09,
-            -1.192140140586091154743934012838306e-09,
-            5.731367241678862251438860640198256e-10,
-            -2.759522885124233355903938000146247e-10,
-            1.330476437424448882003567968775656e-10,
-            -6.422964563838099598864779601491062e-11,
-            3.104424774732227563392725117909481e-11,
-            -1.502138408075414166458594899693007e-11,
-            7.275974480239079174553105387739685e-12,
-            -3.527742476575915065189557262694141e-12,
-            1.711991790559617978544962723971016e-12,
-            -8.315385841420284979314085891293852e-13,
-            4.042200525289440192260810475387441e-13,
-            -1.966475631096616531130082285066565e-13,
-            9.573630387838555566047107836916075e-14,
-            -4.664076026428374443204720937282992e-14,
-            2.273736960065972417354758289435552e-14,
-            -1.109139947083452217591392010101267e-14,
-            5.413659156725363290762788546531038e-15,
-            -2.643880017860994856070349482131646e-15,
-            1.291895906278996649797589034450644e-15,
-            -6.315935504198448152815144632969046e-16,
-            3.089316266963393006209215246728354e-16,
-            -1.511793062810819824159605773901573e-16,
-            7.401486856952319873829631842748242e-17,
-            -3.625218048120653818830185998027823e-17,
-            1.776356842186163324845409513060544e-17,
-            -8.707631574791790545720149092818225e-18,
-            4.270088559227003785123993337350986e-18,
-            -2.094760424794464294026626618224464e-18,
-            1.027984282378792796690424930004006e-18,
-            -5.046468294792953328756217103439225e-19,
-            2.478176394593791693371463023660605e-19,
-            -1.217349807814763706484975874824611e-19,
-            5.981805089941246261735805969233367e-20,
-            -2.940209281436570304360268231936827e-20,
-            1.445602896686655580463396514538343e-20,
-            -7.109522442656804701132425037724369e-21,
-            3.497426362898741485992060556339303e-21,
-            -1.72095582935593863487067784712846e-21,
-            8.470329472588509127992257400804905e-22,
-            -4.170008355728413714644388723288791e-22,
-            2.053413205469873327785584183708884e-22,
-            -1.011382623588834223218438654968751e-22,
-            4.982546748559995181315757428368296e-23,
-            -2.455167963057679834508841583052901e-23,
-            1.210047067506714006249752230839652e-23,
-            -5.965020755313850001484134344489508e-24,
-            2.94108662241138167013854591546344e-24,
-            -1.450398882284963546831250522314545e-24,
-            7.153994486945770419577690483337377e-25,
-            -3.529303946893137333637237751583067e-25,
-            1.741432868532761836030085313783932e-25,
-            -8.594084286265459638986861745574907e-26,
-            4.241951859246373681570791434569361e-26,
-            -2.094128133045665479257451510604071e-26,
-            1.033975765691293117071248028646776e-26,
-            -5.106053163907605937631813496430958e-27,
-            2.521892111442166601145541317372597e-27,
-            -1.245753934567815698050968063924673e-27,
-            6.154617652924322580570211265420193e-28,
-            -3.041105193209663666084926529947082e-28,
-            1.502871752458263590864301387762562e-28,
-            -7.42798682249486384616203701175441e-29,
-            3.671788940665074353595069096429321e-29,
-            -1.815266442575991720790924205550424e-29,
-            8.975484077181291618744867280323167e-30,
-            -4.438426192012726108791433766298579e-30,
-            2.195091214528033139239927553086061e-30,
-            -1.085744041594510958603950819611772e-30,
-            5.370967865334548331937709775945041e-31,
-            -2.657215680744460758305032942851813e-31,
-            1.314768175368352936035252591887088e-31,
-            -6.506069321410406518422414169431261e-32,
-            3.2198404294735178635568008439894e-32,
-            -1.593658394385882570638484064799577e-32,
-            7.888609052210118218331949440892363e-33,
-            -3.905252006044612945499590526006017e-33,
-            1.93348261083581322172569684496361e-33,
-            -9.57355467501227836151640760828371e-34,
-            4.740750632337811593763855621594469e-34,
-            -2.347800313157773442176833948742271e-34,
-            1.162825626799840470674685835501613e-34,
-            -5.759790487887060065479999327489828e-35,
-            2.853229547240349253548725242810521e-35,
-            -1.413526564687879603645186309191011e-35,
-        };
-
-        /* Coefficients of de Moivre's expansion for the digamma function.
-         Each coefficient is B_{2j}/(2j) where B_{2j} are the Bernoulli numbers, starting from j=1
-         Python code to generate this table (must not be indented):
-from __future__ import division
-from sympy import *
-print("            0,")
-for k in range(1,22):
-		print("            " + str(bernoulli(2 * k) / (2 * k)).replace("/", ".0 / ") + ".0,")
-        */
-        private static readonly double[] digammaAsymptoticCoefficients = new[]
-        {
-            0.0,
-            1.0 / 12.0,
-            -1.0 / 120.0,
-            1.0 / 252.0,
-            -1.0 / 240.0,
-            1.0 / 132.0,
-            -691.0 / 32760.0,
-            1.0 / 12.0,
-            -3617.0 / 8160.0,
-            43867.0 / 14364.0,
-            -174611.0 / 6600.0,
-            77683.0 / 276.0,
-            //-236364091.0 / 65520.0,
-            //657931.0 / 12.0,
-            //-3392780147.0 / 3480.0,
-            //1723168255201.0 / 85932.0,
-            //-7709321041217.0 / 16320.0,
-            //151628697551.0 / 12.0,
-            //-26315271553053477373.0 / 69090840.0,
-            //154210205991661.0 / 12.0,
-            //-261082718496449122051.0 / 541200.0,
-            //1520097643918070802691.0 / 75852.0,
-        };
-
-        /* Python code to generate this table (must not be indented):
-from __future__ import division
-from sympy import *
-for k in range(0,20):
-		print("            %0.34g," % ((-1)**k * (k + 1) * zeta(k + 2)).evalf(34))
-         */
-        private static readonly double[] trigammaAt1Coefficients = new[]
-        {
-            1.644934066848226406065691662661266,
-            -2.40411380631918847328165611543227,
-            3.246969701133414432092649803962559,
-            -4.147711020573479956397022760938853,
-            5.086715309922245964457943045999855,
-            -6.050095664291537111978414031909779,
-            7.028541493385610294808429898694158,
-            -8.016067142608656936886291077826172,
-            9.0089511761503633380243627470918,
-            -10.00494188604119472074671648442745,
-            11.00270695208638827011782268527895,
-            -12.00147256017094221647312224376947,
-            13.00079622575576365761662600561976,
-            -14.00042823530829849687506793998182,
-            15.00022923389112960990132705774158,
-            -16.00012219516220568493736209347844,
-            17.00006489398550613145744137000293,
-            -18.00003434782889755183532543014735,
-            19.00001812527864331059390679001808,
-            -20.00000953865973585266146983485669,
-        };
-
-        /* Coefficients of de Moivre's expansion for the digamma function.
-         Each coefficient is B_{2j} where B_{2j} are the Bernoulli numbers, starting from j=1
-         Python code to generate this table (must not be indented):
-from __future__ import division
-from sympy import *
-print("            0,")
-for k in range(1,22):
-		print("            " + str(bernoulli(2 * k)).replace("/", ".0 / ") + ".0,")
-        */
-        private static readonly double[] trigammaAsymptoticCoefficients = new[]
-        {
-            0.0,
-            1.0 / 6.0,
-            -1.0 / 30.0,
-            1.0 / 42.0,
-            -1.0 / 30.0,
-            5.0 / 66.0,
-            -691.0 / 2730.0,
-            7.0 / 6.0,
-            -3617.0 / 510.0,
-            43867.0 / 798.0,
-            //-174611.0 / 330.0,
-            //854513.0 / 138.0,
-            //-236364091.0 / 2730.0,
-            //8553103.0 / 6.0,
-            //-23749461029.0 / 870.0,
-            //8615841276005.0 / 14322.0,
-            //-7709321041217.0 / 510.0,
-            //2577687858367.0 / 6.0,
-            //-26315271553053477373.0 / 1919190.0,
-            //2929993913841559.0 / 6.0,
-            //-261082718496449122051.0 / 13530.0,
-            //1520097643918070802691.0 / 1806.0,
-        };
-
-        /* Python code to generate this table (must not be indented):
-from __future__ import division
-from sympy import *
-for k in range(0,20):
-		print("            %0.34g," % ((-1)**(k + 1) * (k + 1) * (k + 2) * zeta(k + 3)).evalf(34))
-         */
-        private static readonly double[] tetragammaAt1Coefficients = new[]
-        {
-            -2.40411380631918847328165611543227,
-            6.493939402266828864185299607925117,
-            -12.44313306172043986919106828281656,
-            20.34686123968898385783177218399942,
-            -30.25047832145768467171365045942366,
-            42.17124896031366176885057939216495,
-            -56.11246999826060743998823454603553,
-            72.0716094092029067041949019767344,
-            -90.04447697437075248672044835984707,
-            110.0270695208638898066055844537914,
-            -132.0161981618803679339180234819651,
-            156.0095547090691638913995120674372,
-            -182.0055670590078875648032408207655,
-            210.0032092744758074331912212073803,
-            -240.0018329274330994849151466041803,
-            272.0010383037680981033190619200468,
-            -306.0005839130912477230594959110022,
-            342.0003262550155795906903222203255,
-            -380.0001812345349776478542480617762,
-            420.0001001492111640800430905073881,
-        };
-
-        /* Coefficients of de Moivre's expansion for the digamma function.
-         Each coefficient is -(2j+1) B_{2j} where B_{2j} are the Bernoulli numbers, starting from j=0
-         Python code to generate this table (must not be indented):
-from __future__ import division
-from sympy import *
-print("            0.0,")
-for k in range(0,21):
-		print("            " + str(-(2 * k + 1) * bernoulli(2 * k)).replace("/", ".0 / ") + ".0,")
-        */
-        private static readonly double[] tetragammaAsymptoticCoefficients = new[]
-        {
-            0.0,
-            -1.0,
-            -1.0 / 2.0,
-            1.0 / 6.0,
-            -1.0 / 6.0,
-            3.0 / 10.0,
-            -5.0 / 6.0,
-            691.0 / 210.0,
-            -35.0 / 2.0,
-            //3617.0 / 30.0,
-            //-43867.0 / 42.0,
-            //1222277.0 / 110.0,
-            //-854513.0 / 6.0,
-            //1181820455.0 / 546.0,
-            //-76977927.0 / 2.0,
-            //23749461029.0 / 30.0,
-            //-8615841276005.0 / 462.0,
-            //84802531453387.0 / 170.0,
-            //-90219075042845.0 / 6.0,
-            //26315271553053477373.0 / 51870.0,
-            //-38089920879940267.0 / 2.0,
-            //261082718496449122051.0 / 330.0,
-        };
-
-        /* Python code to generate this table (must not be indented):
-from __future__ import division
-from sympy import *
-for k in range(1,21):
-		print("            " + str(bernoulli(2 * k) / (2 * k * (2 * k - 1))).replace("/", ".0 / ") + ".0,")
-        */
-        private static readonly double[] gammalnAsymptoticCoefficients = new[]
-        {
-            1.0 / 12.0,
-            -1.0 / 360.0,
-            1.0 / 1260.0,
-            -1.0 / 1680.0,
-            1.0 / 1188.0,
-            -691.0 / 360360.0,
-            1.0 / 156.0,
-            -3617.0 / 122400.0,
-            43867.0 / 244188.0,
-            -174611.0 / 125400.0,
-            77683.0 / 5796.0,
-            //-236364091.0 / 1506960.0,
-            //657931.0 / 300.0,
-            //-3392780147.0 / 93960.0,
-            //1723168255201.0 / 2492028.0,
-            //-7709321041217.0 / 505920.0,
-            //151628697551.0 / 396.0,
-            //-26315271553053477373.0 / 2418179400.0,
-            //154210205991661.0 / 444.0,
-            //-261082718496449122051.0 / 21106800.0,
-        };
-
-        #endregion
-    }
-}

src/Runtime/Core/Maths/SpecialFunctions.cs

@@ -36,27 +36,6 @@ namespace Microsoft.ML.Probabilistic.Math
     /// </remarks>
     public static class MMath
     {
-        internal static uint NumericalPrecisionBits { get; private set; } = 53;
-
-        internal static void SetNumericalPrecision(uint value)
-        {
-            if (NumericalPrecisionBits != value)
-            {
-                if (double.IsNaN(value) || value <= 0)
-                    throw new ArgumentException(nameof(value));
-
-                var newSeries = new Series(value);
-                lock(precisionSettingLockObject)
-                {
-                    NumericalPrecisionBits = value;
-                    Series = newSeries;
-                }
-            }
-        }
-
-        static object precisionSettingLockObject = new object();
-        static Series Series { get; set; } = new Series(NumericalPrecisionBits);
-
         #region Bessel functions
 
         /// <summary>
@@ -766,7 +745,35 @@ namespace Microsoft.ML.Probabilistic.Math
                 // Use Taylor series at x=2
                 // Reference: https://dlmf.nist.gov/5.7#E3
                 double dx = x - 2;
-                double sum = Series.GammaAt2.Evaluate(dx);
+                double sum =
+                    // Truncated series 1: Gamma at 2
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    dx * (0.322467033424113218236207583323 +
+                    dx * (-0.0673523010531980951332460538371 +
+                    dx * (0.0205808084277845478790009241353 +
+                    dx * (-0.00738555102867398526627309729141 +
+                    dx * (0.00289051033074152328575298829849 +
+                    dx * (-0.00119275391170326097711393569283 +
+                    dx * (0.000509669524743042422335654813582 +
+                    dx * (-0.000223154758453579379761418803601 +
+                    dx * (0.0000994575127818085337145958900319 +
+                    dx * (-0.0000449262367381331417002075024064 +
+                    dx * (0.0000205072127756706915531665039783 +
+                    dx * (-0.00000943948827526839590398742510441 +
+                    dx * (0.00000437486678990748780418179322395 +
+                    dx * (-0.00000203921575380136623678190070967 +
+                    dx * (0.000000955141213040741983285717977295 +
+                    dx * (-0.000000449246919876456604329429033119 +
+                    dx * (0.000000212071848055546658692313590108 +
+                    dx * (-0.000000100432248239680996087208305005 +
+                    dx * (0.0000000476981016936398056576019341725 +
+                    dx * (-0.0000000227110946089431649103199811606 +
+                    dx * (0.0000000108386592148969540910749175797 +
+                    dx * (-0.00000000518347504197004665512124864706 +
+                    dx * (0.00000000248367454380247831718500866399 +
+                    dx * (-0.00000000119214014058609120744254820277 +
+                    dx * 5.73136724167886201333019485796e-10))))))))))))))))))))))))
+                    ;
                 sum = dx * (1 + Digamma1 + sum);
                 result += sum;
                 return result;
@@ -875,10 +882,15 @@ namespace Microsoft.ML.Probabilistic.Math
             // The threshold for applying de Moivre's expansion for the digamma function.
             const double c_digamma_small = 1e-6;
 
-            /* Use Taylor series if argument <= S */
+            /* Shift the argument and use Taylor series near 1 if argument <= S */
             if (x <= c_digamma_small)
             {
-                return (Digamma1 - 1 / x + Zeta2 * x);
+                return - 1 / x +
+                    // Truncated series 2: Digamma at 1
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    -0.577215664901532860606512090082 +
+                    x * 1.64493406684822643647241516665
+                    ;
             }
 
             if (x <= 2.5)
@@ -890,7 +902,35 @@ namespace Microsoft.ML.Probabilistic.Math
                     x++;
                 }
                 double dx = x - 2;
-                double sum2 = Series.DigammaAt2.Evaluate(dx);
+                double sum2 =
+                    // Truncated series 3: Digamma at 2
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    dx * (0.644934066848226436472415166646 +
+                    dx * (-0.202056903159594285399738161511 +
+                    dx * (0.0823232337111381915160036965412 +
+                    dx * (-0.0369277551433699263313654864570 +
+                    dx * (0.0173430619844491397145179297909 +
+                    dx * (-0.00834927738192282683979754984980 +
+                    dx * (0.00407735619794433937868523850865 +
+                    dx * (-0.00200839282608221441785276923241 +
+                    dx * (0.000994575127818085337145958900319 +
+                    dx * (-0.000494188604119464558702282526470 +
+                    dx * (0.000246086553308048298637998047740 +
+                    dx * (-0.000122713347578489146751836526357 +
+                    dx * (0.0000612481350587048292585451051353 +
+                    dx * (-0.0000305882363070204935517285106451 +
+                    dx * (0.0000152822594086518717325714876367 +
+                    dx * (-0.00000763719763789976227360029356303 +
+                    dx * (0.00000381729326499983985646164462194 +
+                    dx * (-0.00000190821271655393892565695779510 +
+                    dx * (0.000000953962033872796113152038683449 +
+                    dx * (-0.000000476932986787806463116719604373 +
+                    dx * (0.000000238450502727732990003648186753 +
+                    dx * (-0.000000119219925965311073067788718882 +
+                    dx * (0.0000000596081890512594796124402079358 +
+                    dx * (-0.0000000298035035146522801860637050694 +
+                    dx * 0.0000000149015548283650412346585066307))))))))))))))))))))))))
+                    ;
                 result2 += sum2;
                 return result2;
             }
@@ -910,7 +950,17 @@ namespace Microsoft.ML.Probabilistic.Math
             double invX = 1 / x;
             result += Math.Log(x) - 0.5 * invX;
             double invX2 = invX * invX;
-            double sum = Series.DigammaAsymptotic.Evaluate(invX2);
+            double sum =
+                // Truncated series 4: Digamma asymptotic
+                // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                invX2 * (1.0 / 12.0 +
+                invX2 * (-1.0 / 120.0 +
+                invX2 * (1.0 / 252.0 +
+                invX2 * (-1.0 / 240.0 +
+                invX2 * (1.0 / 132.0 +
+                invX2 * (-691.0 / 32760.0 +
+                invX2 * 1.0 / 12.0))))))
+                ;
             result -= sum;
             return result;
         }
@@ -947,7 +997,12 @@ namespace Microsoft.ML.Probabilistic.Math
             /* Shift the argument and use Taylor series at 1 if argument <= small */
             if (x <= c_trigamma_small)
             {
-                return (1.0 / (x * x) + Series.TrigammaAt1.Evaluate(x));
+                return 1.0 / (x * x) +
+                    // Truncated series 5: Trigamma at 1
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    1.64493406684822643647241516665 +
+                    x * -2.40411380631918857079947632302
+                    ;
             }
 
             result = 0.0;
@@ -963,7 +1018,18 @@ namespace Microsoft.ML.Probabilistic.Math
             // This expansion can be computed in Maple via asympt(Psi(1,x),x)
             double invX2 = 1 / (x * x);
             result += 0.5 * invX2;
-            double sum = Series.TrigammaAsymptotic.Evaluate(invX2);
+            double sum =
+                    // Truncated series 6: Trigamma asymptotic
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    invX2 * (1.0 / 6.0 +
+                    invX2 * (-1.0 / 30.0 +
+                    invX2 * (1.0 / 42.0 +
+                    invX2 * (-1.0 / 30.0 +
+                    invX2 * (5.0 / 66.0 +
+                    invX2 * (-691.0 / 2730.0 +
+                    invX2 * (7.0 / 6.0 +
+                    invX2 * -3617.0 / 510.0)))))))
+                    ;
             result += (1 + sum) / x;
             return result;
         }
@@ -983,7 +1049,12 @@ namespace Microsoft.ML.Probabilistic.Math
                          c_tetragamma_small = 1e-4;
             /* Use Taylor series if argument <= small */
             if (x < c_tetragamma_small)
-                return -2 / (x * x * x) + Series.TetragammaAt1.Evaluate(x);
+                return -2 / (x * x * x) +
+                    // Truncated series 7: Tetragamma at 1
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    -2.40411380631918857079947632302 +
+                    x * 6.49393940226682914909602217925
+                    ;
             double result = 0;
             /* Reduce to Tetragamma(x+n) where ( X + N ) >= L */
             while (x < c_tetragamma_large)
@@ -996,7 +1067,18 @@ namespace Microsoft.ML.Probabilistic.Math
             // Milton Abramowitz and Irene A. Stegun, Handbook of Mathematical Functions, Section 6.4
             double invX2 = 1 / (x * x);
             result += -invX2 / x;
-            double sum = Series.TetragammaAsymptotic.Evaluate(invX2);
+            double sum =
+                // Truncated series 8: Tetragamma asymptotic
+                // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                invX2 * (-1.0 +
+                invX2 * (-1.0 / 2.0 +
+                invX2 * (1.0 / 6.0 +
+                invX2 * (-1.0 / 6.0 +
+                invX2 * (3.0 / 10.0 +
+                invX2 * (-5.0 / 6.0 +
+                invX2 * (691.0 / 210.0 +
+                invX2 * -35.0 / 2.0)))))))
+                ;
             result += sum;
             return result;
         }
@@ -1185,8 +1267,15 @@ namespace Microsoft.ML.Probabilistic.Math
             double phi;
             if (Math.Abs(xOverAMinus1) < 1e-1)
             {
-                double XMinusLog1PlusCoefficient(int n) => (n <= 1) ? 0.0 : (n % 2 == 0 ? 1.0 : -1.0) / n;
-                phi = new PowerSeries(XMinusLog1PlusCoefficient).Evaluate(xOverAMinus1);
+                phi =
+                    // Truncated series 12: x - log(1 + x)
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    xOverAMinus1 * xOverAMinus1 * (1.0 / 2.0 +
+                    xOverAMinus1 * (-1.0 / 3.0 +
+                    xOverAMinus1 * (1.0 / 4.0 +
+                    xOverAMinus1 * (-1.0 / 5.0 +
+                    xOverAMinus1 * 1.0 / 6.0))))
+                    ;
             }
             else
             {
@@ -1329,40 +1418,28 @@ namespace Microsoft.ML.Probabilistic.Math
         {
             if (x > 0.3)
                 return 1 / MMath.Gamma(x + 1) - 1;
-            double sum = 0;
-            double term = x;
-            for (int i = 0; i < reciprocalFactorialMinus1Coeffs.Length; i++)
-            {
-                sum += reciprocalFactorialMinus1Coeffs[i] * term;
-                term *= x;
-            }
-            return sum;
-        }
 
-        /* using http://sagecell.sagemath.org/ (must not be indented)
-var('x');
-f = 1/gamma(x+1)-1
-[c[0].n(100) for c in f.taylor(x,0,16).coefficients()]
-         */
-        private static readonly double[] reciprocalFactorialMinus1Coeffs =
-        {
-            0.57721566490153286060651209008,
-            -0.65587807152025388107701951515,
-            -0.042002635034095235529003934875,
-            0.16653861138229148950170079510,
-            -0.042197734555544336748208301289,
-            -0.0096219715278769735621149216721,
-            0.0072189432466630995423950103404,
-            -0.0011651675918590651121139710839,
-            -0.00021524167411495097281572996312,
-            0.00012805028238811618615319862640,
-            -0.000020134854780788238655689391375,
-            -1.2504934821426706573453594884e-6,
-            1.1330272319816958823741296196e-6,
-            -2.0563384169776071034501526118e-7,
-            6.1160951044814158178625767024e-9,
-            5.0020076444692229300557063872e-9
-        };
+            return
+                // Truncated series 18: Reciprocal factorial minus 1
+                // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                x * (0.577215664901532860606512090082 +
+                x * (-0.655878071520253881077019515145 +
+                x * (-0.0420026350340952355290039348754 +
+                x * (0.166538611382291489501700795102 +
+                x * (-0.0421977345555443367482083012891 +
+                x * (-0.00962197152787697356211492167231 +
+                x * (0.00721894324666309954239501034044 +
+                x * (-0.00116516759185906511211397108402 +
+                x * (-0.000215241674114950972815729963027 +
+                x * (0.000128050282388116186153198626338 +
+                x * (-0.0000201348547807882386556893913662 +
+                x * (-0.00000125049348214267065734535945301 +
+                x * (0.00000113302723198169588237412961344 +
+                x * (-0.000000205633841697760710345015364430 +
+                x * (0.00000000611609510448141581786252410878 +
+                x * 0.00000000500200764446922293005568923640)))))))))))))))
+                ;
+        }
 
         /// <summary>
         /// Computes <c>x^a e^(-x)/Gamma(a)</c> to high accuracy.
@@ -1446,7 +1523,17 @@ f = 1/gamma(x+1)-1
                 // the series is:  sum_{i=1}^inf B_{2i} / (2i*(2i-1)*x^(2i-1))
                 double invX = 1.0 / x;
                 double invX2 = invX * invX;
-                double sum = invX * Series.GammalnAsymptotic.Evaluate(invX2);
+                double sum = invX * (
+                    // Truncated series 9: GammaLn asymptotic
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    1.0 / 12.0 +
+                    invX2 * (-1.0 / 360.0 +
+                    invX2 * (1.0 / 1260.0 +
+                    invX2 * (-1.0 / 1680.0 +
+                    invX2 * (1.0 / 1188.0 +
+                    invX2 * (-691.0 / 360360.0 +
+                    invX2 * 1.0 / 156.0)))))
+                    );
                 return sum;
             }
         }
@@ -1782,7 +1869,14 @@ f = 1/gamma(x+1)-1
                 // log(NormalCdf(x)) = log(1-NormalCdf(-x))
                 double y = NormalCdf(-x);
                 // log(1-y)
-                return -y * (1 + y * (0.5 + y * (1.0 / 3 + y * (0.25))));
+                return
+                    // Truncated series 11: log(1 - x)
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    y * (-1.0 +
+                    y * (-1.0 / 2.0 +
+                    y * (-1.0 / 3.0 +
+                    y * -1.0 / 4.0)))
+                    ;
             }
             else if (x >= large)
             {
@@ -1794,13 +1888,17 @@ f = 1/gamma(x+1)-1
                 double z = 1 / (x * x);
                 // asymptotic series for log(normcdf)
                 // Maple command: subs(x=-x,asympt(log(normcdf(-x)),x));
-                double s = z * (c_normcdfln_series[0] + z *
-                              (c_normcdfln_series[1] + z *
-                               (c_normcdfln_series[2] + z *
-                                (c_normcdfln_series[3] + z *
-                                 (c_normcdfln_series[4] + z *
-                                  (c_normcdfln_series[5] + z *
-                                   c_normcdfln_series[6]))))));
+                double s =
+                    // Truncated series 16: normcdfln asymptotic
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    z * (-1.0 +
+                    z * (5.0 / 2.0 +
+                    z * (-37.0 / 3.0 +
+                    z * (353.0 / 4.0 +
+                    z * (-4081.0 / 5.0 +
+                    z * (55205.0 / 6.0 +
+                    z * -854197.0 / 7.0))))))
+                    ;
                 return s - LnSqrt2PI - 0.5 * x * x - Math.Log(-x);
             }
         }
@@ -3114,7 +3212,14 @@ rr = mpf('-0.99999824265582826');
             if (x > 1e-4)
                 return 1 - Math.Sqrt(1 - x);
             else
-                return x * (1.0 / 2 + x * (1.0 / 8 + x * (1.0 / 16 + x * 5.0 / 128)));
+                return
+                    // Truncated series 17: 1 - sqrt(1 - x)
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    x * (1.0 / 2.0 +
+                    x * (1.0 / 8.0 +
+                    x * (1.0 / 16.0 +
+                    x * 5.0 / 128.0)))
+                    ;
         }
 
         /// <summary>
@@ -3251,7 +3356,22 @@ rr = mpf('-0.99999824265582826');
             {
                 // use the Taylor series for log(1+x) around x=0
                 // Maple command: series(log(1+x),x);
-                return Series.Log1Plus.Evaluate(x);
+                return
+                    // Truncated series 10: log(1 + x)
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    x * (1.0 +
+                    x * (-1.0 / 2.0 +
+                    x * (1.0 / 3.0 +
+                    x * (-1.0 / 4.0 +
+                    x * (1.0 / 5.0 +
+                    x * (-1.0 / 6.0 +
+                    x * (1.0 / 7.0 +
+                    x * (-1.0 / 8.0 +
+                    x * (1.0 / 9.0 +
+                    x * (-1.0 / 10.0 +
+                    x * (1.0 / 11.0 +
+                    x * -1.0 / 12.0)))))))))))
+                    ;
             }
             else
             {
@@ -3304,7 +3424,20 @@ rr = mpf('-0.99999824265582826');
             if (x < -3.5)
             {
                 double expx = Math.Exp(x);
-                return Series.Log1Minus.Evaluate(expx);
+                return
+                    // Truncated series 11: log(1 - x)
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    expx * (-1.0 +
+                    expx * (-1.0 / 2.0 +
+                    expx * (-1.0 / 3.0 +
+                    expx * (-1.0 / 4.0 +
+                    expx * (-1.0 / 5.0 +
+                    expx * (-1.0 / 6.0 +
+                    expx * (-1.0 / 7.0 +
+                    expx * (-1.0 / 8.0 +
+                    expx * (-1.0 / 9.0 +
+                    expx * -1.0 / 10.0)))))))))
+                    ;
             }
             else
             {
@@ -3325,7 +3458,14 @@ rr = mpf('-0.99999824265582826');
         {
             if (Math.Abs(x) < 2e-3)
             {
-                return x * (1 + x * (0.5 + x * (1.0 / 6 + x * (1.0 / 24))));
+                return
+                    // Truncated series 13: exp(x) - 1
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    x * (1.0 +
+                    x * (1.0 / 2.0 +
+                    x * (1.0 / 6.0 +
+                    x * 1.0 / 24.0)))
+                    ;
             }
             else
             {
@@ -3342,9 +3482,22 @@ rr = mpf('-0.99999824265582826');
         {
             if (Math.Abs(x) < 6e-1)
             {
-                return x * (1.0 / 6 + x * (1.0 / 24 + x * (1.0 / 120 + x * (1.0 / 720 +
-                    x * (1.0 / 5040 + x * (1.0 / 40320 + x * (1.0 / 362880 + x * (1.0 / 3628800 +
-                    x * (1.0 / 39916800 + x * (1.0 / 479001600 + x * (1.0 / 6227020800 + x * (1.0 / 87178291200))))))))))));
+                return
+                    // Truncated series 14: ((exp(x) - 1) / x - 1) / x - 0.5
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    x * (1.0 / 6.0 +
+                    x * (1.0 / 24.0 +
+                    x * (1.0 / 120.0 +
+                    x * (1.0 / 720.0 +
+                    x * (1.0 / 5040.0 +
+                    x * (1.0 / 40320.0 +
+                    x * (1.0 / 362880.0 +
+                    x * (1.0 / 3628800.0 +
+                    x * (1.0 / 39916800.0 +
+                    x * (1.0 / 479001600.0 +
+                    x * (1.0 / 6227020800.0 +
+                    x * 1.0 / 87178291200.0)))))))))))
+                    ;
             }
             else if (double.IsPositiveInfinity(x))
             {
@@ -3370,7 +3523,13 @@ rr = mpf('-0.99999824265582826');
         {
             if (x < 1e-3)
             {
-                return Math.Log(x) + x * (0.5 + x * (1.0 / 24 + x * (-1.0 / 2880)));
+                return Math.Log(x) +
+                    // Truncated series 15: log(exp(x) - 1) / x
+                    // Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py
+                    x * (1.0 / 2.0 +
+                    x * (1.0 / 24.0 +
+                    x * x * -1.0 / 2880.0))
+                    ;
             }
             else if (x > 50)
             {
@@ -4646,14 +4805,6 @@ else if (m < 20.0 - 60.0/11.0 * s) {
         /// </summary>
         private const double Zeta4 = 1.0823232337111381915;
 
-        /// <summary>
-        /// Coefficients of the asymptotic expansion of NormalCdfLn.
-        /// </summary>
-        private static readonly double[] c_normcdfln_series =
-            {
-                -1, 5.0/2, -37.0/3, 353.0/4, -4081.0/5, 55205.0/6, -854197.0/7
-            };
-
         /// <summary>
         /// NormCdf(x)/NormPdf(x) for x = 0, -1, -2, -3, ..., -16
         /// </summary>

src/Tools/GenerateSeries/GenerateSeries.py

@@ -0,0 +1,328 @@
+# Licensed to the .NET Foundation under one or more agreements.
+# The .NET Foundation licenses this file to you under the MIT license.
+# See the LICENSE file in the project root for more information.
+from __future__ import division
+from sympy import zeta, evalf, bernoulli, symbols, Poly, series, factorial, factorial2, S, log, exp, gamma, digamma, sqrt
+
+# configuration
+decimal_precision = 30
+evalf_inner_precision = 500
+
+gamma_at_2_series_length = 26
+gamma_at_2_variable_name = "dx"
+gamma_at_2_indent = "                    "
+
+digamma_at_1_series_length = 2
+digamma_at_1_variable_name = "x"
+digamma_at_1_indent = "                    "
+
+digamma_at_2_series_length = 26
+digamma_at_2_variable_name = "dx"
+digamma_at_2_indent = "                    "
+
+digamma_asymptotic_series_length = 8
+digamma_asymptotic_variable_name = "invX2"
+digamma_asymptotic_indent = "                "
+
+trigamma_at_1_series_length = 2
+trigamma_at_1_variable_name = "x"
+trigamma_at_1_indent = "                    "
+
+trigamma_asymptotic_series_length = 9
+trigamma_asymptotic_variable_name = "invX2"
+trigamma_asymptotic_indent = "                    "
+
+tetragamma_at_1_series_length = 2
+tetragamma_at_1_variable_name = "x"
+tetragamma_at_1_indent = "                    "
+
+tetragamma_asymptotic_series_length = 9
+tetragamma_asymptotic_variable_name = "invX2"
+tetragamma_asymptotic_indent = "                "
+
+gammaln_asymptotic_series_length = 7
+gammaln_asymptotic_variable_name = "invX2"
+gammaln_asymptotic_indent = "                    "
+
+log_1_plus_series_length = 13
+log_1_plus_variable_name = "x"
+log_1_plus_indent = "                    "
+
+log_1_minus_series_length = 11
+log_1_minus_variable_name = "expx"
+log_1_minus_indent = "                    "
+
+x_minus_log_1_plus_series_length = 7
+x_minus_log_1_plus_variable_name = "xOverAMinus1"
+x_minus_log_1_plus_indent = "                    "
+
+exp_minus_1_series_length = 5
+exp_minus_1_variable_name = "x"
+exp_minus_1_indent = "                    "
+
+exp_minus_1_ratio_minus_1_ratio_minus_half_series_length = 13
+exp_minus_1_ratio_minus_1_ratio_minus_half_variable_name = "x"
+exp_minus_1_ratio_minus_1_ratio_minus_half_indent = "                    "
+
+log_exp_minus_1_ratio_series_length = 5
+log_exp_minus_1_ratio_variable_name = "x"
+log_exp_minus_1_ratio_indent = "                    "
+
+normcdfln_asymptotic_series_length = 8
+normcdfln_asymptotic_variable_name = "z"
+normcdfln_asymptotic_indent = "                    "
+
+one_minus_sqrt_one_minus_series_length = 5
+one_minus_sqrt_one_minus_variable_name = "x"
+one_minus_sqrt_one_minus_indent = "                    "
+
+reciprocal_factorial_minus_1_series_length = 17
+reciprocal_factorial_minus_1_variable_name = "x"
+reciprocal_factorial_minus_1_indent = "                "
+
+def print_heading_comment(indent, header):
+    print(f"{indent}// Truncated series {header}")
+    print(f"{indent}// Generated automatically by /src/Tools/GenerateSeries/GenerateSeries.py")
+
+def format_real_coefficient(coefficient):
+    return str(coefficient)
+
+def print_polynomial_with_real_coefficients(varname, coefficients, indent):
+    if len(coefficients) <= 1:
+        print(f"{indent}{format_real_coefficient(coefficients[0])}")
+        return
+    if coefficients[0] != 0.0:
+        print(f"{indent}{format_real_coefficient(coefficients[0])} +")
+    last_non_zero_idx = len(coefficients) - 1
+    while coefficients[last_non_zero_idx] == 0.0:
+        last_non_zero_idx = last_non_zero_idx - 1
+    idx = 1
+    parentheses = 0
+    print(indent, end='')
+    while idx < last_non_zero_idx:
+        print(f"{varname} * ", end='')
+        if coefficients[idx] != 0.0:
+            print(f"({format_real_coefficient(coefficients[idx])} +")
+            print(indent, end='')
+            parentheses = parentheses + 1
+        idx = idx + 1
+    print(f"{varname} * {format_real_coefficient(coefficients[last_non_zero_idx])}", end='')
+    for i in range(0, parentheses):
+        print(")", end='')
+    print()
+
+def format_rational_coefficient(coefficient):
+    return str(coefficient).replace("/", ".0 / ") + ".0"
+
+def print_polynomial_with_rational_coefficients(varname, coefficients, indent):
+    if len(coefficients) <= 1:
+        print(f"{indent}{format_rational_coefficient(coefficients[0])}")
+        return
+    if coefficients[0] != 0:
+        print(f"{indent}{format_rational_coefficient(coefficients[0])} +")
+    last_non_zero_idx = len(coefficients) - 1
+    while coefficients[last_non_zero_idx] == 0:
+        last_non_zero_idx = last_non_zero_idx - 1
+    idx = 1
+    parentheses = 0
+    print(indent, end='')
+    while idx < last_non_zero_idx:
+        print(f"{varname} * ", end='')
+        if coefficients[idx] != 0:
+            print(f"({format_rational_coefficient(coefficients[idx])} +")
+            print(indent, end='')
+            parentheses = parentheses + 1
+        idx = idx + 1
+    print(f"{varname} * {format_rational_coefficient(coefficients[last_non_zero_idx])}", end='')
+    for i in range(0, parentheses):
+        print(")", end='')
+    print()
+
+def gamma_at_2_coefficient(k):
+    if k == 0:
+        return 0.0
+    return ((-1)**(k + 1)*(zeta(k + 1) - 1)/(k + 1)).evalf(decimal_precision, maxn=evalf_inner_precision)
+
+def digamma_at_1_coefficient(k):
+    if k == 0:
+        return digamma(1).evalf(decimal_precision, maxn=evalf_inner_precision)
+    return ((-1)**(k + 1) * zeta(k + 1)).evalf(decimal_precision, maxn=evalf_inner_precision)
+
+def digamma_at_2_coefficient(k):
+    if k == 0:
+        return 0.0
+    return ((-1)**(k + 1)*(zeta(k + 1) - 1)).evalf(decimal_precision, maxn=evalf_inner_precision)
+
+def digamma_asymptotic_coefficient(k):
+    if k == 0:
+        return 0.0
+    return bernoulli(2 * k) / (2 * k)
+
+def trigamma_at_1_coefficient(k):
+    return ((-1)**k * (k + 1) * zeta(k + 2)).evalf(decimal_precision, maxn=evalf_inner_precision)
+
+def trigamma_asymptotic_coefficient(k):
+    if k == 0:
+        return 0.0
+    return bernoulli(2 * k)
+
+def tetragamma_at_1_coefficient(k):
+    return ((-1)**(k + 1) * (k + 1) * (k + 2) * zeta(k + 3)).evalf(decimal_precision, maxn=evalf_inner_precision)
+
+def tetragamma_asymptotic_coefficient(k):
+    if k == 0:
+        return 0.0
+    return -(2 * k - 1) * bernoulli(2 * k - 2)
+
+def gammaln_asymptotic_coefficient(k):
+    return bernoulli(2 * k + 2) / (2 * (k + 1) * (2 * k + 1))
+
+def log_1_plus_coefficient(k):
+    if k == S(0):
+        return 0
+    if k % 2 == 0:
+        return S(-1) / k
+    return S(1) / k
+
+def log_1_minus_coefficient(k):
+    if k == 0:
+        return S(0)
+    return S(-1) / k
+
+def x_minus_log_1_plus_coefficient(k):
+    if k <= 1:
+        return S(0)
+    if k % 2 == 0:
+        return S(1) / k
+    return S(-1) / k
+
+def exp_minus_1_coefficient(k):
+    if k == 0:
+        return S(0)
+    return S(1) / factorial(k)
+
+def exp_minus_1_ratio_minus_1_ratio_minus_half_coefficient(k):
+    if k == 0:
+        return S(0)
+    return S(1) / factorial(k + 2)
+
+def get_log_exp_minus_1_ratio_coefficients(count):
+    x = symbols('x')
+    return list(reversed(Poly(log((exp(x) - 1) / x).series(x, 0, count).removeO()).all_coeffs()))
+
+# Formula for mth coefficient of the normcdfln asymptotic:
+# \sum_{n=1}^m (-1)^{n+m+1} / n * \sum_{l1, l2, ..., ln \in N, l1 + l2 + ... + ln = m} (2 * l1 - 1)!! * (2 * l2 - 1)!! * ... * (2 * ln - 1)!!
+# Can be obtained by composing the Taylor expansion for log(1 + x) and asymtotic expansion for erfc
+def normcdfln_asymptotic_coefficient(m):
+    if m == 0:
+        return 0
+    def next(v):
+        idx = 1
+        while idx < len(v) and v[idx] == 0:
+            idx = idx + 1
+        if idx == len(v): return False
+        v0 = v[0]
+        v[0] = 0
+        v[idx] = v[idx] - 1
+        v[idx - 1] = v0 + 1
+        return True
+    result = S((-1)**(2 + m)) * factorial2(2 * m - 1)
+    for n in range(2,m+1):
+        coef = S((-1)**(n + 1 + m)) / n
+        deltas = []
+        for k in range(0, n):
+            deltas.append(0)
+        deltas[-1] = m - n
+        accSum = S(0)
+        while True:
+            accProd = S(1)
+            for delta in deltas:
+                accProd = accProd * factorial2(2 * (delta + 1) - 1)
+            accSum = accSum + accProd
+            if not next(deltas):
+                break
+        result = result + coef * accSum
+    return result
+
+def get_one_minus_sqrt_one_minus_coefficients(count):
+    x = symbols('x')
+    return list(reversed(Poly((1 - sqrt(1 - x)).series(x, 0, count).removeO()).all_coeffs()))
+
+def get_reciprocal_factorial_minus_1_coefficients(count):
+    x = symbols('x')
+    return list(reversed(Poly((1 / gamma(x + 1) - 1).series(x, 0, count).removeO().evalf(decimal_precision, maxn=evalf_inner_precision), x).all_coeffs()))
+
+def main():
+    print_heading_comment(gamma_at_2_indent, "1: Gamma at 2")
+    gamma_at_2_coefficients = [gamma_at_2_coefficient(k) for k in range(0, gamma_at_2_series_length)]
+    print_polynomial_with_real_coefficients(gamma_at_2_variable_name, gamma_at_2_coefficients, gamma_at_2_indent)
+
+    print_heading_comment(digamma_at_1_indent, "2: Digamma at 1")
+    digamma_at_1_coefficients = [digamma_at_1_coefficient(k) for k in range(0, digamma_at_1_series_length)]
+    print_polynomial_with_real_coefficients(digamma_at_1_variable_name, digamma_at_1_coefficients, digamma_at_1_indent)
+
+    print_heading_comment(digamma_at_2_indent, "3: Digamma at 2")
+    digamma_at_2_coefficients = [digamma_at_2_coefficient(k) for k in range(0, digamma_at_2_series_length)]
+    print_polynomial_with_real_coefficients(digamma_at_2_variable_name, digamma_at_2_coefficients, digamma_at_2_indent)
+
+    print_heading_comment(digamma_asymptotic_indent, "4: Digamma asymptotic")
+    digamma_asymptotic_coefficients = [digamma_asymptotic_coefficient(k) for k in range(0, digamma_asymptotic_series_length)]
+    print_polynomial_with_rational_coefficients(digamma_asymptotic_variable_name, digamma_asymptotic_coefficients, digamma_asymptotic_indent)
+
+    print_heading_comment(trigamma_at_1_indent, "5: Trigamma at 1")
+    trigamma_at_1_coefficients = [trigamma_at_1_coefficient(k) for k in range(0, trigamma_at_1_series_length)]
+    print_polynomial_with_real_coefficients(trigamma_at_1_variable_name, trigamma_at_1_coefficients, trigamma_at_1_indent)
+
+    print_heading_comment(trigamma_asymptotic_indent, "6: Trigamma asymptotic")
+    trigamma_asymptotic_coefficients = [trigamma_asymptotic_coefficient(k) for k in range(0, trigamma_asymptotic_series_length)]
+    print_polynomial_with_rational_coefficients(trigamma_asymptotic_variable_name, trigamma_asymptotic_coefficients, trigamma_asymptotic_indent)
+
+    print_heading_comment(tetragamma_at_1_indent, "7: Tetragamma at 1")
+    tetragamma_at_1_coefficients = [tetragamma_at_1_coefficient(k) for k in range(0, tetragamma_at_1_series_length)]
+    print_polynomial_with_real_coefficients(tetragamma_at_1_variable_name, tetragamma_at_1_coefficients, tetragamma_at_1_indent)
+
+    print_heading_comment(tetragamma_asymptotic_indent, "8: Tetragamma asymptotic")
+    tetragamma_asymptotic_coefficients = [tetragamma_asymptotic_coefficient(k) for k in range(0, tetragamma_asymptotic_series_length)]
+    print_polynomial_with_rational_coefficients(tetragamma_asymptotic_variable_name, tetragamma_asymptotic_coefficients, tetragamma_asymptotic_indent)
+
+    print_heading_comment(gammaln_asymptotic_indent, "9: GammaLn asymptotic")
+    gammaln_asymptotic_coefficients = [gammaln_asymptotic_coefficient(k) for k in range(0, gammaln_asymptotic_series_length)]
+    print_polynomial_with_rational_coefficients(gammaln_asymptotic_variable_name, gammaln_asymptotic_coefficients, gammaln_asymptotic_indent)
+
+    print_heading_comment(log_1_plus_indent, "10: log(1 + x)")
+    log_1_plus_coefficients = [log_1_plus_coefficient(k) for k in range(0, log_1_plus_series_length)]
+    print_polynomial_with_rational_coefficients(log_1_plus_variable_name, log_1_plus_coefficients, log_1_plus_indent)
+
+    print_heading_comment(log_1_minus_indent, "11: log(1 - x)")
+    log_1_minus_coefficients = [log_1_minus_coefficient(k) for k in range(0, log_1_minus_series_length)]
+    print_polynomial_with_rational_coefficients(log_1_minus_variable_name, log_1_minus_coefficients, log_1_minus_indent)
+
+    print_heading_comment(x_minus_log_1_plus_indent, "12: x - log(1 + x)")
+    x_minus_log_1_plus_coefficients = [x_minus_log_1_plus_coefficient(k) for k in range(0, x_minus_log_1_plus_series_length)]
+    print_polynomial_with_rational_coefficients(x_minus_log_1_plus_variable_name, x_minus_log_1_plus_coefficients, x_minus_log_1_plus_indent)
+
+    print_heading_comment(exp_minus_1_indent, "13: exp(x) - 1")
+    exp_minus_1_coefficients = [exp_minus_1_coefficient(k) for k in range(0, exp_minus_1_series_length)]
+    print_polynomial_with_rational_coefficients(exp_minus_1_variable_name, exp_minus_1_coefficients, exp_minus_1_indent)
+
+    print_heading_comment(exp_minus_1_ratio_minus_1_ratio_minus_half_indent, "14: ((exp(x) - 1) / x - 1) / x - 0.5")
+    exp_minus_1_ratio_minus_1_ratio_minus_half_coefficients = [exp_minus_1_ratio_minus_1_ratio_minus_half_coefficient(k) for k in range(0, exp_minus_1_ratio_minus_1_ratio_minus_half_series_length)]
+    print_polynomial_with_rational_coefficients(exp_minus_1_ratio_minus_1_ratio_minus_half_variable_name, exp_minus_1_ratio_minus_1_ratio_minus_half_coefficients, exp_minus_1_ratio_minus_1_ratio_minus_half_indent)
+
+    print_heading_comment(log_exp_minus_1_ratio_indent, "15: log(exp(x) - 1) / x")
+    log_exp_minus_1_ratio_coefficients = get_log_exp_minus_1_ratio_coefficients(log_exp_minus_1_ratio_series_length)
+    print_polynomial_with_rational_coefficients(log_exp_minus_1_ratio_variable_name, log_exp_minus_1_ratio_coefficients, log_exp_minus_1_ratio_indent)
+
+    print_heading_comment(normcdfln_asymptotic_indent, "16: normcdfln asymptotic")
+    normcdfln_asymptotic_coefficients = [normcdfln_asymptotic_coefficient(k) for k in range(0, normcdfln_asymptotic_series_length)]
+    print_polynomial_with_rational_coefficients(normcdfln_asymptotic_variable_name, normcdfln_asymptotic_coefficients, normcdfln_asymptotic_indent)
+    
+    print_heading_comment(one_minus_sqrt_one_minus_indent, "17: 1 - sqrt(1 - x)")
+    one_minus_sqrt_one_minus_coefficients = get_one_minus_sqrt_one_minus_coefficients(one_minus_sqrt_one_minus_series_length)
+    print_polynomial_with_rational_coefficients(one_minus_sqrt_one_minus_variable_name, one_minus_sqrt_one_minus_coefficients, one_minus_sqrt_one_minus_indent)
+    
+    print_heading_comment(reciprocal_factorial_minus_1_indent, "18: Reciprocal factorial minus 1")
+    reciprocal_factorial_minus_1_coefficients = get_reciprocal_factorial_minus_1_coefficients(reciprocal_factorial_minus_1_series_length)
+    print_polynomial_with_real_coefficients(reciprocal_factorial_minus_1_variable_name, reciprocal_factorial_minus_1_coefficients, reciprocal_factorial_minus_1_indent)
+
+if __name__ == '__main__': main()

src/Tools/GenerateSeries/Tools.GenerateSeries.pyproj

@@ -0,0 +1,36 @@
+<Project DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003" ToolsVersion="4.0">
+  <PropertyGroup>
+    <Configuration Condition=" '$(Configuration)' == '' ">Debug</Configuration>
+    <SchemaVersion>2.0</SchemaVersion>
+    <ProjectGuid>d7562cc3-d48a-481d-b40c-07ee51ebbf50</ProjectGuid>
+    <ProjectHome>
+    </ProjectHome>
+    <StartupFile>GenerateSeries.py</StartupFile>
+    <SearchPath>
+    </SearchPath>
+    <WorkingDirectory>.</WorkingDirectory>
+    <OutputPath>.</OutputPath>
+    <Name>Tools.GenerateSeries</Name>
+    <RootNamespace>GenerateSeries</RootNamespace>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)' == 'Debug' ">
+    <DebugSymbols>true</DebugSymbols>
+    <EnableUnmanagedDebugging>false</EnableUnmanagedDebugging>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)' == 'Release' ">
+    <DebugSymbols>true</DebugSymbols>
+    <EnableUnmanagedDebugging>false</EnableUnmanagedDebugging>
+  </PropertyGroup>
+  <ItemGroup>
+    <Compile Include="GenerateSeries.py" />
+  </ItemGroup>
+  <Import Project="$(MSBuildExtensionsPath32)\Microsoft\VisualStudio\v$(VisualStudioVersion)\Python Tools\Microsoft.PythonTools.targets" />
+  <!-- Uncomment the CoreCompile target to enable the Build command in
+       Visual Studio and specify your pre- and post-build commands in
+       the BeforeBuild and AfterBuild targets below. -->
+  <!--<Target Name="CoreCompile" />-->
+  <Target Name="BeforeBuild">
+  </Target>
+  <Target Name="AfterBuild">
+  </Target>
+</Project>

src/Tutorials/Tutorials.csproj

@@ -46,6 +46,12 @@
     <DefineConstants>$(DefineConstants);DEBUG</DefineConstants>
   </PropertyGroup>
 
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|AnyCPU' OR '$(Configuration)|$(Platform)'=='ReleaseFull|AnyCPU' OR '$(Configuration)|$(Platform)'=='ReleaseCore|AnyCPU'">
+    <DebugSymbols>true</DebugSymbols>
+    <DebugType>pdbonly</DebugType>
+    <Optimize>true</Optimize>
+  </PropertyGroup>
+
   <ItemGroup>
     <ProjectReference Include="..\Compiler\Compiler.csproj" />
     <ProjectReference Include="..\Runtime\Runtime.csproj" />

test/Learners/LearnersTests/LearnersTests.csproj

@@ -41,6 +41,7 @@
 
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|AnyCPU' OR '$(Configuration)|$(Platform)'=='ReleaseFull|AnyCPU' OR '$(Configuration)|$(Platform)'=='ReleaseCore|AnyCPU'">
     <DebugType>pdbonly</DebugType>
+    <DebugSymbols>true</DebugSymbols>
     <Optimize>true</Optimize>
   </PropertyGroup>
 

test/Learners/TestApp/TestApp.csproj

@@ -44,6 +44,7 @@
   </PropertyGroup>
 
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|AnyCPU' OR '$(Configuration)|$(Platform)'=='ReleaseFull|AnyCPU' OR '$(Configuration)|$(Platform)'=='ReleaseCore|AnyCPU'">
+    <DebugSymbols>true</DebugSymbols>
     <DebugType>pdbonly</DebugType>
     <Optimize>true</Optimize>
   </PropertyGroup>

test/TestPublic/TestPublic.csproj

@@ -41,6 +41,12 @@
     <DefineConstants>$(DefineConstants);DEBUG</DefineConstants>
   </PropertyGroup>
 
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|AnyCPU' OR '$(Configuration)|$(Platform)'=='ReleaseFull|AnyCPU' OR '$(Configuration)|$(Platform)'=='ReleaseCore|AnyCPU'">
+    <DebugType>pdbonly</DebugType>
+    <DebugSymbols>true</DebugSymbols>
+    <Optimize>true</Optimize>
+  </PropertyGroup>
+
   <ItemGroup>
     <Compile Include="..\..\src\Shared\SharedAssemblyFileVersion.cs" />
     <Compile Include="..\..\src\Shared\SharedAssemblyInfo.cs" />

test/Tests/Tests.csproj

@@ -42,6 +42,12 @@
     <DefineConstants>$(DefineConstants);DEBUG</DefineConstants>
   </PropertyGroup>
 
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|AnyCPU' OR '$(Configuration)|$(Platform)'=='ReleaseFull|AnyCPU' OR '$(Configuration)|$(Platform)'=='ReleaseCore|AnyCPU'">
+    <DebugType>pdbonly</DebugType>
+    <DebugSymbols>true</DebugSymbols>
+    <Optimize>true</Optimize>
+  </PropertyGroup>
+
   <ItemGroup>
     <PackageReference Include="Microsoft.NET.Test.Sdk" Version="15.8.0" />
     <PackageReference Include="Newtonsoft.Json" Version="11.0.2" />
@@ -89,4 +95,4 @@
       <Visible>false</Visible>
     </EmbeddedResource>
   </ItemGroup>
-</Project>
+</Project>