Special functions improvements (#253)

- Ensured argument consistency in NormalCdfIntegralTest
- Fixed some test values
- Exp-Sinh quadrature in LogisticGaussian
- Arithmetic-precision-based constants in logistic gaussian
- Named const for Ulp(1)
- Generating series for gamma(x) - 1/x
- Fixed integer overflow in BGRat
- GenerateSeries can print arrays of bigfloats.

src/Runtime/Core/Maths/Quadrature.cs

@@ -15,6 +15,8 @@ namespace Microsoft.ML.Probabilistic.Math
     /// </summary>
     public static class Quadrature
     {
+        private const int AdaptiveQuadratureMaxNodes = 10000;
+
         /// <summary>
         /// Integrate the function f from -Infinity to Infinity
         /// </summary>
@@ -35,7 +37,103 @@ namespace Microsoft.ML.Probabilistic.Math
                 double sinX = System.Math.Sin(x);
                 return f(scale / System.Math.Tan(x)) / (sinX * sinX);
             }
-            return scale * AdaptiveTrapeziumRule(fInvTan, nodeCount, 0, System.Math.PI, relTol, 10000);
+            return scale * AdaptiveTrapeziumRule(fInvTan, nodeCount, 0, System.Math.PI, relTol, AdaptiveQuadratureMaxNodes);
+        }
+
+        /// <summary>
+        /// Integrate the function f from 0 to +Infinity using exp-sinh quadrature.
+        /// </summary>
+        /// <param name="f">The function to integrate</param>
+        /// <param name="scale">A positive tuning parameter.
+        /// <paramref name="f"/> is assumed to be continuous and at least one of <paramref name="f"/>(0) and
+        /// cosh(arsinh(2 * ln(scale)/pi)) * scale * <paramref name="f"/>(<paramref name="scale"/>) to be non-zero.</param>
+        /// <param name="relTol">A threshold to stop subdividing</param>
+        /// <returns></returns>
+        public static double AdaptiveExpSinh(Converter<double, double> f, double scale, double relTol)
+        {
+            const double halfPi = System.Math.PI / 2;
+            const double ln2ByHalfPi = MMath.Ln2 / halfPi;
+            // Let g(x) = f(exp(pi/2 * sinh(x))) * cosh(x) * exp(pi/2 * sinh(x)),
+            // fExpSinh(x) = g(x) + g(-x), when x != 0, fExpSinh(0) = g(0),
+            // then int_-inf^+inf g(x) dx = int_0^+inf fExpSinh(x) dx
+            double fExpSinh(double x)
+            {
+                if (x == 0)
+                    return f(1);
+                double halfPiSinhX = halfPi * System.Math.Sinh(x);
+                double expHalfPiSinhX = System.Math.Exp(halfPiSinhX);
+                double expHalfPiSinhNegX = System.Math.Exp(-halfPiSinhX);
+                return System.Math.Cosh(x) * (expHalfPiSinhX * f(expHalfPiSinhX) + expHalfPiSinhNegX * f(expHalfPiSinhNegX));
+            }
+            double invHalfPiLnScale = System.Math.Log(scale) / halfPi;
+            double rescale = System.Math.Abs(System.Math.Log(invHalfPiLnScale + System.Math.Sqrt(invHalfPiLnScale * invHalfPiLnScale + 1))); // abs . arsinh
+            while (fExpSinh(rescale) == 0)
+            {
+                invHalfPiLnScale -= ln2ByHalfPi;
+                rescale = System.Math.Abs(System.Math.Log(invHalfPiLnScale + System.Math.Sqrt(invHalfPiLnScale * invHalfPiLnScale + 1)));
+            }
+            return halfPi * AdaptivePositiveHalfAxisTrapeziumRule(fExpSinh, rescale, relTol, AdaptiveQuadratureMaxNodes / 2 + 1);
+        }
+
+        /// <summary>
+        /// Integrate the function f from 0 to +Infinity
+        /// </summary>
+        /// <param name="f">The function to integrate. Must have at most one extremum.</param>
+        /// <param name="scale">A positive tuning parameter.
+        /// <paramref name="f"/> is assumed to be continuous on (0, + Infinity) and non-negligible somewhere on (0, <paramref name="scale"/>]</param>
+        /// <param name="relTol">A threshold to stop subdividing</param>
+        /// <param name="maxNodes">Another threshold to stop subdividing</param>
+        /// <returns></returns>
+        public static double AdaptivePositiveHalfAxisTrapeziumRule(Converter<double, double> f, double scale, double relTol, int maxNodes)
+        {
+            double intervalWidth = 1.0 / 8;
+            double sumf2 = 0;
+            double sumf1 = f(0);
+            double x = 0;
+            double oldSum;
+            int usedNodes = 1;
+            do
+            {
+                ++usedNodes;
+                x += intervalWidth;
+                oldSum = sumf1;
+                sumf1 += f(x);
+            }
+            while (!MMath.AreEqual(oldSum, sumf1) || x < scale);
+
+            if (x > scale)
+                scale = x;
+
+            while (usedNodes < maxNodes)
+            {
+                intervalWidth /= 2;
+                sumf2 = sumf1;
+                x = intervalWidth;
+                sumf2 += f(x);
+                do
+                {
+                    ++usedNodes;
+                    if (x < scale)
+                        x += 2 * intervalWidth;
+                    else
+                        x += intervalWidth;
+                    oldSum = sumf2;
+                    sumf2 += f(x);
+                }
+                while (!MMath.AreEqual(oldSum, sumf2) || x < scale);
+                if (x > scale)
+                    scale = x;
+                double i1 = sumf1 * intervalWidth * 2;
+                double i2 = sumf2 * intervalWidth;
+                double err_est = System.Math.Abs((i1 - i2) / i2);
+                if (err_est < relTol)
+                {
+                    return i2;
+                }
+                sumf1 = sumf2;
+            }
+
+            return sumf1 * intervalWidth;
         }
 
         /// <summary>

