Dev/working (#60)

* Improve special function performance and increase resiliance, especially to NaNs. Improve Poisson and Beta random deviate performance.

* UInt128, Int128 working. DoubleDouble improved. Complete Elliptic Inegral of 3rd kind. Improved Gamma speed and accuracy using series aound 1 and 2. Explicit ScaledModifiedBessel and SphericalBessel.

* Mostly tests and documentation. RegressionResult layer.

* Add tests, remove unused code, add increment and decriment operators to Int128 and UInt128.

* Int128 % and cast from double. Documentation improvements.

Numerics/Extended/Int128.cs

@@ -258,6 +258,23 @@ public static explicit operator Int128 (BigInteger b) {
             return new Int128(u);
         }
 
+        /// <summary>
+        /// Converts a floating-point value into a 128-bit integer.
+        /// </summary>
+        /// <param name="x">A floating-point number.</param>
+        /// <returns>The lower 128 bits of the integer part of the floating-point number.</returns>
+        /// <exception cref="InvalidCastException"><paramref name="x"/> is NaN, or infinite.</exception>
+        public static explicit operator Int128 (double x) {
+            DoubleInfo s = new DoubleInfo(Math.Truncate(x));
+            if (!s.IsFinite) throw new InvalidCastException();
+            int e = s.Exponent;
+            if (e < 0) return Int128.Zero;
+            UInt128 u = (ulong) s.Mantissa;
+            u = u << e;
+            if (s.IsNegative) u = u.Negate();
+            return new Int128(u);
+        }
+
         // Serialization and deserialization
 
         /// <summary>
@@ -339,6 +356,28 @@ public static bool TryParse (string text, out Int128 value) {
             return new Int128(x.u - y.u);
         }
 
+        /// <summary>
+        /// Increments a 128-bit integer.
+        /// </summary>
+        /// <param name="x">The integer.</param>
+        /// <returns>One more than <paramref name="x"/>.</returns>
+        public static Int128 operator ++ (Int128 x) {
+            UInt128 v = x.u;
+            v++;
+            return new Int128(v);
+        }
+
+        /// <summary>
+        /// Decrements a 128-bit unsigned integer.
+        /// </summary>
+        /// <param name="x">The integer.</param>
+        /// <returns>One less than <paramref name="x"/>.</returns>
+        public static Int128 operator -- (Int128 x) {
+            UInt128 v = x.u;
+            v--;
+            return new Int128(v);
+        }
+
         // Multiplication
 
         /// <summary>
@@ -416,6 +455,17 @@ public static Int128 DivRem (Int128 x, Int128 y, out Int128 r) {
 
         }
 
+        /// <summary>
+        /// Computes the remainder of two 128 bit integers.
+        /// </summary>
+        /// <param name="x">The dividend.</param>
+        /// <param name="y">The divisor.</param>
+        /// <returns>The remainder of <paramref name="x"/> divided by <paramref name="y"/>.</returns>
+        public static Int128 operator % (Int128 x, Int128 y) {
+            Int128.DivRem(x, y, out Int128 r);
+            return r;
+        }
+
         // Absolute Value
 
         /// <summary>

Numerics/Extended/Int128Calculator.cs

@@ -5,7 +5,6 @@ namespace Meta.Numerics.Extended {
 
 
     // These routines are adapted from
-
     // Warren, Henry, "Hacker's Delight", 2nd edition
 
 
@@ -35,15 +34,15 @@ public static void Subtract128From128 (ulong x1, ulong x0, ulong y1, ulong y0, o
             z1 = x1 - y1 - b;
         }
 
-        public static void Increment128 (ref ulong x1, ref ulong x0) {
-            x0++;
-            if (x0 == 0UL) x1++; 
+        public static void Increment128 (ulong x1, ulong x0, out ulong y1, out ulong y0) {
+            y0 = x0 + 1UL;
+            // This is a specialization of the trick above to avoid branching on (y0 == 0UL).
+            ulong c = (x0 & ~y0) >> 63;
+            y1 = x1 + c;
         }
 
         public static void TwosComplement (ulong x1, ulong x0, out ulong y1, out ulong y0) {
-            y1 = ~x1;
-            y0 = ~x0;
-            Increment128(ref y1, ref y0);
+            Increment128(~x1, ~x0, out y1, out y0);
         }
 
         public static void Decompose (ulong u, out ulong u1, out ulong u0) {

Numerics/Extended/UInt128.cs

@@ -170,7 +170,7 @@ public static explicit operator UInt128 (BigInteger b) {
         /// <returns>The lower 128 bits of the integer part of the floating-point number.</returns>
         /// <exception cref="InvalidCastException"><paramref name="x"/> is negative, or NaN, or infinite.</exception>
         public static explicit operator UInt128 (double x) {
-            DoubleInfo s = new DoubleInfo(Math.Floor(x));
+            DoubleInfo s = new DoubleInfo(Math.Truncate(x));
             if (s.IsNegative || !s.IsFinite) throw new InvalidCastException();
             int e = s.Exponent;
             if (e < 0) return UInt128.Zero;
@@ -248,11 +248,7 @@ public override int GetHashCode () {
         internal bool IsNegative => ((s1 >> 63) != 0UL);
 
         internal UInt128 Negate () {
-            //ulong y1 = ~s1;
-            //ulong y0 = ~s0;
-            //Int128Calculator.Increment128(ref y1, ref y0);
-            //return new UInt128(y1, y0);
-            Int128Calculator.Add128To128(~s1, ~s0, 0UL, 1UL, out ulong y1, out ulong y0);
+            Int128Calculator.TwosComplement(s1, s0, out ulong y1, out ulong y0);
             return new UInt128(y1, y0);
         }
 
@@ -343,6 +339,16 @@ public int CompareTo (UInt128 u) {
             return new UInt128(s1, s0);
         }
 
+        /// <summary>
+        /// Increments a 128-bit unsigned integer.
+        /// </summary>
+        /// <param name="x">The integer.</param>
+        /// <returns>One more than <paramref name="x"/> (or <see cref="UInt128.Zero"/>, if <paramref name="x"/> is <see cref="UInt128.MaxValue"/>).</returns>
+        public static UInt128 operator ++ (UInt128 x) {
+            Int128Calculator.Increment128(x.s1, x.s0, out ulong s1, out ulong s0);
+            return new UInt128(s1, s0);
+        }
+
         /// <summary>
         /// Subtracts one 128 bit unsigned integer from another.
         /// </summary>
@@ -354,6 +360,16 @@ public int CompareTo (UInt128 u) {
             return new UInt128(s1, s0);
         }
 
+        /// <summary>
+        /// Decrements a 128-bit unsigned integer.
+        /// </summary>
+        /// <param name="x">The integer.</param>
+        /// <returns>One less than <paramref name="x"/> (or <see cref="UInt128.MaxValue"/>, if <paramref name="x"/> is <see cref="UInt128.Zero"/>).</returns>
+        public static UInt128 operator-- (UInt128 x) {
+            Int128Calculator.Subtract128From128(x.s1, x.s0, 0UL, 1UL, out ulong s1, out ulong s0);
+            return new UInt128(s1, s0);
+        }
+
         // Multiplication
 
         /// <summary>

Numerics/Functions/AdvancedMath_Airy.cs

@@ -379,6 +379,7 @@ private static SolutionPair Airy_Asymptotic_Negative (double x) {
         /// <param name="k">The index of the zero.</param>
         /// <returns>The <paramref name="k"/>th value of x for which Ai(x) = 0.</returns>
         /// <exception cref="ArgumentOutOfRangeException"><paramref name="k"/> is less than 1.</exception>
+        /// <seealso cref="AiryAi(double)"/>
         public static double AiryAiZero (int k) {
             return (BesselMath.AiryAiZero(k));
         }
@@ -389,6 +390,7 @@ public static double AiryAiZero (int k) {
         /// <param name="k">The index of the zero.</param>
         /// <returns>The <paramref name="k"/>th value of x for which Bi(x) = 0.</returns>
         /// <exception cref="ArgumentOutOfRangeException"><paramref name="k"/> is less than 1.</exception>
+        /// <seealso cref="AiryBi(double)"/>
         public static double AiryBiZero (int k) {
             return (BesselMath.AiryBiZero(k));
         }

Numerics/Functions/AdvancedMath_Bessel.cs

@@ -19,6 +19,7 @@ public static partial class AdvancedMath {
         /// <para>For information on the cylindrical Bessel functions, see <see cref="AdvancedMath.Bessel"/>.</para>
         /// </remarks>
         /// <seealso href="http://en.wikipedia.org/wiki/Bessel_function"/>
+        /// <seealso href="https://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html"/>
         public static double BesselJ (int n, double x) {
 
             // Relate negative n to positive n.
@@ -82,6 +83,7 @@ public static double BesselJ (int n, double x) {
         /// <para>For information on the cylindrical Bessel functions, see <see cref="AdvancedMath.Bessel"/>.</para>
         /// </remarks>
         /// <seealso href="http://en.wikipedia.org/wiki/Bessel_function"/>
+        /// <seealso href="https://mathworld.wolfram.com/BesselFunctionoftheSecondKind.html"/>
         public static double BesselY (int n, double x) {
 
             // Relate negative n to positive n.
@@ -563,6 +565,7 @@ private static Complex Bessel_CF2 (double nu, double x) {
         /// </remarks>
         /// <exception cref="ArgumentOutOfRangeException"><paramref name="x"/> is negative.</exception>
         /// <seealso href="http://en.wikipedia.org/wiki/Bessel_function"/>
+        /// <seealso href="https://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html"/>
         public static double BesselJ (double nu, double x) {
 
             if (x < 0.0) throw new ArgumentOutOfRangeException(nameof(x));
@@ -631,6 +634,7 @@ public static double BesselJ (double nu, double x) {
         /// <exception cref="ArgumentOutOfRangeException"><paramref name="x"/> is negative.</exception>
         /// <seealso cref="Bessel"/>
         /// <seealso href="http://en.wikipedia.org/wiki/Bessel_function"/>
+        /// <seealso href="https://mathworld.wolfram.com/BesselFunctionoftheSecondKind.html"/>
         public static double BesselY (double nu, double x) {
 
             if (x < 0.0) throw new ArgumentOutOfRangeException(nameof(x));
@@ -989,7 +993,7 @@ private static SolutionPair Bessel_Zero (double nu, int s) {
 
         // **** Spherical Bessel functions ****
 
-        
+
 
         // Functions needed for uniform asymptotic expansions
 
@@ -998,7 +1002,7 @@ private static SolutionPair Bessel_Zero (double nu, int s) {
         /// </summary>
         /// <param name="nu">The order, which must be non-negative.</param>
         /// <param name="k">The index of the zero, which must be positive.</param>
-        /// <returns>The <paramref name="k"/>th value of x for which J<sub>nu</sub>(x) = 0.</returns>
+        /// <returns>The <paramref name="k"/>th value of x for which J<sub>&#x3BD;</sub>(x) = 0.</returns>
         /// <exception cref="ArgumentOutOfRangeException"><paramref name="nu"/> is negative or <paramref name="k"/> is non-positive.</exception>
         public static double BesselJZero (double nu, int k) {
             if (nu < 0.0) throw new ArgumentOutOfRangeException(nameof(nu));
@@ -1062,22 +1066,6 @@ private static double ApproximateBesselJZero_LargeOrder (double nu, int k) {
             double z = ZFromZeta(zeta, out double c1);
             return (z * (nu + c1 / nu));
         }
-
-        private static double ZetaFromZ (double z) {
-            Debug.Assert(z > 0.0);
-            if (z < 0.75) {
-                double s = Math.Sqrt((1.0 - z) * (1.0 + z));
-                return (Math.Pow(3.0 / 2.0 * (Math.Log((1.0 + s) / z) - s), 2.0 / 3.0));
-            } else if (z < 1.25) {
-                double y = 1.0 - z;
-                double c = Math.Pow(2.0, 1.0 / 3.0);
-                // Need more terms
-                return (c * y * (1.0 + 3.0 / 10.0 * y + 32.0 / 175.0 * y * y + 1037.0 / 7875 * y * y * y));
-            } else {
-                double s = Math.Sqrt((z - 1.0) * (z + 1.0));
-                return (-Math.Pow(3.0 / 2.0 * (s - Math.Acos(1.0 / z)), 2.0 / 3.0));
-            }
-        }
 
         // This is an inversion of the zeta-from-z function.
         // Since it is used only in the approximate root expressions, we

Numerics/Functions/AdvancedMath_Exponential.cs

@@ -227,11 +227,14 @@ public static double IntegralCi (double x) {
         /// <param name="x">The argument.</param>
         /// <returns>The value of Cin(x).</returns>
         /// <remarks>
-        /// <para>The entine cosine integral Cin(x) is related to the conventional cosine integral Ci(x) by:</para>
+        /// <para>The entire cosine integral Cin(x) is related to the conventional cosine integral Ci(x) by:</para>
+        /// <img src="../images/IntegralCin.png"/>
         /// <para>Unlike Ci(x), Cin(x) is regular at the origin, and may therefore be more useful in applications
-        /// that need the non-divergent part of a cosine integral. But, again unlike Ci(x), Cin(x) does diverge
-        /// (logarithmicaly) for large x.</para>
+        /// that need the non-divergent part of a cosine integral. In fact, Cin(x) is entire, meaning it has no poles
+        /// or cuts anywhere in the complex plane. But, unlike Ci(x), Cin(x) does diverge (logarithmicaly) for large x.</para>
         /// </remarks>
+        /// <seealso href="https://en.wikipedia.org/wiki/Trigonometric_integral"/>
+        /// <seealso href="https://dlmf.nist.gov/6"/>
         /// <seealso cref="IntegralCi"/>
         public static double IntegralCin (double x) {
             if (x < 0.0) {
@@ -278,6 +281,8 @@ public static double IntegralSi (double x) {
         /// </summary>
         /// <param name="x">The argument.</param>
         /// <returns>The value of Shi(x).</returns>
+        /// <seealso href="https://en.wikipedia.org/wiki/Trigonometric_integral#Hyperbolic_sine_integral"/>
+        /// <seealso href="https://mathworld.wolfram.com/Shi.html"/>
         public static double IntegralShi (double x) {
             if (x < 0.0) {
                 return (-IntegralShi(-x));
@@ -300,7 +305,6 @@ public static double IntegralShi (double x) {
         // This is written so as to be usable for both Si and Shi.
         // Converges to full accuracy in about 30 terms for x~4, less for smaller x
         private static double IntegralSi_Series (double xSquared) {
-            //double xx = x * x;
             double dy = 1.0;
             double y = dy;
             for (int k=3; k<Global.SeriesMax; k+= 2) {

Numerics/Functions/AdvancedMath_Gamma.cs

@@ -842,23 +842,6 @@ public static Complex LogGamma (Complex z) {
             }
         }
 
-        private static Complex LogGammaReflectionTerm (Complex z) {
-
-            // s will overflow for large z.Im, but we are taking its log, so we should
-            // re-formulate this 
-            Complex s = new Complex(MoreMath.SinPi(z.Re) * Math.Cosh(Math.PI * z.Im), MoreMath.CosPi(z.Re) * Math.Sinh(Math.PI * z.Im));
-            Complex logs = ComplexMath.Log(s);
-
-            double m1 = Math.Log(MoreMath.Hypot(MoreMath.SinPi(z.Re), Math.Sinh(Math.PI * z.Im)));
-
-            double t = Math.Abs(Math.PI * z.Im);
-            double c1 = Math.Exp(-2.0 * t);
-            double c2 = MoreMath.SinPi(z.Re) / Math.Sinh(t);
-            double m2 = t + MoreMath.LogOnePlus(-c1) + 0.5 * MoreMath.LogOnePlus(c2 * c2) - Math.Log(2.0);
-
-            return (Math.Log(Math.PI) - logs);
-        }
-
         private static Complex LogGamma_Stirling (Complex z) {
 
             // work in the upper complex plane; i think this isn't actually necessary

Numerics/Functions/AdvancedMath_Harmonic.cs

@@ -15,6 +15,7 @@ public static partial class AdvancedMath {
         /// <param name="phi">The cylindrical angle &#x3C6;. This angle is usually expressed as between 0 and 2&#x3C0;, measured counter-clockwise (as seen from above) from the positive x-axis. It is also possible to use negative values to represent clockwise movement. </param>
         /// <returns>The value of Y<sub>l,m</sub>(&#x3B8;,&#x3C6;).</returns>
         /// <exception cref="ArgumentOutOfRangeException"><paramref name="l"/> is negative, or <paramref name="m"/> lies outside the range [-l, l].</exception>
+        /// <seealso href="https://mathworld.wolfram.com/SphericalHarmonic.html"/>
         public static Complex SphericalHarmonic (int l, int m, double theta, double phi) {
             if (l < 0) throw new ArgumentOutOfRangeException(nameof(l));
             if ((m > l) || (m < -l)) throw new ArgumentOutOfRangeException(nameof(m));

Numerics/Functions/AdvancedMath_Hypergeometric.cs

@@ -18,7 +18,7 @@ public static partial class AdvancedMath {
         /// <returns>The value of <sub>2</sub>F<sub>1</sub>(a, b; c; x).</returns>
         /// <remarks>
         /// <para>The Gauss Hypergeometric function is defined by the hypergeometric series:</para>
-        /// <img src="Hypergeometric2F1.png" />
+        /// <img src="../images/Hypergeometric2F1.png" />
         /// <para>For generic values of a, b, and c, the Gauss hypergeometric function becomes complex for x > 1.
         /// However, there are specific cases, most commonly for negative integer values of a and b, for which the
         /// function remains real in this range.</para>

Numerics/Functions/AdvancedMath_Lambert.cs

@@ -15,8 +15,10 @@ public static partial class AdvancedMath {
         /// The function appears in a number of contexts, including the solution of differential
         /// equations and the enumeration of trees.</para>
         /// </remarks>
+        /// <exception cref="ArgumentOutOfRangeException"><paramref name="x"/> is less than -1/e.</exception>
         /// <seealso href="http://en.wikipedia.org/wiki/Lambert_W_function" />
-        /// <seealso href="http://www.apmaths.uwo.ca/~djeffrey/Offprints/W-adv-cm.pdf"/>
+        /// <seealso href="https://mathworld.wolfram.com/LambertW-Function.html"/>
+        /// <seealso href="https://dlmf.nist.gov/4.13"/>
         public static double LambertW (double x) {
 
             if (x < -EI) throw new ArgumentOutOfRangeException(nameof(x));