Dev/start5 (#68)

* Change to .NET Standard 2.0. Make structs readonly. Start to rip out statistics container classes. Add ChebyshevU. Add Benford and SkewNormal distributions.

* Beta use unit interval. no Sample, and other cleanup.

* Many mostly small improvements to accuracy and added tests.

* Final changes to prep for release.

Author: dcwuser

Numerics/Core/ComplexMath.cs

@@ -592,5 +592,38 @@ public static Complex Pow (Complex z, int n) {
 
         }
 
+
+
+
+        /// <summary>
+        /// Computes the nth roots of unity.
+        /// </summary>
+        /// <param name="n">The order of roots to be computed, which must be positive.</param>
+        /// <returns>The n complex numbers x such that x<sup>n</sup> = 1.</returns>
+        /// <remarks>
+        /// <para>The roots are listed in the usual anti-clockwise order, starting from 1.</para>
+        /// <para>This method ensures that 1 and -1 (when it appears) are exact, and that conjugate roots are exactly conjugate.</para>
+        /// </remarks>
+        public static Complex[] RootsOfUnity(int n) {
+            if (n < 1) throw new ArgumentOutOfRangeException(nameof(n));
+            Complex[] roots = new Complex[n];
+            roots[0] = Complex.One;
+            double delta = Global.TwoPI / n;
+            int half_n = (n + 1) / 2; // This computes n/2 rounded up rather than down.
+            for (int i = 1; i < half_n; i++) {
+                double theta = delta * i;
+                roots[i] = new Complex(MoreMath.Cos(theta), MoreMath.Sin(theta));
+                roots[n - i] = roots[i].Conjugate;
+            }
+            if (n % 2 == 0) roots[half_n] = -Complex.One;
+            return roots;
+        }
+
+        // This computation of the complex roots of unity is careful to ensure
+        //   +1 and -1 (when it appears) are always exact
+        //   conjugate roots are always exactly conjugate
+        // Using CosPi and SinPi is not necessary since the arguments are never reducable.
+        // I notice for n=8 not all components are exactly equal; look into this.
+
     }
 }

Numerics/Core/DiscreteInterval.cs