src/Runtime/Core/Maths/SpecialFunctions.cs

@@ -156,7 +156,7 @@ namespace Microsoft.ML.Probabilistic.Math
         private static double BGRat(double x, double a, double b, double epsilon)
         {
             List<double> p = new List<double>();
-            List<int> twoNPlus1Factorial = new List<int>();
+            List<double> twoNPlus1Factorial = new List<double>();
             p.Add(1.0);
             twoNPlus1Factorial.Add(1);
             double T = a + 0.5 * (b - 1);
@@ -180,6 +180,8 @@ namespace Microsoft.ML.Probabilistic.Math
                 bPlus2n += 2;
                 bPlus2nPlus1 += 2;
                 lnXTerm *= lnXOver2Sq;
+                // This can overflow, which will result in infinity. This is OK, because
+                // it will mean only that some of the terms added to currP will be zero.
                 twoNPlus1Factorial.Add(twoNPlus1Factorial.Last() * (2 * n) * (2 * n + 1));
 
                 double currP = 0.0;
@@ -499,6 +501,8 @@ namespace Microsoft.ML.Probabilistic.Math
             return mult * result;
         }
 
+        private const double DefaultBetaEpsilon = 1e-15;
+
         /// <summary>
         /// Computes the regularized incomplete beta function: int_0^x t^(a-1) (1-t)^(b-1) dt / Beta(a,b)
         /// </summary>
@@ -508,7 +512,7 @@ namespace Microsoft.ML.Probabilistic.Math
         /// <param name="epsilon">A tolerance for terminating the series calculation.</param>
         /// <returns>The incomplete beta function at (<paramref name="x"/>, <paramref name="a"/>, <paramref name="b"/>).</returns>
         /// <remarks>The beta function is obtained by setting x to 1.</remarks>
-        public static double Beta(double x, double a, double b, double epsilon = 1e-15)
+        public static double Beta(double x, double a, double b, double epsilon = DefaultBetaEpsilon)
         {
             double result = 0.0;
             bool swap = false;
@@ -1168,6 +1172,8 @@ namespace Microsoft.ML.Probabilistic.Math
             return result;
         }
 
+        private static readonly double Ulp1 = Ulp(1.0);
+
         /// <summary>
         /// Compute the regularized upper incomplete Gamma function: int_x^inf t^(a-1) exp(-t) dt / Gamma(a)
         /// </summary>
@@ -1209,7 +1215,7 @@ namespace Microsoft.ML.Probabilistic.Math
                 return GammaAsympt(a, x, true);
             else if (x > 1.5)
                 return GammaUpperConFrac(a, x);
-            else if (a <= 1e-16)
+            else if (a <= Ulp1)
                 // Gamma(a) = 1/a for a <= 1e-16
                 return a * GammaUpperSeries(a, x, false);
             else
