New Interval.Contains method. Documentation. (#62)

Examples/Program.cs

@@ -20,59 +20,8 @@ private static MethodInfo[] GetExampleMethods () {
             return(methods);
         }
 
-        private static double EllipticPi(double n, double k) {
-
-            if (n > 1.0) throw new ArgumentOutOfRangeException(nameof(n));
-            if (k < 1.0 || k > 1.0) throw new ArgumentOutOfRangeException(nameof(k));
-
-            if (n == 1.0 || k == 1.0) return Double.PositiveInfinity;
-
-            // DLMF 19.8 describes how to compute \Pi(n, k) via AGM plus some auxiluary
-            // calculations. Here a and g, along with p,  converge to AGM(1,k') in the
-            // usual way; the sum of auxiluary variables q computed along the way gives \Pi(n, k).
-            // This method appears to have been derived by Carlson in "Three Improvements in
-            // Reduction and Computation of Elliptic Integrals", Journal of Research of the
-            // National Institute of Standards and Technology, 2002 Sep-Oct 107(5): 413-418
-            // (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4861378/)
-            // as a specialization of his method to calcuate R_J.
-
-            double a = 1.0;
-            double g = Math.Sqrt((1.0 - k) * (1.0 + k));
-            double n1 = 1.0 - n;
-            double pSquared = n1;
-            double p = Math.Sqrt(pSquared);
-            double q = 1.0;
-            double s = q;
-
-            for (int i = 1; i < 100; i++) {
-
-                double s_old = s;
-                double ag = a * g;
-                double e = (pSquared - ag) / (pSquared + ag);
-                q = 0.5 * q * e;
-                s += q;
-
-                double p_old = p;
-                p = (pSquared + ag ) / (2.0 * p);
-
-                if (p == p_old && s == s_old) {
-                    return Math.PI / 4.0 / p * (2.0 + n / n1 * s);
-                }
-
-                pSquared = p * p;
-                a = 0.5 * (a + g);
-                g = Math.Sqrt(ag);
-
-            }
-
-            throw new Exception();
-
-        }
-
-
         static void Main(string[] args)
         {
-            double e3 = EllipticPi(-0.25,0.99);
 
             MethodInfo[] methods = GetExampleMethods();
             Dictionary<string, MethodInfo> index = new Dictionary<string, MethodInfo>();
@@ -81,6 +30,7 @@ static void Main(string[] args)
             }
 
             if (args.Length == 0) {
+                Console.WriteLine("The following examples are defined:");
                 foreach(string key in index.Keys) {
                     Console.WriteLine(key);
                 }

Numerics/Core/Interval.cs

@@ -39,6 +39,46 @@ public double RightEndpoint {
             }
         }
 
+        /// <summary>
+        /// Determines whether a given point is contained in the closed interval.
+        /// </summary>
+        /// <param name="x">The point.</param>
+        /// <returns><see langword="true"/> if <paramref name="x"/> is contained in the closed interval [a,b],
+        /// otherwise <see langword="false"/>.</returns>
+        public bool Contains (double x) {
+            return (x >= a) && (x <= b);
+        }
+
+        /// <summary>
+        /// Determines whether a given point is contained in the closed or open interval.
+        /// </summary>
+        /// <param name="x">The point.</param>
+        /// <param name="type">Specifier of whether the interval should be considered closed or open.</param>
+        /// <returns><see langword="true"/> if <paramref name="x"/> is contained in the interval as specified,
+        /// otherwise <see langword="false"/></returns>
+        public bool Contains (double x, IntervalType type) {
+            if (type == IntervalType.Closed) {
+                return (x >= a) && (x <= b);
+            } else {
+                return (x > a) && (x < b);
+            }
+        }
+
+        /// <summary>
+        /// Determines whether a given point is contained in the interval, with the left and right endpoints
+        /// separately specififed as closed or open.
+        /// </summary>
+        /// <param name="x">The point.</param>
+        /// <param name="leftEndpointType">Specifier of whether the interval should be considered as including its left endpoint.</param>
+        /// <param name="rightEndpointType">Specifier of whether the interval should be considered as including its right endpoint.</param>
+        /// <returns><see langword="true"/> if <paramref name="x"/> is contained in the interval as specified,
+        /// otherwise <see langword="false"/></returns>
+        public bool Contains (double x, IntervalType leftEndpointType, IntervalType rightEndpointType) {
+            bool leftSatisfied = leftEndpointType == IntervalType.Closed ? x >= a : x > a;
+            bool rightSatisfied = rightEndpointType == IntervalType.Closed ? x <= b : x < b;
+            return leftSatisfied && rightSatisfied;
+        }
+
         /// <summary>
         /// Determines whether the argument lies in the open interval.
         /// </summary>
@@ -70,7 +110,7 @@ public double Width {
         /// </summary>
         public double Midpoint {
             get {
-                return ((a + b) / 2.0);
+                return 0.5 * (a + b);
             }
         }
 
@@ -214,4 +254,20 @@ internal IEnumerable<int> GetContainedIntegers () {
 
     }
 
+    /// <summary>
+    /// Indicates whether an interval should be considered closed or open.
+    /// </summary>
+    public enum IntervalType {
+
+        /// <summary>
+        /// Endpoints should be considered within the interval.
+        /// </summary>
+        Closed,
+
+        /// <summary>
+        /// Endpoints should be considered outside the interval.
+        /// </summary>
+        Open
+    }
+
 }

Numerics/Core/MoreMath.cs

@@ -140,7 +140,7 @@ public static double Hypot (double x, double y) {
         // Beebe, "Computation of expm1(x) = exp(x)  - 1", 2002 (http://www.math.utah.edu/~beebe/reports/expm1.pdf)
         // makes some good points about e^x - 1.
         //   * He shows that the point e^x = 1/2 and e^x = 3/2 are the relevent limits where Math.Exp(x) - 1.0
-        //     looses a one bit of accuracy.
+        //     looses one bit of accuracy.
         //   * He measures that the maximum number of terms in the Taylor series required in this region is 17.
         //   * He measures that the RMS error of the Taylor series in this region is ~0.8 bits and it's maximum
         //     relative error is ~ 2.7 bits.
@@ -177,9 +177,9 @@ private static double ReducedExpm1Series (double x) {
         /// <returns>The value of e<sup>x</sup>-1.</returns>
         /// <remarks>
         /// <para>If x is close to 0, then e<sup>x</sup> is close to 1, and computing e<sup>x</sup>-1 by
-        /// <tt>Math.Exp(x) - 1.0</tt> will be subject to severe loss of significance due to cancelation.
-        /// This method maintains full precision for all values of x by switching to a series expansion for values of
-        /// x near zero.</para>
+        /// <c>Math.Exp(x) - 1.0</c> will be subject to severe loss of significance due to cancelation.
+        /// This method maintains full precision for all values of x by switching to a series expansion
+        /// when x is near zero.</para>
         /// </remarks>
         public static double ExpMinusOne (double x) {
             if ((expm1SeriesLowerLimit < x) && (x < expm1SeriesUpperLimit)) {

Numerics/Extended/AdvancedDoubleDoubleMath.cs

@@ -10,20 +10,21 @@ namespace Meta.Numerics.Extended {
     /// <see cref="Double"/>-based functions, you must ensure that the arguments you provide are also more
     /// accurate than <see cref="Double"/>. Because decimal numbers expressed in code are automatically
     /// intrepreted by the compiler as <see cref="Double"/>, it's easier than you might think do this wrong.
-    /// For example, invoking Gamma(0.2) will <i>not</i> compute &#x393;(1/5) to <see cref="DoubleDouble"/>
+    /// For example, invoking <c>Gamma(0.2)</c> will <i>not</i> compute &#x393;(1/5) to <see cref="DoubleDouble"/>
     /// precision. The reason is that 0.2 is intrepreted by the compiler as a <see cref="Double"/>, and there
     /// is no <see cref="Double"/> value that is precisely 1/5. Instead, 0.2 parsed as a <see cref="Double"/>,
     /// is stored as 3602879701896397 X 2<sup>-54</sup> = 0.20000000000000001110223024625157... This equals
-    /// 0.2 within the 16 decimal-place accuracy of <see cref="Double"/>, but clearly not within the 32
+    /// 0.2 to within the 16 decimal-place accuracy of <see cref="Double"/>, but clearly not to within the 32
     /// decimal-place accuracy of <see cref="DoubleDouble"/>.</para>
     /// <para>There are a number of ways to ensure that you are providing an argument to <see cref="DoubleDouble"/> accuracy.
-    /// One possibility is to use the <see cref="DoubleDouble"/> text parser, for example by invoking new DoubleDouble("0.2").
-    /// Another is to produce tha argument as the result of calculation from exact integers, e.g. DoubleDouble.One / 5 or
-    /// ((DoubleDouble) 1) / 5. 
+    /// One possibility is to use the <see cref="DoubleDouble"/> text parser, for example by invoking <c>new DoubleDouble("0.2")</c>.
+    /// Another is to produce tha argument as the result of calculation from exact integers, <c>(DoubleDouble) 1 / 5</c>, which works
+    /// because 1 and 5 (like all integers within range) are represented exactly and the because of the cast the division operation
+    /// is <see cref="DoubleDouble.op_Division"/>.
     /// (The stored value in these cases is again not precisely 1/5, but is equal within the accuracy of <see cref="DoubleDouble"/>.)
     /// Finally, if you know that the argument you want is precisely represetable as a <see cref="Double"/>, you can safely
     /// use the compiler's parser. For example, invoking Gamma(0.25) does compute &#x393;(1/4) to <see cref="DoubleDouble"/>
-    /// to <see cref="DoubleDouble"/> accuracy, because 1/4 is exactly representable via <see cref="Double"/>.</para>
+    /// precision, because 1/4 is exactly representable by <see cref="Double"/>.</para>
     /// </remarks>
     public static partial class AdvancedDoubleDoubleMath {
 

Numerics/Extended/AdvancedDoubleDoubleMath_Gamma.cs

@@ -17,10 +17,12 @@ namespace Meta.Numerics.Extended {
     public static partial class AdvancedDoubleDoubleMath {
 
         /// <summary>
-        /// Computes the logarithm of the Gamma function double double precision.
+        /// Computes the logarithm of the Gamma function with double double precision.
         /// </summary>
         /// <param name="x">The argument, which must be non-negative.</param>
-        /// <returns>The value of Gamma(x).</returns>
+        /// <returns>The value of ln(&#x393;(x)).</returns>
+        /// <exception cref="ArgumentOutOfRangeException"><paramref name="x"/> is less than zero.</exception>
+        /// <seealso cref="Meta.Numerics.Functions.AdvancedMath.LogGamma(double)"/>
         public static DoubleDouble LogGamma (DoubleDouble x) {
 
             if (x < DoubleDouble.Zero) {

Numerics/Extended/DoubleDouble.cs

@@ -9,16 +9,38 @@ namespace Meta.Numerics.Extended {
     /// Represents a floating point number with quadruple precision.
     /// </summary>
     /// <remarks>
-    /// <para>The double double format uses two <see cref="Double"/> values to effectively
-    /// double the precision with which a number can be stored and manipulated as compared to
-    /// to the <see cref="Double"/> structure, i.e. to approximately 31 decimal digits of accuracy.</para>
+    /// <para>The <see cref="DoubleDouble"/> structure uses two <see cref="Double"/> values to achieve
+    /// twice the precision with which a floating point number can be stored and manipulated as compared to
+    /// to the <see cref="Double"/> structure, approximately 31 decimal digits.</para>
     /// <para>Of all the extended precision floating point systems, double double is the
     /// fastest when implemented in software. A typical floating point operation using
-    /// <see cref="DoubleDouble"/>s is just 3-4 times slower than on <see cref="Double"/>s.</para>
+    /// <see cref="DoubleDouble"/>s is just 3-4 times slower than using <see cref="Double"/>s.</para>
     /// <para>To instantiate a <see cref="DoubleDouble"/>, you can use <see cref="DoubleDouble.TryParse(string, out DoubleDouble)"/>,
     /// or <see cref="DoubleDouble.Parse(string)"/>, or the constructor <see cref="DoubleDouble.DoubleDouble(string)"/>
     /// to parse the text representation of the decimal value you want. If the value you want can be represented as a <see cref="Double"/>
-    /// or <see cref="Int32"/> or other built in type, you can cast that value to a <see cref="DoubleDouble"/>.</para>
+    /// or <see cref="Int32"/> or other built in numeric type, you can cast that value to a <see cref="DoubleDouble"/>.</para>
+    /// <para>When casting a <see cref="Double"/> to a <see cref="DoubleDouble"/>, there is a gotcha that you must be careful
+    /// to avoid. Suppose you write <c>DoubleDouble x = 0.2</c> or <c>DoubleDouble x = 1.0 / 5.0</c>. You might think that this produces the 
+    /// <see cref="DoubleDouble"/> representation of 1/5, but you would be wrong. The problem is that the compiler intreprets 0.2 or 1.0/5.0
+    /// as <see cref="Double"/>s, and 1/5th is not exactly representable as a double, since it is not a rational number with a power-of-two denominator.
+    /// Double's best attempt at 1/5th is 3602879701896397 X 2^<sup>-54</sup> = 0.20000000000000001110223024625157..., which
+    /// is accurate to 16 decimal digits, but not to 32. Therefore when it is cast to a <see cref="DoubleDouble"/> it is much
+    /// farther away from 1/5th than <see cref="DoubleDouble"/> can achieve. To obtain 1/5th to the accuracy of a <see cref="DoubleDouble"/>,
+    /// you must write <c>DoubleDouble x = new DoubleDouble("0.2")</c> or <c>DoubleDouble x = (DoubleDouble) 1 / 5</c>. (The latter works
+    /// because 1 and 5 are exactly representable and the division is performed as <see cref="DoubleDouble"/> division. All integers in range,
+    /// indeed all rational numbers with in-range numerators and power-of-two denominators, are exactly representable. So, for example,
+    /// <c>DoubleDouble x = 0.25</c> <i>does</i> work as expected, because 1/4 is exactly representable. But to avoid the gotcha
+    /// it's best to simply train yourself to avoid assigning <see cref="DoubleDouble"/> variables from factional <see cref="Double"/> values.)</para>
+    /// <para>Many of the mathematical functions which are implemented for <see cref="Double"/> arguments by static methods of the <see cref="Math"/> class
+    /// are implemented for <see cref="DoubleDouble"/> arguments by static methods of the <see cref="DoubleDouble"/> type itself,
+    /// for example <see cref="DoubleDouble.Sqrt(DoubleDouble)"/> and <see cref="DoubleDouble.Log(DoubleDouble)"/>.
+    /// Some of the advanced functions which are implemented for <see cref="Double"/> arguments by static methods of the <see cref="Meta.Numerics.Functions.AdvancedMath"/> class
+    /// are implemented for <see cref="DoubleDouble"/> arguments by static methods of the <see cref="AdvancedDoubleDoubleMath"/> class.
+    /// </para>
+    /// <para>You may wonder why <see cref="DoubleDouble"/> is not simply named "Quad". The reason is that "Quad" would propertly refer to an
+    /// implementation of the <see href="https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format">IEEE 754 quadruple-precision binary floating
+    /// point format</see>, which would have not only the extended presion of <see cref="DoubleDouble"/>, but also an extended range (up to 10<sup>4932</sup>).
+    /// </para>
     /// </remarks>
     public struct DoubleDouble : IEquatable<DoubleDouble>, IComparable<DoubleDouble> {