@@ -63,7 +63,7 @@ public int RightEndpoint {
         /// This is one less than the number of integers in the interval. Thus an interval with equal left and right endpoints
         /// has width 0, not width 1. And and the interval {0 .. MaxValue} has width MaxValue, not MaxValue + 1.</para>
         /// </remarks>
-        internal uint Width {
+        public uint Width {
             get {
                 return (uint) (max - min);
             }

Numerics/Core/MoreMath.cs

@@ -238,7 +238,7 @@ internal static double ReducedExpMinusOne (double x, double e) {
         /// <param name="x">The argument.</param>
         /// <returns>The value of log(1+x).</returns>
         /// <remarks>
-        /// <para>If x is close to 0, computing log(1+x) by <tt>Math.Log(1.0 + x)</tt> can result in a
+        /// <para>If x is close to 0, computing log(1+x) by <c>Math.Log(1.0 + x)</c> can result in a
         /// significant loss of accuracy.
         /// This function maintains full precision of all values of x by switching to a series expansion
         /// for values of x near zero.</para>
@@ -416,9 +416,8 @@ public static double Cos (double x) {
         /// that arise from the finite precision of the stored constant <see cref="Math.PI"/>.</para>
         /// <para>Suppose, for example, x = 1.0E6. Since x is an integer, sin(&#x3C0;x) = 0.0.
         /// However, due to the finite accuracy of Math.PI, Math.PI * x is not a perfect multiple of &#x3C0;, and
-        /// Math.Sin(Math.PI * x) = -2.2318717360358953E-10. But MoreMath.SinPi(x) = 0.0 exactly. Even for
-        /// arguments that are not exact integers, the accurary of MoreMath.SinPi will be better (because
-        /// it is possible to do argument reductions by integers exactly.)
+        /// <c>Math.Sin(Math.PI * x)</c> = -2.2318717360358953E-10. But <c>MoreMath.SinPi(x)</c> = 0.0 exactly. Even for
+        /// arguments that are not exact integers, the accurary of MoreMath.SinPi will be better.
         /// </para>
         /// </remarks>
         /// <seealso cref="CosPi(double)"/>

Numerics/Core/NamespaceDoc.cs

@@ -0,0 +1,11 @@
+using System.Runtime.CompilerServices;
+
+namespace Meta.Numerics {
+
+    /// <summary>
+    /// Contains types used throughout the Meta.Numerics library.
+    /// </summary>
+    [CompilerGenerated]
+    internal class NamespaceDoc {
+    }
+}

Numerics/Core/Point.cs

@@ -8,7 +8,7 @@ namespace Meta.Numerics {
     /// <summary>
     /// Represents a two-dimensional point.
     /// </summary>
-    public struct XY : IEquatable<XY> {
+    public readonly struct XY : IEquatable<XY> {
 
         /// <summary>
         /// Initializes a new point with the given coordinates.
@@ -20,7 +20,7 @@ public XY (double x, double y) {
             this.y = y;
         }
 
-        private double x, y;
+        private readonly double x, y;
 
         /// <summary>
         /// Gets the X-coordinate of the point.

Numerics/Data/FrameView.cs

@@ -50,7 +50,7 @@ public enum SortOrder {
     /// <see cref="this[string]"/> accessor to get a column, together with the <see cref="FrameColumn.As{T}"/>
     /// caster to expose it as a collection of the required type.
     /// For example, to obtain a estimate of the mean of the population from the sample in the
-    /// column named "heights", write <tt>view["height"].As&lt;double&gt;().PopulationMean()</tt>.
+    /// column named "heights", write <c>view["height"].As&lt;double&gt;().PopulationMean()</c>.
     /// Note that, for this to succeed, the underlying storage type of the heights column need not be double. As
     /// long as the data are convertible to the target type, no problems will arise. For example,
     /// the underlying storage type might be int, or double? as long as no null values are present in the view.</para>

Numerics/Functions/AdvancedMath_Elliptic.cs

@@ -558,6 +558,13 @@ public static double EllipticPi (double n, double k) {
         /// </summary>
         /// <param name="k">The elliptic modulus.</param>
         /// <returns>The corresponding Jacobi nome q.</returns>
+        /// <remarks>
+        /// <para>The elliptic nome q(k) = e<sup>&#x3C0; K'/ K</sup>, where K is the elliptic integral of the first kind
+        /// for the modulus k, and K' is the elliptic integral of the first kind for its complement k'.</para>
+        /// </remarks>
+        /// <seealso cref="EllipticK(double)"/>
+        /// <seealso href="https://en.wikipedia.org/wiki/Nome_(mathematics)"/>
+        /// <seealso href="https://mathworld.wolfram.com/Nome.html"/>
         public static double EllipticNome (double k) {
 
             if (k < 0.0 || k > 1.0) throw new ArgumentOutOfRangeException(nameof(k));

Numerics/Functions/AdvancedMath_Exponential.cs

@@ -1,5 +1,6 @@
 using System;
 using System.Diagnostics;
+using static System.Net.WebRequestMethods;
 
 
 namespace Meta.Numerics.Functions {
@@ -260,6 +261,8 @@ public static double IntegralCin (double x) {
         /// </summary>
         /// <param name="x">The argument.</param>
         /// <returns>The value of Chi(x).</returns>
+        /// <seealso href="https://en.wikipedia.org/wiki/Trigonometric_integral#Hyperbolic_cosine_integral"/>
+        /// <seealso href="https://mathworld.wolfram.com/Chi.html"/>
         public static double IntegralChi (double x) {
             if (x < 0.0) {
                 throw new ArgumentOutOfRangeException(nameof(x));
@@ -285,6 +288,7 @@ public static double IntegralChi (double x) {
         /// <img src="../images/SiIntegral.png" />
         /// <para>The sine integral is zero at the origin and executes a damped oscillation around &#x3C0;/2 as its argument increases.</para>
         /// </remarks>
+        /// <seealso cref="IntegralLittleSi(double)"/>
         /// <seealso href="https://en.wikipedia.org/wiki/Trigonometric_integral"/>
         /// <seealso href="http://mathworld.wolfram.com/SineIntegral.html"/>
         /// <seealso href="https://dlmf.nist.gov/6"/>
@@ -308,6 +312,15 @@ public static double IntegralSi (double x) {
         /// </summary>
         /// <param name="x">The argument.</param>
         /// <returns>The value of si(x).</returns>
+        /// <remarks>
+        /// <para>The value of the little sine integral is offset from the value of the big sine integral by &#x3C0;/2.
+        /// Thus at the orgin it has the value -&#x3C0;/2, but asymptotically it oscialtes around zero. Thus its
+        /// better to use this function for large values when you want to accurate know the amplitude of these oscilations.</para>
+        /// </remarks>
+        /// <seealso cref="IntegralSi(double)"/>
+        /// <seealso href="https://en.wikipedia.org/wiki/Trigonometric_integral"/>
+        /// <seealso href="http://mathworld.wolfram.com/SineIntegral.html"/>
+        /// <seealso href="https://dlmf.nist.gov/6"/>
         public static double IntegralLittleSi (double x) {
             if (x < 0.0) {
                 throw new ArgumentOutOfRangeException(nameof(x));

Numerics/Functions/OrthogonalPolynomials.cs

@@ -353,7 +353,7 @@ public static double ChebyshevT (int n, double x) {
         /// <remarks>
         /// <para>Chebyshev polynomials of the second kind are orthogonal on the interval [-1,1] with the weight (1-x<sup>2</sup>)<sup>1/2</sup>.</para>
         /// <para>Values returned are fully accurate (14-16 decimal digits) over the full range of argument for orders up to thousands.
-        /// By orders up to a million, ~11 decimal digits remain accurate.</para>
+        /// For orders up to a million, ~11 decimal digits remain accurate.</para>
         /// </remarks>
         /// <exception cref="ArgumentOutOfRangeException"><paramref name="n"/> is negative, or <paramref name="x"/> lies outside [-1,+1].</exception>
         /// <seealso href="https://en.wikipedia.org/wiki/Chebyshev_polynomials"/>

Numerics/Functions/Permutation.cs

@@ -2,6 +2,7 @@
 // Copyright 2014 by David Wright.
 // All rights reserved.
 
+using Meta.Numerics.Matrices;
 using System;
 using System.Collections.Generic;
 using System.Diagnostics;
@@ -119,8 +120,7 @@ public static bool TryParse (string text, out PermutationAsCycles result) {
                 int[] cycle = new int[cycleText.Length];
                 dimension += cycle.Length;
                 for (int k = 0; k < cycle.Length; k++) {
-                    int value;
-                    if (!Int32.TryParse(cycleText[k], out value)) return (false);
+                    if (!Int32.TryParse(cycleText[k], out int value)) return (false);
                     cycle[k] = value;
                 }
                 cyclesList.Add(cycle);
@@ -217,8 +217,7 @@ public static bool TryParse (string text, out PermutationAsMap result) {
             int[] map = new int[mapText.Length];
             bool[] flags = new bool[mapText.Length];
             for (int i = 0; i < map.Length; i++) {
-                int value;
-                if (!Int32.TryParse(mapText[i], out value)) return (false);
+                if (!Int32.TryParse(mapText[i], out int value)) return (false);
                 if ((value < 0) || (value > (map.Length - 1))) return (false);
                 if (flags[value]) return (false);
                 map[i] = value;
@@ -411,9 +410,8 @@ public string ToString (string format, IFormatProvider formatProvider) {
         /// <exception cref="FormatException"><paramref name="text"/> is not a valid text representation of a permutation.</exception>
         public static Permutation Parse (string text) {
             if (text == null) throw new ArgumentNullException(nameof(text));
-            Permutation output;
-            if (TryParse(text, out output)) {
-                return (output);
+            if (TryParse(text, out Permutation output)) {
+                return output;
             } else {
                 throw new FormatException();
             }
@@ -437,13 +435,11 @@ public static bool TryParse (string text, out Permutation output) {
             bool success = false;
             if ((text.Length > 0) && (text[0] == '[')) {
                 // try to parse as ordering
-                PermutationAsMap map;
-                success = PermutationAsMap.TryParse(text, out map);
+                success = PermutationAsMap.TryParse(text, out PermutationAsMap map);
                 if (success) output = new Permutation(map);
             } else {
                 // try to parse as cycles
-                PermutationAsCycles cycles;
-                success = PermutationAsCycles.TryParse(text, out cycles);
+                success = PermutationAsCycles.TryParse(text, out PermutationAsCycles cycles);
                 if (success) output = new Permutation(cycles);
             }
 
@@ -494,6 +490,19 @@ private void ComputeMap () {
             map = new PermutationAsMap(FromCyclesToMap(cycles.cycles));
         }
 
+        internal int[] Map {
+            get {
+                if (map == null) ComputeMap();
+                return map.map;
+            }
+        }
+
+        internal int[][] Cycles {
+            get {
+                if (cycles == null) ComputeCycles();
+                return cycles.cycles;
+            }
+        }
 
         /// <summary>
         /// Applies the permutation to a list.
@@ -798,6 +807,13 @@ public static Permutation GetRandomPermutation (int dimension, Random rng) {
             return (new Permutation(PermutationAsMap.GetRandomPermutation(dimension, rng)));
         }
 
+        /// <summary>
+        /// Get a matrix representation of the permuation.
+        /// </summary>
+        /// <returns>A matrix representation of the permutation.</returns>
+        public PermutationMatrix ToMatrix () {
+            return new PermutationMatrix(this);
+        }
     }
 
 }