@@ -1430,7 +1436,7 @@ namespace Microsoft.ML.Probabilistic.Math
             }
             else
             {
-                if (Math.Abs(alogx) <= 1e-16) offset = GammaSeries(a) - logx;
+                if (Math.Abs(alogx) <= Ulp1) offset = GammaSeries(a) - logx;
                 else offset = GammaSeries(a) - xaMinus1 / a;
                 scale = 1 + xaMinus1;
                 term = x;
@@ -1465,14 +1471,17 @@ var('x');
 f = gamma(x)-1/x
 [c[0].n(100) for c in f.taylor(x,0,16).coefficients()]
                  */
-                return Digamma1
-                    + x * (0.98905599532797255539539565150
-                    + x * (-0.90747907608088628901656016736
-                    + x * (0.98172808683440018733638029402
-                    + x * (-0.98199506890314520210470141379
-                    + x * (0.99314911462127619315386725333
-                    + x * (-0.99600176044243153397007841966
-                    ))))));
+                return
+                    // Truncated series 19: Gamma(x) - 1/x
+                    // Generated automatically by /src/Tools/PythonScripts/GenerateSeries.py
+                    -0.577215664901532860606512090082 +
+                    x * (0.989055995327972555395395651501 +
+                    x * (-0.907479076080886289016560167356 +
+                    x * (0.981728086834400187336380294022 +
+                    x * (-0.981995068903145202104701413791 +
+                    x * (0.993149114621276193153867253329 +
+                    x * -0.996001760442431533970078419665)))))
+                    ;
         }
 
         /// <summary>
@@ -2181,7 +2190,7 @@ f = gamma(x)-1/x
             const double smallestNormalized = 1e-308;
             const double smallestNormalizedOverEpsilon = smallestNormalized / double.Epsilon;
             // the number of iterations required may grow with n, so we need to explicitly test for convergence.
-            for (int i = 0; i < 1000; i++)
+            for (int i = 0; i < 10000; i++)
             {
                 double numerNew = numer + a * numerPrev;
                 double denomNew = denom + a * denomPrev;
@@ -2607,7 +2616,8 @@ f = gamma(x)-1/x
             // Second term of the Taylor series
             sum += Qderiv * rPowerN;
             double sumOld = sum;
-            for (int n = 2; n <= 100; n++)
+            int maxIterations = 1000;
+            for (int n = 2; n <= maxIterations; n++)
             {
                 //Console.WriteLine($"n = {n - 1} sum = {sum:g17}");
                 double dlogphiOverFactorial;
@@ -2623,7 +2633,7 @@ f = gamma(x)-1/x
                 Qderivs.Add(QderivOverFactorial);
                 rPowerN *= r;
                 sum += QderivOverFactorial * rPowerN;
-                if ((sum > double.MaxValue) || double.IsNaN(sum) || n >= 100)
+                if ((sum > double.MaxValue) || double.IsNaN(sum) || n >= maxIterations)
                     throw new Exception($"NormalCdfRatioTaylor not converging for x={x:g17}, y={y:g17}, r={r:g17}");
                 if (AreEqual(sum, sumOld)) break;
                 sumOld = sum;
@@ -2708,6 +2718,7 @@ rr = mpf('-0.99999824265582826');
         }
 
         internal static bool TraceConFrac;
+        private static readonly double NormalCdfRatioConfracTolerance = Ulp1 * 1024;
 
         /// <summary>
         /// Returns NormalCdf divided by N(x;0,1) N((y-rx)/sqrt(1-r^2);0,1), multiplied by scale.
@@ -2921,7 +2932,7 @@ rr = mpf('-0.99999824265582826');
                 b /= x;
             }
             double bIncr = sqrtomr2 / x;
-            const int iterationCount = 1000;
+            const int iterationCount = 10000;
             for (int i = 1; i <= iterationCount; i++)
             {
                 double numerNew, denomNew;
@@ -3013,7 +3024,7 @@ rr = mpf('-0.99999824265582826');
                         Trace.WriteLine($"iter {i}: result={result:g17} c={c:g17} cOdd={cOdd:g17} numer={numer:g17} numer2={numer2:g17} denom={denom:g17} numerPrev={numerPrev:g17}");
                     if ((result > double.MaxValue) || double.IsNaN(result) || result < 0 || i >= iterationCount - 1)
                         throw new Exception($"NormalCdfRatioConFrac2 not converging for x={x:g17} y={y:g17} r={r:g17} sqrtomr2={sqrtomr2:g17} scale={scale:g17}");
-                    if (AreEqual(result, resultPrev) || AbsDiff(result, resultPrev, 0) < 1e-13)
+                    if (AreEqual(result, resultPrev) || AbsDiff(result, resultPrev, 0) < NormalCdfRatioConfracTolerance)
                         break;
                     resultPrev = result;
                 }
@@ -3514,7 +3525,7 @@ rr = mpf('-0.99999824265582826');
             //   log(1+exp(x))-x = log(1+exp(-x)) <= exp(-x)
             // Thus we should use the approximation when exp(-x) <= ulp(x)/2
             // ulp(1) < ulp(x)  therefore x > -log(ulp(1)/2)
-            if (x > -logEpsilon)
+            if (x > -logHalfUlpPrev1)
                 return x;
             else
                 return Log1Plus(Math.Exp(x));
@@ -3965,10 +3976,13 @@ rr = mpf('-0.99999824265582826');
             vf = Ex2 - mf * mf;
         }
 
-        // Math.Exp(-745.14) == 0
-        private const double log0 = -745.14;
-        // 1-Math.Exp(-38) == 1
-        private const double logEpsilon = -38;
+        // Math.Exp(log0) == 0
+        private static readonly double log0 = Math.Log(double.Epsilon) - Ln2;
+        // 1-Math.Exp(logHalfUlpPrev1) == 1
+        private static readonly double logHalfUlpPrev1 = Math.Log((1.0 - PreviousDouble(1.0)) / 2);
+        private static readonly double logisticGaussianQuadratureRelativeTolerance = 512 * Ulp1;
+        private static readonly double logisticGaussianDerivativeQuadratureRelativeTolerance = 1024 * 512 * Ulp1;
+        private static readonly double logisticGaussianSeriesApproximmationThreshold = 1e-8;
 
         /// <summary>
         /// Calculate sigma(m,v) = \int N(x;m,v) logistic(x) dx
@@ -3990,19 +4004,19 @@ rr = mpf('-0.99999824265582826');
             // use the upper bound exp(m+v/2) to prune cases that must be zero or one
             if (mean + halfVariance < log0)
                 return 0.0;
-            if (-mean + halfVariance < logEpsilon)
+            if (-mean + halfVariance < logHalfUlpPrev1)
                 return 1.0;
 
             // use the upper bound 0.5 exp(-0.5 m^2/v) to prune cases that must be zero or one
             double q = -0.5 * mean * mean / variance - MMath.Ln2;
             if (mean <= 0 && mean + variance >= 0 && q < log0)
                 return 0.0;
-            if (mean >= 0 && variance - mean >= 0 && q < logEpsilon)
+            if (mean >= 0 && variance - mean >= 0 && q < logHalfUlpPrev1)
                 return 1.0;
             // sigma(|m|,v) <= 0.5 + |m| sigma'(0,v)
             // sigma'(0,v) <= N(0;0,v+8/pi)
             //double d0Upper = MMath.InvSqrt2PI / Math.Sqrt(variance + 8 / Math.PI);
-            if (mean * mean / (variance + 8 / Math.PI) < 1e-8)
+            if (mean * mean / (variance + 8 / Math.PI) < logisticGaussianSeriesApproximmationThreshold)
             {
                 double deriv = LogisticGaussianDerivative(0, variance);
                 return 0.5 + mean * deriv;
@@ -4010,14 +4024,14 @@ rr = mpf('-0.99999824265582826');
 
             // Handle tail cases using the following exact formulas:
             // sigma(m,v) = 1 - exp(-m+v/2) + exp(-2m+2v) - exp(-3m+9v/2) sigma(m-3v,v)
-            if (2 * (variance - mean) < logEpsilon)
+            if (2 * (variance - mean) < logHalfUlpPrev1)
                 return 1.0 - Math.Exp(halfVariance - mean);
-            if (-3 * mean + 9 * halfVariance < logEpsilon)
+            if (-3 * mean + 9 * halfVariance < logHalfUlpPrev1)
                 return 1.0 - Math.Exp(halfVariance - mean) + Math.Exp(2 * (variance - mean));
             // sigma(m,v) = exp(m+v/2) - exp(2m+2v) + exp(3m + 9v/2) (1 - sigma(m+3v,v))
-            if (mean + 1.5 * variance < logEpsilon)
+            if (mean + 1.5 * variance < logHalfUlpPrev1)
                 return Math.Exp(mean + halfVariance);
-            if (2 * mean + 4 * variance < logEpsilon)
+            if (2 * mean + 4 * variance < logHalfUlpPrev1)
                 return Math.Exp(mean + halfVariance) * (1 - Math.Exp(mean + 1.5 * variance));
 
             if (variance > LogisticGaussianVarianceThreshold)
@@ -4032,7 +4046,8 @@ rr = mpf('-0.99999824265582826');
                 }
                 double upperBound = mean + Math.Sqrt(variance);
                 double scale = Math.Max(upperBound, 10) / sqrtv;
-                return new ExtendedDouble(Quadrature.AdaptiveClenshawCurtis(f, scale, 32, 1e-13), shift).ToDouble();
+                return new ExtendedDouble(Quadrature.AdaptiveExpSinh(f, scale, logisticGaussianQuadratureRelativeTolerance / 2) +
+                    Quadrature.AdaptiveExpSinh(x => f(-x), scale, logisticGaussianQuadratureRelativeTolerance / 2), shift).ToDouble();
             }
             else
             {
@@ -4078,9 +4093,9 @@ rr = mpf('-0.99999824265582826');
 
             // Handle the tail cases using the following exact formula:
             // sigma'(m,v) = exp(-m+v/2) -2 exp(-2m+2v) +3 exp(-3m+9v/2) sigma(m-3v,v) - exp(-3m+9v/2) sigma'(m-3v,v)
-            if (-mean + 1.5 * variance < logEpsilon)
+            if (-mean + 1.5 * variance < logHalfUlpPrev1)
                 return Math.Exp(halfVariance - mean);
-            if (-2 * mean + 4 * variance < logEpsilon)
+            if (-2 * mean + 4 * variance < logHalfUlpPrev1)
                 return Math.Exp(halfVariance - mean) - 2 * Math.Exp(2 * (variance - mean));
 
             if (variance > LogisticGaussianVarianceThreshold)
@@ -4090,7 +4105,7 @@ rr = mpf('-0.99999824265582826');
                 {
                     return Math.Exp(MMath.LogisticLn(x) + MMath.LogisticLn(-x) + Gaussian.GetLogProb(x, mean, variance) - shift);
                 }
-                return new ExtendedDouble(Quadrature.AdaptiveClenshawCurtis(f, 10, 32, 1e-10), shift).ToDouble();
+                return new ExtendedDouble(Quadrature.AdaptiveClenshawCurtis(f, 10, 32, logisticGaussianDerivativeQuadratureRelativeTolerance), shift).ToDouble();
             }
             else
             {
@@ -4139,14 +4154,14 @@ rr = mpf('-0.99999824265582826');
 
             // Handle the tail cases using the following exact formulas:
             // sigma''(m,v) = -exp(-m+v/2) +4 exp(-2m+2v) -9 exp(-3m+9v/2) sigma(m-3v,v) +6 exp(-3m+9v/2) sigma'(m-3v,v) - exp(-3m+9v/2) sigma''(m-3v,v)
-            if (-mean + 1.5 * variance < logEpsilon)
+            if (-mean + 1.5 * variance < logHalfUlpPrev1)
                 return -Math.Exp(halfVariance - mean);
-            if (-2 * mean + 4 * variance < logEpsilon)
+            if (-2 * mean + 4 * variance < logHalfUlpPrev1)
                 return -Math.Exp(halfVariance - mean) + 4 * Math.Exp(2 * (variance - mean));
             // sigma''(m,v) = exp(m+v/2) -4 exp(2m+2v) +9 exp(3m + 9v/2) (1 - sigma(m+3v,v)) - 6 exp(3m+9v/2) sigma'(m+3v,v) - exp(3m + 9v/2) sigma''(m+3v,v)
-            if (mean + 1.5 * variance < logEpsilon)
+            if (mean + 1.5 * variance < logHalfUlpPrev1)
                 return Math.Exp(mean + halfVariance);
-            if (2 * mean + 4 * variance < logEpsilon)
+            if (2 * mean + 4 * variance < logHalfUlpPrev1)
                 return Math.Exp(mean + halfVariance) * (1 - 4 * Math.Exp(mean + 1.5 * variance));
 
             if (variance > LogisticGaussianVarianceThreshold)
@@ -4159,7 +4174,7 @@ rr = mpf('-0.99999824265582826');
                     double OneMinus2Sigma = -Math.Tanh(x / 2);
                     return OneMinus2Sigma * Math.Exp(logSigma + log1MinusSigma + Gaussian.GetLogProb(x, mean, variance) - shift);
                 }
-                return new ExtendedDouble(Quadrature.AdaptiveClenshawCurtis(f, 10, 32, 1e-10), shift).ToDouble();
+                return new ExtendedDouble(Quadrature.AdaptiveClenshawCurtis(f, 10, 32, logisticGaussianDerivativeQuadratureRelativeTolerance), shift).ToDouble();
             }
             else
             {

src/Tools/PythonScripts/ComputeSpecialFunctionsTestValues.py

@@ -65,7 +65,7 @@ def normal_cdf2(x, y, r):
         y = z
 
     # Avoid quadrature with r < 0 since it is sometimes inaccurate.
-    if r < 0 and x - y <= 0:
+    if r < 0 and x - y < 0:
         # phi(x,y,r) = phi(inf,y,r) - phi(-x,y,-r)
         # phi(x,y,r) = phi(x,inf,r) - phi(x,-y,-r)
         return ncdf(x) - normal_cdf2(x, -y, -r)
@@ -196,7 +196,7 @@ def logistic_gaussian_deriv(m, v):
     return result
 
 def logistic_gaussian_deriv2(m, v):
-    if m == inf or m == -inf or v == inf:
+    if m == inf or m == -inf or v == inf or m == mpf(0):
         return mpf(0)
     # The integration routine below is obtained by substituting x = atanh(t)
     # into the definition of logistic_gaussian''

src/Tools/PythonScripts/GenerateSeries.py

@@ -5,8 +5,8 @@ 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
+default_decimal_precision = 30
+default_evalf_inner_precision = 500
 
 gamma_at_2_series_length = 26
 gamma_at_2_variable_name = "dx"
@@ -80,12 +80,16 @@ reciprocal_factorial_minus_1_series_length = 17
 reciprocal_factorial_minus_1_variable_name = "x"
 reciprocal_factorial_minus_1_indent = "                "
 
+gamma_minus_reciprocal_series_length = 7
+gamma_minus_reciprocal_variable_name = "x"
+gamma_minus_reciprocal_indent = "                    "
+
 def print_heading_comment(indent, header):
     print(f"{indent}// Truncated series {header}")
     print(f"{indent}// Generated automatically by /src/Tools/PythonScripts/GenerateSeries.py")
 
-def format_real_coefficient(coefficient):
-    return str(coefficient)
+def format_real_coefficient(coefficient, decimal_precision = default_decimal_precision, evalf_inner_precision = default_evalf_inner_precision):
+    return str(coefficient.evalf(decimal_precision, maxn=evalf_inner_precision))
 
 def print_polynomial_with_real_coefficients(varname, coefficients, indent):
     if len(coefficients) <= 1:
@@ -111,6 +115,19 @@ def print_polynomial_with_real_coefficients(varname, coefficients, indent):
         print(")", end='')
     print()
 
+def print_big_float_array(coefficients, decimal_precision, evalf_inner_precision):
+    print("new BigFloat[]")
+    print("{")
+    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 = 0
+    while idx < last_non_zero_idx:
+        print(f'    BigFloatFactory.Create("{format_real_coefficient(coefficients[idx], decimal_precision, evalf_inner_precision)}"),')
+        idx = idx + 1
+    print(f'    BigFloatFactory.Create("{format_real_coefficient(coefficients[last_non_zero_idx], decimal_precision, evalf_inner_precision)}")')
+    print("};")
+
 def format_rational_coefficient(coefficient):
     return str(coefficient).replace("/", ".0 / ") + ".0"
 
@@ -143,38 +160,38 @@ def print_polynomial_with_rational_coefficients(varname, coefficients, indent):
 
 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)
+        return S(0)
+    return ((-1)**(k + 1)*(zeta(k + 1) - 1)/(k + 1))
 
 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)
+        return digamma(1)
+    return ((-1)**(k + 1) * zeta(k + 1))
 
 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)
+        return S(0)
+    return ((-1)**(k + 1)*(zeta(k + 1) - 1))
 
 def digamma_asymptotic_coefficient(k):
     if k == 0:
-        return 0.0
+        return S(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)
+    return ((-1)**k * (k + 1) * zeta(k + 2))
 
 def trigamma_asymptotic_coefficient(k):
     if k == 0:
-        return 0.0
+        return S(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)
+    return ((-1)**(k + 1) * (k + 1) * (k + 2) * zeta(k + 3))
 
 def tetragamma_asymptotic_coefficient(k):
     if k == 0:
-        return 0.0
+        return S(0)
     return -(2 * k - 1) * bernoulli(2 * k - 2)
 
 def gammaln_asymptotic_coefficient(k):
@@ -182,7 +199,7 @@ def gammaln_asymptotic_coefficient(k):
 
 def log_1_plus_coefficient(k):
     if k == S(0):
-        return 0
+        return S(0)
     if k % 2 == 0:
         return S(-1) / k
     return S(1) / k
@@ -218,7 +235,7 @@ def get_log_exp_minus_1_ratio_coefficients(count):
 # 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
+        return S(0)
     def next(v):
         idx = 1
         while idx < len(v) and v[idx] == 0:
@@ -253,7 +270,11 @@ def get_one_minus_sqrt_one_minus_coefficients(count):
 
 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()))
+    return list(reversed(Poly((1 / gamma(x + 1) - 1).series(x, 0, count).removeO(), x).all_coeffs()))
+
+def get_gamma_minus_reciprocal_coefficients(count):
+    x = symbols('x')
+    return list(reversed(Poly((gamma(x) - 1 / x).series(x, 0, count).removeO(), x).all_coeffs()))
 
 def main():
     print_heading_comment(gamma_at_2_indent, "1: Gamma at 2")
@@ -327,5 +348,54 @@ def main():
     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)
+    
+    print_heading_comment(gamma_minus_reciprocal_indent, "19: Gamma(x) - 1/x")
+    gamma_minus_reciprocal_coefficients = get_gamma_minus_reciprocal_coefficients(gamma_minus_reciprocal_series_length)
+    print_polynomial_with_real_coefficients(gamma_minus_reciprocal_variable_name, gamma_minus_reciprocal_coefficients, gamma_minus_reciprocal_indent)
 
-if __name__ == '__main__': main()
+def big_float_main():
+    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, 10)]
+    print_big_float_array(trigamma_at_1_coefficients, 50, 500)
+
+    print_heading_comment(trigamma_asymptotic_indent, "6: Trigamma asymptotic")
+    trigamma_asymptotic_coefficients = [trigamma_asymptotic_coefficient(k) for k in range(0, 32)]
+    print_big_float_array(trigamma_asymptotic_coefficients, 50, 500)
+
+    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, 11)]
+    print_big_float_array(tetragamma_at_1_coefficients, 50, 500)
+
+    print_heading_comment(tetragamma_asymptotic_indent, "8: Tetragamma asymptotic")
+    tetragamma_asymptotic_coefficients = [tetragamma_asymptotic_coefficient(k) for k in range(0, 32)]
+    print_big_float_array(tetragamma_asymptotic_coefficients, 50, 500)
+
+    print_heading_comment(gammaln_asymptotic_indent, "9: GammaLn asymptotic")
+    gammaln_asymptotic_coefficients = [gammaln_asymptotic_coefficient(k) for k in range(0, 31)]
+    print_big_float_array(gammaln_asymptotic_coefficients, 50, 500)
+
+    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, 50)]
+    print_big_float_array(log_1_minus_coefficients, 50, 500)
+
+    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, 26)]
+    print_big_float_array(x_minus_log_1_plus_coefficients, 50, 500)
+
+    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, 19)]
+    print_big_float_array(exp_minus_1_ratio_minus_1_ratio_minus_half_coefficients, 50, 500)
+
+    print_heading_comment(normcdfln_asymptotic_indent, "16: normcdfln asymptotic")
+    normcdfln_asymptotic_coefficients = [normcdfln_asymptotic_coefficient(k) for k in range(0, 19)]
+    print_big_float_array(normcdfln_asymptotic_coefficients, 50, 500)
+    
+    print_heading_comment(reciprocal_factorial_minus_1_indent, "18: Reciprocal factorial minus 1")
+    reciprocal_factorial_minus_1_coefficients = get_reciprocal_factorial_minus_1_coefficients(22)
+    print_big_float_array(reciprocal_factorial_minus_1_coefficients, 50, 500)
+    
+    print_heading_comment(gamma_minus_reciprocal_indent, "19: Gamma(x) - 1/x")
+    gamma_minus_reciprocal_coefficients = get_gamma_minus_reciprocal_coefficients(30)
+    print_big_float_array(gamma_minus_reciprocal_coefficients, 50, 500)
+
+if __name__ == '__main__': big_float_main()

test/Tests/Core/SpecialFunctionsTest.cs

@@ -791,9 +791,9 @@ exp(x*x/4)*pcfu(0.5+n,-x)
         [Fact]
         public void NormalCdfIntegralTest()
         {
-            Assert.True(0 <= MMath.NormalCdfIntegral(1.9018718309533485E+77, -1.9018718309533485E+77, -1, 8.17880416082724E-79).Mantissa);
-            Assert.True(0 <= MMath.NormalCdfIntegral(213393529.2046707, -213393529.2046707, -1, 7.2893668811495072E-10).Mantissa);
-            Assert.True(0 < MMath.NormalCdfIntegral(-0.42146853220760722, 0.42146843802130329, -0.99999999999999989, 6.2292398855983019E-09).Mantissa);
+            Assert.True(0 <= NormalCdfIntegral(190187183095334850882507750944849586799124505055478568794871547478488387682304.0, -190187183095334850882507750944849586799124505055478568794871547478488387682304.0, -1, 0.817880416082724044547388352452631856079457366800004151664125953519049673808376291470533145141236089924006896061006277409614237094627499958581030715374379576478204968748786874650796450332240045653919846557755590765736997127532958984375e-78).Mantissa);
+            Assert.True(0 <= NormalCdfIntegral(213393529.2046706974506378173828125, -213393529.2046706974506378173828125, -1, 0.72893668811495072384656764856902984306419313043079455383121967315673828125e-9).Mantissa);
+            Assert.True(0 < NormalCdfIntegral(-0.421468532207607216033551367218024097383022308349609375, 0.42146843802130329326161017888807691633701324462890625, -0.99999999999999989, 0.62292398855983019004972723654291189010479001808562316000461578369140625e-8).Mantissa);
  
             Parallel.ForEach (OperatorTests.Doubles(), x =>
             {
@@ -807,6 +807,24 @@ exp(x*x/4)*pcfu(0.5+n,-x)
             });
         }
 
+        // Same as MMath.NormalCdfIntegral but avoids inconsistent values of r and sqrtomr2 when using arbitrary precision.
+        public static ExtendedDouble NormalCdfIntegral(double x, double y, double r, double sqrtomr2)
+        {
+            if (sqrtomr2 < 0.618)
+            {
+                // In this regime, it is more accurate to compute r from sqrtomr2.
+                // See NormalCdfRatioLn
+                r = System.Math.Sign(r) * System.Math.Sqrt((1 - sqrtomr2) * (1 + sqrtomr2));
+            }
+            else
+            {
+                // In this regime, it is more accurate to compute sqrtomr2 from r.
+                double omr2 = 1 - r * r;
+                sqrtomr2 = System.Math.Sqrt(omr2);
+            }
+            return MMath.NormalCdfIntegral(x, y, r, sqrtomr2);
+        }
+
         internal void NormalCdfIntegralTest2()
         {
             double x = 0.0093132267868981222;

test/Tests/Data/SpecialFunctionsValues/NormalCdfLn2.csv

@@ -1,6 +1,6 @@
 arg0,arg1,arg2,expectedresult
 -Infinity,Infinity,0.5,-Infinity
--9945842.62906782515347003936767578125,9945822.06800303049385547637939453125,-0.9891958110248051383450729190371930599212646484375,-49459892801285.500323383881351469927253362271460929
+-9945842.62906782515347003936767578125,9945822.06800303049385547637939453125,-0.9891958110248051383450729190371930599212646484375,-49459892801108.425726931886363214188970559217641419
 -312498.3686245033168233931064605712890625,312498.2982211210182867944240570068359375,-0.9999893339082690513208717675297521054744720458984375,-48827615210.059272573857776712857226747641701183369
 -63,63,-0.46374946374946379723525069493916817009449005126953125,-1989.5623250537256917339818749521250583973045698991
 -20,-21,0.5,-287.59599647372995645990323264954196323654372871123

test/Tests/Data/SpecialFunctionsValues/logisticGaussianDeriv2.csv

@@ -117,7 +117,7 @@ arg0,arg1,expectedresult
 -0.1000000000000000055511151231257827021181583404541015625,0.79999999999999993338661852249060757458209991455078125,0.0069381918986754304427733173894720655777488875979834
 -0.1000000000000000055511151231257827021181583404541015625,0.90000000000000002220446049250313080847263336181640625,0.0065658283515188578988999317891629206501062940418998
 -0.1000000000000000055511151231257827021181583404541015625,1,0.0062289286083992064760315583970153118959636706295122
-0,1,1.6187994390755941681207614628194839967692740493738e-510
+0,1,0.0
 1,Infinity,0.0
 1.1999999999999999555910790149937383830547332763671875,2.5,-0.035252643710161957083496782969994794254464616517505
 1.930314781134816026764156049466691911220550537109375,0.32785733257471194601606612195610068738460540771484375,-0.077215104413877430947396768192089921234041942366547