diff --git a/Numerics/Analysis/MultiFunctionMath_Extrema_ModelTrust.cs b/Numerics/Analysis/MultiFunctionMath_Extrema_ModelTrust.cs index e420497..5c4b6ae 100644 --- a/Numerics/Analysis/MultiFunctionMath_Extrema_ModelTrust.cs +++ b/Numerics/Analysis/MultiFunctionMath_Extrema_ModelTrust.cs @@ -1,757 +1,742 @@ -using System; -using System.Collections.Generic; -using System.Diagnostics; - -using Meta.Numerics.Functions; -using Meta.Numerics.Matrices; - -namespace Meta.Numerics.Analysis { - - - - public static partial class MultiFunctionMath { - - /// - /// Finds a local minimum of a multi-dimensional function in the vicinity of the given starting location. - /// - /// The multi-dimensional function to minimize. - /// The starting location for the search. - /// The local minimum. - /// - /// For two dimensions, the default evaluation settings target a precision of about 10-8 - /// and allow up to about 1000 function evaluations. For increasing dimensions, these requirements are - /// gradually relaxed. - /// - /// , is . - /// The number of function evaluations required exceeded the evaluation budget. - public static MultiExtremum FindLocalMinimum (Func, double> function, IReadOnlyList start) { - return (FindLocalMinimum(function, start, new MultiExtremumSettings())); - } - - private static void SetDefaultOptimizationSettings (MultiExtremumSettings settings, int d) { - if (settings.RelativePrecision < 0.0) settings.RelativePrecision = Math.Pow(10.0, -(10.0 + 4.0 / d)); - if (settings.AbsolutePrecision < 0.0) settings.AbsolutePrecision = Math.Pow(10.0, -(10.0 + 4.0 / d)); - if (settings.EvaluationBudget < 0) settings.EvaluationBudget = 16 * (d + 1) * (d + 2) * (d + 3); - } - - /// - /// Finds a local minimum of a multi-dimensional function in the vicinity of the given starting location, - /// subject to the given evaluation constraints. - /// - /// The multi-dimensional function to minimize. - /// The starting location for the search. - /// The evaluation settings that govern the search for the minimum. - /// The local minimum. - /// - /// The Hessian (matrix of second derivatives) returned with the minimum is an approximation that is constructed in the course of search. - /// It should be considered a crude approximation, and may not even be that if the minimum is highly non-quadratic. - /// If you have a constrained minimization problem, require a high-precision solution, and do not have a good initial guess, consider first feeding - /// your constrained problem into , - /// which supports constraints but gives relatively lower precision solutions, then - /// feeding the result of that method into this method, which finds relatively precise solutions but does not support constraints. - /// - /// , , or is . - /// The number of function evaluations required exceeded the evaluation budget. - public static MultiExtremum FindLocalMinimum (Func, double> function, IReadOnlyList start, MultiExtremumSettings settings) { - if (function == null) throw new ArgumentNullException(nameof(function)); - if (start == null) throw new ArgumentNullException(nameof(start)); - if (settings == null) throw new ArgumentNullException(nameof(settings)); - return (FindLocalExtremum(function, start, settings, false)); - } - - /// - /// Finds a local maximum of a multi-dimensional function in the vicinity of the given starting location. - /// - /// The multi-dimensional function to maximize. - /// The starting location for the search. - /// The local maximum. - /// , is . - /// The number of function evaluations required exceeded the evaluation budget. - public static MultiExtremum FindLocalMaximum (Func, double> function, IReadOnlyList start) { - return (FindLocalMaximum(function, start, new MultiExtremumSettings())); - } - - /// - /// Finds a local maximum of a multi-dimensional function in the vicinity of the given starting location, - /// subject to the given evaluation constraints. - /// - /// The multi-dimensional function to maximize. - /// The starting location for the search. - /// The evaluation settings that govern the search for the maximum. - /// The local maximum. - public static MultiExtremum FindLocalMaximum (Func, double> function, IReadOnlyList start, MultiExtremumSettings settings) { - if (function == null) throw new ArgumentNullException(nameof(function)); - if (start == null) throw new ArgumentNullException(nameof(start)); - if (settings == null) throw new ArgumentNullException(nameof(settings)); - return (FindLocalExtremum(function, start, settings, true)); - } - - private static MultiExtremum FindLocalExtremum (Func, double> function, IReadOnlyList start, MultiExtremumSettings settings, bool negate) { - MultiFunctor f = new MultiFunctor(function, negate); - - // Pick an initial trust region radius; we need to do this better. - double s = 0.0; - foreach (double x in start) s += (Math.Abs(x) + 1.0 / 4.0) / 4.0; - s = s / start.Count; - Debug.WriteLine("s={0}", s); - - SetDefaultOptimizationSettings(settings, start.Count); - - return (FindMinimum_ModelTrust(f, start, s, settings)); - } - - // This method is due to Powell (http://en.wikipedia.org/wiki/Michael_J._D._Powell), but it is not what - // is usually called Powell's Method (http://en.wikipedia.org/wiki/Powell%27s_method); Powell - // developed that method in the 1960s, it was included in Numerical Recipes and is very popular. - // This is a model trust algorithm developed by Powell in the 2000s. It typically uses many - // fewer function evaluations, but does more intensive calculations between each evaluation. - - // This is basically the UOBYQA variant of Powell's new methods. It maintains a quadratic model - // that interpolates between (d + 1) (d + 2) / 2 points. The model is trusted - // within a given radius. At each step, it moves to the minimum of the model (or the boundary of - // the trust region in that direction) and evaluates the function. The new value is incorporated - // into the model and the trust region expanded or contracted depending on how accurate its - // prediction of the function value was. - - // Papers on these methods are collected at http://mat.uc.pt/~zhang/software.html#powell_software. - // The UOBYQA paper is here: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.28.1756. - // The NEWUOA paper is here: http://www.damtp.cam.ac.uk/user/na/NA_papers/NA2004_08.pdf. - // The CONDOR system (http://www.applied-mathematics.net/optimization/CONDORdownload.html) is based on these same ideas. - // The thesis of CONDOR's author (http://www.applied-mathematics.net/mythesis/index.html) was also helpful. - - // It should be very easy to extend this method to constrained optimization, either by incorporating the bounds into - // the step limits or by mapping hyper-space into a hyper-cube. - - private static MultiExtremum FindMinimum_ModelTrust (MultiFunctor f, IReadOnlyList x, double s, MultiExtremumSettings settings) { - - // Construct an initial model. - QuadraticInterpolationModel model = QuadraticInterpolationModel.Construct(f, x, s); - double trustRadius = s; - - while (f.EvaluationCount < settings.EvaluationBudget) { - - // Find the minimum point of the model within the trust radius - double[] z = model.FindMinimum(trustRadius); - double expectedValue = model.Evaluate(z); - - double deltaExpected = model.MinimumValue - expectedValue; - - // Evaluate the function at the suggested minimum - double[] point = model.ConvertPoint(z); - double value = f.Evaluate(point); - - double delta = model.MinimumValue - value; - double tol = settings.ComputePrecision(Math.Min(model.MinimumValue, value)); - // Note value can be way off, so use better of old best and new value to compute tol. - // When we didn't do this before, we got value = infinity, so tol = infinity, and thus terminated! - - if (delta > 0.0 && settings.Listener != null) { - MultiExtremum report = new MultiExtremum(f.EvaluationCount, settings, point, value, Math.Max(Math.Abs(delta), 0.75 * tol), model.GetHessian()); - settings.Listener(report); - } - - // To terminate, we demand: a reduction, that the reduction be small, that the reduction be in line with - // its expected value, that we have not run up against the trust boundary, and that the gradient is small. - // I had wanted to demand delta > 0, but we run into some cases where delta keeps being very slightly - // negative, typically orders of magnitude less than tol, causing the trust radius to shrink in an - // endless cycle that causes our approximation to ultimately go sour, even though terminating on the original - // very slightly negative delta would have produced an accurate estimate. So we tolerate this case for now. - if ((delta <= tol) && (-0.25 * tol <= delta)) { - // We demand that the model be decent, i.e. that the expected delta was within tol of the measured delta. - if (Math.Abs(delta - deltaExpected) <= tol) { - // We demand that the step not just be small because it ran up against the trust radius. - // If it ran up against the trust radius, there is probably more to be hand by continuing. - double zm = Blas1.dNrm2(z, 0, 1, z.Length); - if (zm < trustRadius) { - // Finally, we demand that the gradient be small. You might think this was obvious since - // z was small, but if the Hessian is not positive definite - // the interplay of the Hessian and the gradient can produce a small z even if the model looks nothing like a quadratic minimum. - double gm = Blas1.dNrm2(model.GetGradient(), 0, 1, z.Length); - if (gm * zm <= tol) { - if (f.IsNegated) value = -value; - return (new MultiExtremum(f.EvaluationCount, settings, point, value, Math.Max(Math.Abs(delta), 0.75 * tol), model.GetHessian())); - } - } - } - } - - // There are now three decisions to be made: - // 1. How to change the trust radius - // 2. Whether to accept the new point - // 3. Which existing point to replace - - // If the actual change was very far from the expected change, reduce the trust radius. - // If the expected change did a good job of predicting the actual change, increase the trust radius. - if ((delta < 0.25 * deltaExpected) /*|| (8.0 * deltaExpected < delta)*/) { - trustRadius = 0.5 * trustRadius; - } else if ((0.75 * deltaExpected <= delta) /*&& (delta <= 2.0 * deltaExpected)*/) { - trustRadius = 2.0 * trustRadius; - } - // It appears that the limits on delta being too large don't help, and even hurt if made too stringent. - - // Replace an old point with the new point. - int iMax = 0; double fMax = model.values[0]; - int iBad = 0; double fBad = model.ComputeBadness(0, z, point, value); - for (int i = 1; i < model.values.Length; i++) { - if (model.values[i] > fMax) { iMax = i; fMax = model.values[i]; } - double bad = model.ComputeBadness(i, z, point, value); - if (bad > fBad) { iBad = i; fBad = bad; } - } - // Use the new point as long as it is better than our worst existing point. - if (value < fMax) { - Debug.Assert(!Double.IsPositiveInfinity(value) && !Double.IsNaN(value)); - model.ReplacePoint(iBad, point, z, value); - } - // There is some question about how best to choose which point to replace. - // The largest value? The furthest away? The one closest to new min? - - } - - throw new NonconvergenceException(); - } - - } - - // Represents a quadratic multinomial, i.e. all terms up to squares. - // E.g. in two dimensions 1, x, y, x^2, xy, and y^2. - - internal class QuadraticModel { - - public QuadraticModel (int d) { - g = new double[d]; - h = new double[d][]; - for (int i = 0; i < d; i++) h[i] = new double[i + 1]; - } - - - // The constant coefficient. - public double f; - - // Linear coefficients. - public double[] g; - - // Quadratic coefficients. - public double[][] h; - - public void ShiftOrigin (double[] z) { - for (int i = 0; i < g.Length; i++) { - f += g[i] * z[i]; - for (int j = 0; j <= i; j++) { - f += z[i] * h[i][j] * z[j]; - g[i] += h[i][j] * z[j]; - } - for (int j = i; j < g.Length; j++) { - g[i] += h[j][i] * z[j]; - } - } - } - - public double Evaluate (double[] z) { - double f = this.f; - for (int i = 0; i < g.Length; i++) { - f += g[i] * z[i]; - for (int j = 0; j <= i; j++) { - f += h[i][j] * z[i] * z[j]; - } - } - return (f); - } - - // Finds the minimum of the quadratic multinomial. - public double[] FindMinimum (double trustRadius) { - - // Our model is y = f + g^T z + z^T H z / 2. Differentiating and setting equal to zero - // gives g + H z = 0 or z = H^{-1} (-g). We could solve this linear system by direct - // inversion, but that could take us to a maximum or H could be non-invertable. - - // Instead we use the truncated conjugate gradient method. - - // The location starts at the origin. - double[] z = new double[g.Length]; - - // The step direction starts as the negative gradient. - double[] u = new double[g.Length]; - for (int i = 0; i < g.Length; i++) u[i] = -g[i]; - - // We will need a vector to hold the deficit r = (-g) - H z. - double[] r = new double[g.Length]; - - // At first we need only its magnitude, and since we start with z = 0, this is just |g|. - double old_r2 = 0.0; - for (int i = 0; i < g.Length; i++) { - old_r2 += MoreMath.Sqr(g[i]); - } - if (old_r2 == 0.0) return (z); - - // Since we have a perfect quadratic form, we need to loop at most d times - // to obtain an exact solution. - for (int k = 0; k < g.Length; k++) { - - // We will step z -> z + \alpha u. - - // First compute the maximum \alpha that will keep us within the trust radius, - // i.e. \alpha s.t. |z + \alpha u| = \Delta. - double z2 = 0.0, u2 = 0.0, zu = 0.0; - for (int i = 0; i < z.Length; i++) { - z2 += MoreMath.Sqr(z[i]); - zu += z[i] * u[i]; - u2 += MoreMath.Sqr(u[i]); - } - double t = trustRadius * trustRadius - z2; - Debug.Assert(t >= 0.0); - double alphaMax = t / (Math.Sqrt(zu * zu + u2 * t) + zu); - - // Now compute \alpha as prescribed by the conjugate gradient method. - double uhu = 0.0; - for (int i = 0; i < u.Length; i++) { - for (int j = 0; j <= i; j++) { - uhu += 2.0 * u[i] * h[i][j] * u[j]; - } - } - - // If we have negative curvature or alpha exceed the maximum, use the maximum. - // Otherwise use the computed value. - double alpha; - if (uhu > 0.0) { - alpha = old_r2 / uhu; - if (alpha > alphaMax) alpha = alphaMax; - } else { - alpha = alphaMax; - } - - // Take the step. - for (int i = 0; i < z.Length; i++) z[i] += alpha * u[i]; - - // If we stepped to the boundary, terminate. - if (alpha == alphaMax) return (z); - - // If the step was very small, terminate? - - // Compute the new deficit. - double r2 = ComputeDeficit(z, ref r); - - // If we don't terminate here, in the next round - // u = 0 and \delta z = 0/0 = NaN. - // Can we terminate on some small value test? - if (r2 == 0.0) return (z); - - // Compute the next direction. - double beta = r2 / old_r2; - for (int i = 0; i < u.Length; i++) { - u[i] = r[i] + beta * u[i]; - } - - old_r2 = r2; - - } - - return (z); - - } - - private double ComputeDeficit (double[] z, ref double[] r) { - double r2 = 0.0; - for (int i = 0; i < r.Length; i++) { - r[i] = -g[i]; - for (int j = 0; j <= i; j++) { - r[i] -= h[i][j] * z[j]; - } - for (int j = i; j < z.Length; j++) { - r[i] -= h[j][i] * z[j]; - } - r2 += MoreMath.Sqr(r[i]); - } - return (r2); - } - - public void DivideBy (double q) { - f /= q; - for (int i = 0; i < g.Length; i++) { - g[i] /= q; - for (int j = 0; j <= i; j++) { - h[i][j] /= q; - } - } - } - - public void Subtract (double p, QuadraticModel Q) { - f -= p * Q.f; - for (int i = 0; i < g.Length; i++) { - g[i] -= p * Q.g[i]; - for (int j = 0; j <= i; j++) { - h[i][j] -= p * Q.h[i][j]; - } - } - } - - public double[] GetGradient () { - return (g); - } - - public double[][] GetHessian () { - double[][] H = new double[g.Length][]; - for (int i = 0; i < g.Length; i++) { - H[i] = new double[i + 1]; - for (int j = 0; j < i; j++) { - H[i][j] = h[i][j]; - } - H[i][i] = 2.0 * h[i][i]; - } - return (H); - } - - } - - // Represents a quadratic multinomial interpolation on a set - // of given points. - - internal class QuadraticInterpolationModel { - - private int d; - - public double[] origin; - -#if PAST - // delete these when finished - private double f0; - private double[] g; - private double[][] h; -#endif - - // The Lagrange polynomials for each interpolation point. - - private QuadraticModel[] polynomials; - - // The total interpolating polynomial. - - private QuadraticModel total; - - // The interpolation points and the function values at those points. - - private double[][] points; - - public double[] values; - - // We keep track of a few - - private int minValueIndex; - - //private int maxBadnessIndex; - - //public double[] badnesses; - - public double[] GetGradient () { - return (total.GetGradient()); - } - - public double[][] GetHessian () { - return (total.GetHessian()); - } - - public double MinimumValue { - get { - return (values[minValueIndex]); - } - } - - private void Initialize (MultiFunctor f, IReadOnlyList x, double s) { - - // Allocate storage - d = x.Count; - int m = (d + 1) * (d + 2) / 2; - origin = new double[d]; - points = new double[m][]; - for (int i = 0; i < m; i++) points[i] = new double[d]; - values = new double[m]; - polynomials = new QuadraticModel[m]; - for (int i = 0; i < m; i++) polynomials[i] = new QuadraticModel(d); - - // Start with x as the origin. - x.CopyTo(origin, 0); - - // The first interpolation point is the origin. - x.CopyTo(points[0], 0); - values[0] = f.Evaluate(points[0]); - - // Compute 2d more interpolation points one step along each axis. - int k = 0; - for (int i = 0; i < d; i++) { - k++; - x.CopyTo(points[k], 0); - points[k][i] += s; - double plusValue = f.Evaluate(points[k]); - values[k] = plusValue; - k++; - x.CopyTo(points[k], 0); - points[k][i] -= s; - double minusValue = f.Evaluate(points[k]); - values[k] = minusValue; - } - - // Compute d(d+1)/2 more interpolation points at the corners. - for (int i = 0; i < d; i++) { - for (int j = 0; j < i; j++) { - k++; - x.CopyTo(points[k], 0); - points[k][i] += s; - points[k][j] += s; - double cornerValue = f.Evaluate(points[k]); - values[k] = cornerValue; - } - } - - double s1 = 1.0 / s; - double s2 = s1 * s1; - - // Compute the Lagrange polynomial for each point - - k = 0; - - for (int i = 0; i < d; i++) { - k++; - polynomials[2 * i + 1].g[i] = 0.5 * s1; - polynomials[2 * i + 1].h[i][i] = 0.5 * s2; - k++; - polynomials[2 * i + 2].g[i] = -0.5 * s1; - polynomials[2 * i + 2].h[i][i] = 0.5 * s2; - } - - for (int i = 0; i < d; i++) { - for (int j = 0; j < i; j++) { - k++; - polynomials[k].h[i][j] = s2; - polynomials[2 * i + 1].h[i][j] -= s2; - polynomials[2 * j + 1].h[i][j] -= s2; - } - } - - polynomials[0].f = 1.0; - for (int l = 1; l < m; l++) { - for (int i = 0; i < d; i++) { - polynomials[0].g[i] -= polynomials[l].g[i]; - for (int j = 0; j <= i; j++) { - polynomials[0].h[i][j] -= polynomials[l].h[i][j]; - } - } - } - - // Compute the total interpolating polynomial. - - total = new QuadraticModel(d); - for (int l = 0; l < m; l++) { - total.f += values[l] * polynomials[l].f; - for (int i = 0; i < d; i++) { - total.g[i] += values[l] * polynomials[l].g[i]; - for (int j = 0; j <= i; j++) { - total.h[i][j] += values[l] * polynomials[l].h[i][j]; - } - } - } - - // Find the minimum point. - - minValueIndex = 0; - for (int i = 1; i < m; i++) { - if (values[i] < values[minValueIndex]) minValueIndex = i; - } - - // move the origin to the minimum point - double[] z = new double[d]; - for (int i = 0; i < z.Length; i++) z[i] = points[minValueIndex][i] - origin[i]; - ShiftOrigin(z); - - // compute badnesses - //badnesses = new double[points.Length]; - //for (int i = 0; i < points.Length; i++) { - // badnesses[i] = ComputeBadressAbsolute(points[i], values[i]); - //} - - } - -#if PAST - public static QuadraticInterpolationModel Construct (double[][] points, double[] values) { - - int m = points.Length; - int d = points[0].Length; - - // find the minimum point, use it as the origin - int iMin = 0; double fMin = values[0]; - for (int i = 1; i < values.Length; i++) { - if (values[i] < fMin) { iMin = i; fMin = values[i]; } - } - - SquareMatrix A = new SquareMatrix(m); - int c = 0; - for (int r = 0; r < m; r++) { - A[r, 0] = 1.0; - } - for (int i = 0; i < d; i++) { - c++; - for (int r = 0; r < m; r++) { - A[r, c] = points[r][i] - points[iMin][i]; - } - } - for (int i = 0; i < d; i++) { - for (int j = 0; j <= i; j++) { - c++; - for (int r = 0; r < m; r++) { - A[r, c] = (points[r][i] - points[iMin][i]) * (points[r][j] - points[iMin][j]); - } - } - } - ColumnVector b = new ColumnVector(values); - - SquareQRDecomposition QR = A.QRDecomposition(); - ColumnVector a = QR.Solve(b); - - QuadraticInterpolationModel model = new QuadraticInterpolationModel(); - model.d = d; - model.origin = points[iMin]; - model.f0 = a[0]; - model.g = new double[d]; - c = 0; - for (int i = 0; i < d; i++) { - c++; - model.g[i] = a[c]; - } - model.h = new double[d][]; - for (int i = 0; i < d; i++) { - model.h[i] = new double[d]; - } - for (int i = 0; i < d; i++) { - for (int j = 0; j <= i; j++) { - c++; - if (i == j) { - model.h[i][j] = 2.0 * a[c]; - } else { - model.h[i][j] = a[c]; - model.h[j][i] = a[c]; - } - } - } - - return (model); - - } - - public static QuadraticInterpolationModel Construct (Func, double> f, double[] x, double s) { - MultiFunctor mf = new MultiFunctor(f); - return (Construct(mf, x, s)); - } -#endif - - internal static QuadraticInterpolationModel Construct (MultiFunctor f, IReadOnlyList x, double s) { - QuadraticInterpolationModel model = new QuadraticInterpolationModel(); - model.Initialize(f, x, s); - return (model); - } - - // Shift to x0' = x0 + z, i.e. measure all x from point x0 + z instead of x0. - // By multiplying out quadratic form f0 + g^T (x - x0) + (x - x0)^T h (x - x0) with and without primes and demanding equality, we find - // f0' = f0 + g^T z + 1/2 z^T h z - // g' = g + h z^T - // h' = h - - public void ShiftOrigin (double[] z) { - for (int i = 0; i < origin.Length; i++) origin[i] += z[i]; - for (int i = 0; i < polynomials.Length; i++) polynomials[i].ShiftOrigin(z); - total.ShiftOrigin(z); - - } - - public double[] FindMinimum (double trustRadius) { - return (total.FindMinimum(trustRadius)); - } - - public void ReplacePoint (int index, double[] point, double[] z, double value) { - - // There are ~ m Lagrange polynomials, changing each requires ~ d^2 work, so - // the total work to replace a point is ~ m * d^2. - - // For m ~ d^2, that is ~ d^4 work total. That is a lot, but still a lot - // less than the direct inversion of a d^2 X d^2 design matrix, which would - // be ~ d^6. - - // Rescale the polynomial being replaced so that it is 1 at the new point. - // It is still zero at all the other points. - double a = polynomials[index].Evaluate(z); - polynomials[index].DivideBy(a); - - // Subtract from each of the other polynomials its value at the new point times the new polynomial. - // The result is zero at the new point, and does not affect the values at the other points, - // because the new polynomial is zero at those points. - for (int i = 0; i < polynomials.Length; i++) { - if (i == index) continue; - double b = polynomials[i].Evaluate(z); - polynomials[i].Subtract(b, polynomials[index]); - } - - // Adjust the total polynomial to reflect the changes in all the Lagrange polynomials. - double c = total.Evaluate(z) - value; - total.Subtract(c, polynomials[index]); - - // Update the minimum. - if (value < values[minValueIndex]) { - minValueIndex = index; - ShiftOrigin(z); - } - - points[index] = point; - values[index] = value; - - //badnesses[index] = ComputeBadressAbsolute(points[index], value); - } - - - public double Evaluate (double[] z) { - return (total.Evaluate(z)); - } - - public double[] ConvertPoint (double[] z) { - double[] x = new double[d]; - for (int i = 0; i < d; i++) x[i] = origin[i] + z[i]; - return (x); - } - - // Badness is used as a criteria for which point to throw out of collection. - // Basic idea is to compute a number which (i) increases with value above - // minimum and also (ii) increases with distance from minimum. - - public double ComputeBadness (int index, double[] z, double[] point, double value) { - double s = 0.0; - if (value < MinimumValue) { - // If new candidate is best, compute badness relative to it. - for (int i = 0; i < point.Length; i++) { - s += MoreMath.Sqr(points[index][i] - point[i]); - } - s = Math.Pow(s, 3.0 / 2.0) * (values[index] - value); - } else { - // Otherwise compute badness relative to best existing point. - if (index == minValueIndex) return (0.0); - for (int i = 0; i < point.Length; i++) { - s += MoreMath.Sqr(points[index][i] - origin[i]); - } - s = Math.Pow(s, 3.0 / 2.0) * (values[index] - values[minValueIndex]); - } - Debug.Assert(s >= 0.0); - return s; - } - - public double FindMinimumSeperation () { - double min = Double.MaxValue; - for (int i = 0; i < points.Length; i++) { - for (int j = 0; j < i; j++) { - double s = 0.0; - for (int k = 0; k < d; k++) { - s += MoreMath.Sqr(points[i][k] - points[j][k]); - } - s = Math.Sqrt(s); - if (s < min) min = s; - } - } - return (min); - } - - } - -} +using System; +using System.Collections.Generic; +using System.Diagnostics; + +using Meta.Numerics.Functions; +using Meta.Numerics.Matrices; + +namespace Meta.Numerics.Analysis { + + + + public static partial class MultiFunctionMath { + + /// + /// Finds a local minimum of a multi-dimensional function in the vicinity of the given starting location. + /// + /// The multi-dimensional function to minimize. + /// The starting location for the search. + /// The local minimum. + /// + /// For two dimensions, the default evaluation settings target a precision of about 10-8 + /// and allow up to about 1000 function evaluations. For increasing dimensions, these requirements are + /// gradually relaxed. + /// + /// , is . + /// The number of function evaluations required exceeded the evaluation budget. + public static MultiExtremum FindLocalMinimum (Func, double> function, IReadOnlyList start) { + return (FindLocalMinimum(function, start, new MultiExtremumSettings())); + } + + private static void SetDefaultOptimizationSettings (MultiExtremumSettings settings, int d) { + if (settings.RelativePrecision < 0.0) settings.RelativePrecision = Math.Pow(10.0, -(10.0 + 4.0 / d)); + if (settings.AbsolutePrecision < 0.0) settings.AbsolutePrecision = Math.Pow(10.0, -(10.0 + 4.0 / d)); + if (settings.EvaluationBudget < 0) settings.EvaluationBudget = 16 * (d + 1) * (d + 2) * (d + 3); + } + + /// + /// Finds a local minimum of a multi-dimensional function in the vicinity of the given starting location, + /// subject to the given evaluation constraints. + /// + /// The multi-dimensional function to minimize. + /// The starting location for the search. + /// The evaluation settings that govern the search for the minimum. + /// The local minimum. + /// + /// The Hessian (matrix of second derivatives) returned with the minimum is an approximation that is constructed in the course of search. + /// It should be considered a crude approximation, and may not even be that if the minimum is highly non-quadratic. + /// If you have a constrained minimization problem, require a high-precision solution, and do not have a good initial guess, consider first feeding + /// your constrained problem into , + /// which supports constraints but gives relatively lower precision solutions, then + /// feeding the result of that method into this method, which finds relatively precise solutions but does not support constraints. + /// + /// , , or is . + /// The number of function evaluations required exceeded the evaluation budget. + public static MultiExtremum FindLocalMinimum (Func, double> function, IReadOnlyList start, MultiExtremumSettings settings) { + if (function == null) throw new ArgumentNullException(nameof(function)); + if (start == null) throw new ArgumentNullException(nameof(start)); + if (settings == null) throw new ArgumentNullException(nameof(settings)); + return (FindLocalExtremum(function, start, settings, false)); + } + + /// + /// Finds a local maximum of a multi-dimensional function in the vicinity of the given starting location. + /// + /// The multi-dimensional function to maximize. + /// The starting location for the search. + /// The local maximum. + /// , is . + /// The number of function evaluations required exceeded the evaluation budget. + public static MultiExtremum FindLocalMaximum (Func, double> function, IReadOnlyList start) { + return (FindLocalMaximum(function, start, new MultiExtremumSettings())); + } + + /// + /// Finds a local maximum of a multi-dimensional function in the vicinity of the given starting location, + /// subject to the given evaluation constraints. + /// + /// The multi-dimensional function to maximize. + /// The starting location for the search. + /// The evaluation settings that govern the search for the maximum. + /// The local maximum. + public static MultiExtremum FindLocalMaximum (Func, double> function, IReadOnlyList start, MultiExtremumSettings settings) { + if (function == null) throw new ArgumentNullException(nameof(function)); + if (start == null) throw new ArgumentNullException(nameof(start)); + if (settings == null) throw new ArgumentNullException(nameof(settings)); + return (FindLocalExtremum(function, start, settings, true)); + } + + private static MultiExtremum FindLocalExtremum (Func, double> function, IReadOnlyList start, MultiExtremumSettings settings, bool negate) { + MultiFunctor f = new MultiFunctor(function, negate); + + // Pick an initial trust region radius; we need to do this better. + double s = 0.0; + foreach (double x in start) s += (Math.Abs(x) + 1.0 / 4.0) / 4.0; + s = s / start.Count; + Debug.WriteLine("s={0}", s); + + SetDefaultOptimizationSettings(settings, start.Count); + + return (FindMinimum_ModelTrust(f, start, s, settings)); + } + + // This method is due to Powell (http://en.wikipedia.org/wiki/Michael_J._D._Powell), but it is not what + // is usually called Powell's Method (http://en.wikipedia.org/wiki/Powell%27s_method); Powell + // developed that method in the 1960s, it was included in Numerical Recipes and is very popular. + // This is a model trust algorithm developed by Powell in the 2000s. It typically uses many + // fewer function evaluations, but does more intensive calculations between each evaluation. + + // This is basically the UOBYQA variant of Powell's new methods. It maintains a quadratic model + // that interpolates between (d + 1) (d + 2) / 2 points. The model is trusted + // within a given radius. At each step, it moves to the minimum of the model (or the boundary of + // the trust region in that direction) and evaluates the function. The new value is incorporated + // into the model and the trust region expanded or contracted depending on how accurate its + // prediction of the function value was. + + // Papers on these methods are collected at http://mat.uc.pt/~zhang/software.html#powell_software. + // The UOBYQA paper is here: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.28.1756. + // The NEWUOA paper is here: http://www.damtp.cam.ac.uk/user/na/NA_papers/NA2004_08.pdf. + // The CONDOR system (http://www.applied-mathematics.net/optimization/CONDORdownload.html) is based on these same ideas. + // The thesis of CONDOR's author (http://www.applied-mathematics.net/mythesis/index.html) was also helpful. + + // It should be very easy to extend this method to constrained optimization, either by incorporating the bounds into + // the step limits or by mapping hyper-space into a hyper-cube. + + private static MultiExtremum FindMinimum_ModelTrust (MultiFunctor f, IReadOnlyList x, double s, MultiExtremumSettings settings) { + + // Construct an initial model. + QuadraticInterpolationModel model = QuadraticInterpolationModel.Construct(f, x, s); + double trustRadius = s; + + while (f.EvaluationCount < settings.EvaluationBudget) { + + // Find the minimum point of the model within the trust radius + double[] z = model.FindMinimum(trustRadius); + double expectedValue = model.Evaluate(z); + + double deltaExpected = model.MinimumValue - expectedValue; + + // Evaluate the function at the suggested minimum + double[] point = model.ConvertPoint(z); + double value = f.Evaluate(point); + + double delta = model.MinimumValue - value; + double tol = settings.ComputePrecision(Math.Min(model.MinimumValue, value)); + // Note value can be way off, so use better of old best and new value to compute tol. + // When we didn't do this before, we got value = infinity, so tol = infinity, and thus terminated! + + if (delta > 0.0 && settings.Listener != null) { + MultiExtremum report = new MultiExtremum(f.EvaluationCount, settings, point, value, Math.Max(Math.Abs(delta), 0.75 * tol), model.GetHessian()); + settings.Listener(report); + } + + // To terminate, we demand: a reduction, that the reduction be small, that the reduction be in line with + // its expected value, that we have not run up against the trust boundary, and that the gradient is small. + // I had wanted to demand delta > 0, but we run into some cases where delta keeps being very slightly + // negative, typically orders of magnitude less than tol, causing the trust radius to shrink in an + // endless cycle that causes our approximation to ultimately go sour, even though terminating on the original + // very slightly negative delta would have produced an accurate estimate. So we tolerate this case for now. + if ((delta <= tol) && (-0.25 * tol <= delta)) { + // We demand that the model be decent, i.e. that the expected delta was within tol of the measured delta. + if (Math.Abs(delta - deltaExpected) <= tol) { + // We demand that the step not just be small because it ran up against the trust radius. + // If it ran up against the trust radius, there is probably more to be hand by continuing. + double zm = Blas1.dNrm2(z, 0, 1, z.Length); + if (zm < trustRadius) { + // Finally, we demand that the gradient be small. You might think this was obvious since + // z was small, but if the Hessian is not positive definite + // the interplay of the Hessian and the gradient can produce a small z even if the model looks nothing like a quadratic minimum. + double gm = Blas1.dNrm2(model.GetGradient(), 0, 1, z.Length); + if (gm * zm <= tol) { + if (f.IsNegated) value = -value; + return (new MultiExtremum(f.EvaluationCount, settings, point, value, Math.Max(Math.Abs(delta), 0.75 * tol), model.GetHessian())); + } + } + } + } + + // There are now three decisions to be made: + // 1. How to change the trust radius + // 2. Whether to accept the new point + // 3. Which existing point to replace + + // If the actual change was very far from the expected change, reduce the trust radius. + // If the expected change did a good job of predicting the actual change, increase the trust radius. + if ((delta < 0.25 * deltaExpected) /*|| (8.0 * deltaExpected < delta)*/) { + trustRadius = 0.5 * trustRadius; + } else if ((0.75 * deltaExpected <= delta) /*&& (delta <= 2.0 * deltaExpected)*/) { + trustRadius = 2.0 * trustRadius; + } + // It appears that the limits on delta being too large don't help, and even hurt if made too stringent. + + // Replace an old point with the new point. + int iMax = 0; double fMax = model.values[0]; + int iBad = 0; double fBad = model.ComputeBadness(0, z, point, value); + for (int i = 1; i < model.values.Length; i++) { + if (model.values[i] > fMax) { iMax = i; fMax = model.values[i]; } + double bad = model.ComputeBadness(i, z, point, value); + if (bad > fBad) { iBad = i; fBad = bad; } + } + // Use the new point as long as it is better than our worst existing point. + if (value < fMax) { + Debug.Assert(!Double.IsPositiveInfinity(value) && !Double.IsNaN(value)); + model.ReplacePoint(iBad, point, z, value); + } + // There is some question about how best to choose which point to replace. + // The largest value? The furthest away? The one closest to new min? + + } + + throw new NonconvergenceException(); + } + + } + + // Represents a quadratic multinomial, i.e. all terms up to squares. + // E.g. in two dimensions 1, x, y, x^2, xy, and y^2. + + internal class QuadraticModel { + + public QuadraticModel (int d) { + g = new double[d]; + h = new double[d][]; + for (int i = 0; i < d; i++) h[i] = new double[i + 1]; + } + + + // The constant coefficient. + public double f; + + // Linear coefficients. + public double[] g; + + // Quadratic coefficients. + public double[][] h; + + public void ShiftOrigin (double[] z) { + for (int i = 0; i < g.Length; i++) { + f += g[i] * z[i]; + for (int j = 0; j <= i; j++) { + f += z[i] * h[i][j] * z[j]; + g[i] += h[i][j] * z[j]; + } + for (int j = i; j < g.Length; j++) { + g[i] += h[j][i] * z[j]; + } + } + } + + public double Evaluate (double[] z) { + double f = this.f; + for (int i = 0; i < g.Length; i++) { + f += g[i] * z[i]; + for (int j = 0; j <= i; j++) { + f += h[i][j] * z[i] * z[j]; + } + } + return (f); + } + + // Finds the minimum of the quadratic multinomial. + public double[] FindMinimum (double trustRadius) { + + // Our model is y = f + g^T z + z^T H z / 2. Differentiating and setting equal to zero + // gives g + H z = 0 or z = H^{-1} (-g). We could solve this linear system by direct + // inversion, but that could take us to a maximum or H could be non-invertable. + + // Instead we use the truncated conjugate gradient method. + + // The location starts at the origin. + double[] z = new double[g.Length]; + + // The step direction starts as the negative gradient. + double[] u = new double[g.Length]; + for (int i = 0; i < g.Length; i++) u[i] = -g[i]; + + // We will need a vector to hold the deficit r = (-g) - H z. + double[] r = new double[g.Length]; + + // At first we need only its magnitude, and since we start with z = 0, this is just |g|. + double old_r2 = 0.0; + for (int i = 0; i < g.Length; i++) { + old_r2 += MoreMath.Sqr(g[i]); + } + if (old_r2 == 0.0) return (z); + + // Since we have a perfect quadratic form, we need to loop at most d times + // to obtain an exact solution. + for (int k = 0; k < g.Length; k++) { + + // We will step z -> z + \alpha u. + + // First compute the maximum \alpha that will keep us within the trust radius, + // i.e. \alpha s.t. |z + \alpha u| = \Delta. + double z2 = 0.0, u2 = 0.0, zu = 0.0; + for (int i = 0; i < z.Length; i++) { + z2 += MoreMath.Sqr(z[i]); + zu += z[i] * u[i]; + u2 += MoreMath.Sqr(u[i]); + } + double t = trustRadius * trustRadius - z2; + Debug.Assert(t >= 0.0); + double alphaMax = t / (Math.Sqrt(zu * zu + u2 * t) + zu); + + // Now compute \alpha as prescribed by the conjugate gradient method. + double uhu = 0.0; + for (int i = 0; i < u.Length; i++) { + for (int j = 0; j <= i; j++) { + uhu += 2.0 * u[i] * h[i][j] * u[j]; + } + } + + // If we have negative curvature or alpha exceed the maximum, use the maximum. + // Otherwise use the computed value. + double alpha; + if (uhu > 0.0) { + alpha = old_r2 / uhu; + if (alpha > alphaMax) alpha = alphaMax; + } else { + alpha = alphaMax; + } + + // Take the step. + for (int i = 0; i < z.Length; i++) z[i] += alpha * u[i]; + + // If we stepped to the boundary, terminate. + if (alpha == alphaMax) return (z); + + // If the step was very small, terminate? + + // Compute the new deficit. + double r2 = ComputeDeficit(z, ref r); + + // If we don't terminate here, in the next round + // u = 0 and \delta z = 0/0 = NaN. + // Can we terminate on some small value test? + if (r2 == 0.0) return (z); + + // Compute the next direction. + double beta = r2 / old_r2; + for (int i = 0; i < u.Length; i++) { + u[i] = r[i] + beta * u[i]; + } + + old_r2 = r2; + + } + + return (z); + + } + + private double ComputeDeficit (double[] z, ref double[] r) { + double r2 = 0.0; + for (int i = 0; i < r.Length; i++) { + r[i] = -g[i]; + for (int j = 0; j <= i; j++) { + r[i] -= h[i][j] * z[j]; + } + for (int j = i; j < z.Length; j++) { + r[i] -= h[j][i] * z[j]; + } + r2 += MoreMath.Sqr(r[i]); + } + return (r2); + } + + public void DivideBy (double q) { + f /= q; + for (int i = 0; i < g.Length; i++) { + g[i] /= q; + for (int j = 0; j <= i; j++) { + h[i][j] /= q; + } + } + } + + public void Subtract (double p, QuadraticModel Q) { + f -= p * Q.f; + for (int i = 0; i < g.Length; i++) { + g[i] -= p * Q.g[i]; + for (int j = 0; j <= i; j++) { + h[i][j] -= p * Q.h[i][j]; + } + } + } + + public double[] GetGradient () { + return (g); + } + + public double[][] GetHessian () { + double[][] H = new double[g.Length][]; + for (int i = 0; i < g.Length; i++) { + H[i] = new double[i + 1]; + for (int j = 0; j < i; j++) { + H[i][j] = h[i][j]; + } + H[i][i] = 2.0 * h[i][i]; + } + return (H); + } + + } + + // Represents a quadratic multinomial interpolation on a set + // of given points. + + internal class QuadraticInterpolationModel { + + private int d; + + public double[] origin; + +#if PAST + // delete these when finished + private double f0; + private double[] g; + private double[][] h; +#endif + + // The Lagrange polynomials for each interpolation point. + + private QuadraticModel[] polynomials; + + // The total interpolating polynomial. + + private QuadraticModel total; + + // The interpolation points and the function values at those points. + + private double[][] points; + + public double[] values; + + // We keep track of a few + + private int minValueIndex; + + //private int maxBadnessIndex; + + //public double[] badnesses; + + public double[] GetGradient () { + return (total.GetGradient()); + } + + public double[][] GetHessian () { + return (total.GetHessian()); + } + + public double MinimumValue { + get { + return (values[minValueIndex]); + } + } + + private void Initialize (MultiFunctor f, IReadOnlyList x, double s) { + + // Allocate storage + d = x.Count; + int m = (d + 1) * (d + 2) / 2; + origin = new double[d]; + points = new double[m][]; + for (int i = 0; i < m; i++) points[i] = new double[d]; + values = new double[m]; + polynomials = new QuadraticModel[m]; + for (int i = 0; i < m; i++) polynomials[i] = new QuadraticModel(d); + + // Start with x as the origin. + x.CopyTo(origin, 0); + + // The first interpolation point is the origin. + x.CopyTo(points[0], 0); + values[0] = f.Evaluate(points[0]); + + // Compute 2d more interpolation points one step along each axis. + int k = 0; + for (int i = 0; i < d; i++) { + k++; + x.CopyTo(points[k], 0); + points[k][i] += s; + double plusValue = f.Evaluate(points[k]); + values[k] = plusValue; + k++; + x.CopyTo(points[k], 0); + points[k][i] -= s; + double minusValue = f.Evaluate(points[k]); + values[k] = minusValue; + } + + // Compute d(d+1)/2 more interpolation points at the corners. + for (int i = 0; i < d; i++) { + for (int j = 0; j < i; j++) { + k++; + x.CopyTo(points[k], 0); + points[k][i] += s; + points[k][j] += s; + double cornerValue = f.Evaluate(points[k]); + values[k] = cornerValue; + } + } + + double s1 = 1.0 / s; + double s2 = s1 * s1; + + // Compute the Lagrange polynomial for each point + + k = 0; + + for (int i = 0; i < d; i++) { + k++; + polynomials[2 * i + 1].g[i] = 0.5 * s1; + polynomials[2 * i + 1].h[i][i] = 0.5 * s2; + k++; + polynomials[2 * i + 2].g[i] = -0.5 * s1; + polynomials[2 * i + 2].h[i][i] = 0.5 * s2; + } + + for (int i = 0; i < d; i++) { + for (int j = 0; j < i; j++) { + k++; + polynomials[k].h[i][j] = s2; + polynomials[2 * i + 1].h[i][j] -= s2; + polynomials[2 * j + 1].h[i][j] -= s2; + } + } + + polynomials[0].f = 1.0; + for (int l = 1; l < m; l++) { + for (int i = 0; i < d; i++) { + polynomials[0].g[i] -= polynomials[l].g[i]; + for (int j = 0; j <= i; j++) { + polynomials[0].h[i][j] -= polynomials[l].h[i][j]; + } + } + } + + // Compute the total interpolating polynomial. + + total = new QuadraticModel(d); + for (int l = 0; l < m; l++) { + total.f += values[l] * polynomials[l].f; + for (int i = 0; i < d; i++) { + total.g[i] += values[l] * polynomials[l].g[i]; + for (int j = 0; j <= i; j++) { + total.h[i][j] += values[l] * polynomials[l].h[i][j]; + } + } + } + + // Find the minimum point. + + minValueIndex = 0; + for (int i = 1; i < m; i++) { + if (values[i] < values[minValueIndex]) minValueIndex = i; + } + + // move the origin to the minimum point + double[] z = new double[d]; + for (int i = 0; i < z.Length; i++) z[i] = points[minValueIndex][i] - origin[i]; + ShiftOrigin(z); + + // compute badnesses + //badnesses = new double[points.Length]; + //for (int i = 0; i < points.Length; i++) { + // badnesses[i] = ComputeBadressAbsolute(points[i], values[i]); + //} + + } + +#if PAST + public static QuadraticInterpolationModel Construct (double[][] points, double[] values) { + + int m = points.Length; + int d = points[0].Length; + + // find the minimum point, use it as the origin + int iMin = 0; double fMin = values[0]; + for (int i = 1; i < values.Length; i++) { + if (values[i] < fMin) { iMin = i; fMin = values[i]; } + } + + SquareMatrix A = new SquareMatrix(m); + int c = 0; + for (int r = 0; r < m; r++) { + A[r, 0] = 1.0; + } + for (int i = 0; i < d; i++) { + c++; + for (int r = 0; r < m; r++) { + A[r, c] = points[r][i] - points[iMin][i]; + } + } + for (int i = 0; i < d; i++) { + for (int j = 0; j <= i; j++) { + c++; + for (int r = 0; r < m; r++) { + A[r, c] = (points[r][i] - points[iMin][i]) * (points[r][j] - points[iMin][j]); + } + } + } + ColumnVector b = new ColumnVector(values); + + SquareQRDecomposition QR = A.QRDecomposition(); + ColumnVector a = QR.Solve(b); + + QuadraticInterpolationModel model = new QuadraticInterpolationModel(); + model.d = d; + model.origin = points[iMin]; + model.f0 = a[0]; + model.g = new double[d]; + c = 0; + for (int i = 0; i < d; i++) { + c++; + model.g[i] = a[c]; + } + model.h = new double[d][]; + for (int i = 0; i < d; i++) { + model.h[i] = new double[d]; + } + for (int i = 0; i < d; i++) { + for (int j = 0; j <= i; j++) { + c++; + if (i == j) { + model.h[i][j] = 2.0 * a[c]; + } else { + model.h[i][j] = a[c]; + model.h[j][i] = a[c]; + } + } + } + + return (model); + + } + + public static QuadraticInterpolationModel Construct (Func, double> f, double[] x, double s) { + MultiFunctor mf = new MultiFunctor(f); + return (Construct(mf, x, s)); + } +#endif + + internal static QuadraticInterpolationModel Construct (MultiFunctor f, IReadOnlyList x, double s) { + QuadraticInterpolationModel model = new QuadraticInterpolationModel(); + model.Initialize(f, x, s); + return (model); + } + + // Shift to x0' = x0 + z, i.e. measure all x from point x0 + z instead of x0. + // By multiplying out quadratic form f0 + g^T (x - x0) + (x - x0)^T h (x - x0) with and without primes and demanding equality, we find + // f0' = f0 + g^T z + 1/2 z^T h z + // g' = g + h z^T + // h' = h + + public void ShiftOrigin (double[] z) { + for (int i = 0; i < origin.Length; i++) origin[i] += z[i]; + for (int i = 0; i < polynomials.Length; i++) polynomials[i].ShiftOrigin(z); + total.ShiftOrigin(z); + + } + + public double[] FindMinimum (double trustRadius) { + return (total.FindMinimum(trustRadius)); + } + + public void ReplacePoint (int index, double[] point, double[] z, double value) { + + // There are ~ m Lagrange polynomials, changing each requires ~ d^2 work, so + // the total work to replace a point is ~ m * d^2. + + // For m ~ d^2, that is ~ d^4 work total. That is a lot, but still a lot + // less than the direct inversion of a d^2 X d^2 design matrix, which would + // be ~ d^6. + + // Rescale the polynomial being replaced so that it is 1 at the new point. + // It is still zero at all the other points. + double a = polynomials[index].Evaluate(z); + polynomials[index].DivideBy(a); + + // Subtract from each of the other polynomials its value at the new point times the new polynomial. + // The result is zero at the new point, and does not affect the values at the other points, + // because the new polynomial is zero at those points. + for (int i = 0; i < polynomials.Length; i++) { + if (i == index) continue; + double b = polynomials[i].Evaluate(z); + polynomials[i].Subtract(b, polynomials[index]); + } + + // Adjust the total polynomial to reflect the changes in all the Lagrange polynomials. + double c = total.Evaluate(z) - value; + total.Subtract(c, polynomials[index]); + + // Update the minimum. + if (value < values[minValueIndex]) { + minValueIndex = index; + ShiftOrigin(z); + } + + points[index] = point; + values[index] = value; + + //badnesses[index] = ComputeBadressAbsolute(points[index], value); + } + + + public double Evaluate (double[] z) { + return (total.Evaluate(z)); + } + + public double[] ConvertPoint (double[] z) { + double[] x = new double[d]; + for (int i = 0; i < d; i++) x[i] = origin[i] + z[i]; + return (x); + } + + // Badness is used as a criteria for which point to throw out of collection. + // Basic idea is to compute a number which (i) increases with value above + // minimum and also (ii) increases with distance from minimum. + + public double ComputeBadness (int index, double[] z, double[] point, double value) { + double s = 0.0; + if (value < MinimumValue) { + // If new candidate is best, compute badness relative to it. + for (int i = 0; i < point.Length; i++) { + s += MoreMath.Sqr(points[index][i] - point[i]); + } + s = Math.Pow(s, 3.0 / 2.0) * (values[index] - value); + } else { + // Otherwise compute badness relative to best existing point. + if (index == minValueIndex) return (0.0); + for (int i = 0; i < point.Length; i++) { + s += MoreMath.Sqr(points[index][i] - origin[i]); + } + s = Math.Pow(s, 3.0 / 2.0) * (values[index] - values[minValueIndex]); + } + Debug.Assert(s >= 0.0); + return s; + } + + } + +} diff --git a/Numerics/Extended/Int128.cs b/Numerics/Extended/Int128.cs index 49f19e8..24d861f 100644 --- a/Numerics/Extended/Int128.cs +++ b/Numerics/Extended/Int128.cs @@ -258,6 +258,23 @@ public static explicit operator Int128 (BigInteger b) { return new Int128(u); } + /// + /// Converts a floating-point value into a 128-bit integer. + /// + /// A floating-point number. + /// The lower 128 bits of the integer part of the floating-point number. + /// is NaN, or infinite. + 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 /// @@ -339,6 +356,28 @@ public static bool TryParse (string text, out Int128 value) { return new Int128(x.u - y.u); } + /// + /// Increments a 128-bit integer. + /// + /// The integer. + /// One more than . + public static Int128 operator ++ (Int128 x) { + UInt128 v = x.u; + v++; + return new Int128(v); + } + + /// + /// Decrements a 128-bit unsigned integer. + /// + /// The integer. + /// One less than . + public static Int128 operator -- (Int128 x) { + UInt128 v = x.u; + v--; + return new Int128(v); + } + // Multiplication /// @@ -416,6 +455,17 @@ public static Int128 DivRem (Int128 x, Int128 y, out Int128 r) { } + /// + /// Computes the remainder of two 128 bit integers. + /// + /// The dividend. + /// The divisor. + /// The remainder of divided by . + public static Int128 operator % (Int128 x, Int128 y) { + Int128.DivRem(x, y, out Int128 r); + return r; + } + // Absolute Value /// diff --git a/Numerics/Extended/Int128Calculator.cs b/Numerics/Extended/Int128Calculator.cs index 50ea83d..115736c 100644 --- a/Numerics/Extended/Int128Calculator.cs +++ b/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) { diff --git a/Numerics/Extended/UInt128.cs b/Numerics/Extended/UInt128.cs index f0c4b41..2641004 100644 --- a/Numerics/Extended/UInt128.cs +++ b/Numerics/Extended/UInt128.cs @@ -170,7 +170,7 @@ public static explicit operator UInt128 (BigInteger b) { /// The lower 128 bits of the integer part of the floating-point number. /// is negative, or NaN, or infinite. 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); } + /// + /// Increments a 128-bit unsigned integer. + /// + /// The integer. + /// One more than (or , if is ). + public static UInt128 operator ++ (UInt128 x) { + Int128Calculator.Increment128(x.s1, x.s0, out ulong s1, out ulong s0); + return new UInt128(s1, s0); + } + /// /// Subtracts one 128 bit unsigned integer from another. /// @@ -354,6 +360,16 @@ public int CompareTo (UInt128 u) { return new UInt128(s1, s0); } + /// + /// Decrements a 128-bit unsigned integer. + /// + /// The integer. + /// One less than (or , if is ). + 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 /// diff --git a/Numerics/Functions/AdvancedMath_Airy.cs b/Numerics/Functions/AdvancedMath_Airy.cs index f2171a1..dd77a6d 100644 --- a/Numerics/Functions/AdvancedMath_Airy.cs +++ b/Numerics/Functions/AdvancedMath_Airy.cs @@ -379,6 +379,7 @@ private static SolutionPair Airy_Asymptotic_Negative (double x) { /// The index of the zero. /// The th value of x for which Ai(x) = 0. /// is less than 1. + /// public static double AiryAiZero (int k) { return (BesselMath.AiryAiZero(k)); } @@ -389,6 +390,7 @@ public static double AiryAiZero (int k) { /// The index of the zero. /// The th value of x for which Bi(x) = 0. /// is less than 1. + /// public static double AiryBiZero (int k) { return (BesselMath.AiryBiZero(k)); } diff --git a/Numerics/Functions/AdvancedMath_Bessel.cs b/Numerics/Functions/AdvancedMath_Bessel.cs index 5e7c6ee..bc81402 100644 --- a/Numerics/Functions/AdvancedMath_Bessel.cs +++ b/Numerics/Functions/AdvancedMath_Bessel.cs @@ -19,6 +19,7 @@ public static partial class AdvancedMath { /// For information on the cylindrical Bessel functions, see . /// /// + /// 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) { /// For information on the cylindrical Bessel functions, see . /// /// + /// 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) { /// /// is negative. /// + /// 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) { /// is negative. /// /// + /// 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) { /// /// The order, which must be non-negative. /// The index of the zero, which must be positive. - /// The th value of x for which Jnu(x) = 0. + /// The th value of x for which Jν(x) = 0. /// is negative or is non-positive. 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 diff --git a/Numerics/Functions/AdvancedMath_Exponential.cs b/Numerics/Functions/AdvancedMath_Exponential.cs index ff46798..2ea4bd0 100644 --- a/Numerics/Functions/AdvancedMath_Exponential.cs +++ b/Numerics/Functions/AdvancedMath_Exponential.cs @@ -227,11 +227,14 @@ public static double IntegralCi (double x) { /// The argument. /// The value of Cin(x). /// - /// The entine cosine integral Cin(x) is related to the conventional cosine integral Ci(x) by: + /// The entire cosine integral Cin(x) is related to the conventional cosine integral Ci(x) by: + /// /// 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. + /// 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. /// + /// + /// /// public static double IntegralCin (double x) { if (x < 0.0) { @@ -278,6 +281,8 @@ public static double IntegralSi (double x) { /// /// The argument. /// The value of Shi(x). + /// + /// 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; kThe cylindrical angle φ. This angle is usually expressed as between 0 and 2π, measured counter-clockwise (as seen from above) from the positive x-axis. It is also possible to use negative values to represent clockwise movement. /// The value of Yl,m(θ,φ). /// is negative, or lies outside the range [-l, l]. + /// 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)); diff --git a/Numerics/Functions/AdvancedMath_Hypergeometric.cs b/Numerics/Functions/AdvancedMath_Hypergeometric.cs index 185f8d6..4dae535 100644 --- a/Numerics/Functions/AdvancedMath_Hypergeometric.cs +++ b/Numerics/Functions/AdvancedMath_Hypergeometric.cs @@ -18,7 +18,7 @@ public static partial class AdvancedMath { /// The value of 2F1(a, b; c; x). /// /// The Gauss Hypergeometric function is defined by the hypergeometric series: - /// + /// /// 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. diff --git a/Numerics/Functions/AdvancedMath_Lambert.cs b/Numerics/Functions/AdvancedMath_Lambert.cs index a784905..9f3c492 100644 --- a/Numerics/Functions/AdvancedMath_Lambert.cs +++ b/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. /// + /// is less than -1/e. /// - /// + /// + /// public static double LambertW (double x) { if (x < -EI) throw new ArgumentOutOfRangeException(nameof(x)); diff --git a/Numerics/Functions/AdvancedMath_ModifiedBessel.cs b/Numerics/Functions/AdvancedMath_ModifiedBessel.cs index 8ca3d82..2ef3ee7 100644 --- a/Numerics/Functions/AdvancedMath_ModifiedBessel.cs +++ b/Numerics/Functions/AdvancedMath_ModifiedBessel.cs @@ -1,665 +1,670 @@ -using System; -using System.Diagnostics; - -using Meta.Numerics; - -namespace Meta.Numerics.Functions { - - // mostly direct and straightforward adaptations of the Bessel routines - - public static partial class AdvancedMath { - - /// - /// Computes modified cylindrical Bessel functions. - /// - /// The order, which must be non-negative. - /// The argument, which must be non-negative. - /// The values of I, I', K, and K' for the given order and argument. - /// - /// The modified Bessel functions fulfill a differential equation similar to the Bessel differential equation. - /// - /// - public static SolutionPair ModifiedBessel (double nu, double x) { - if (nu < 0.0) throw new ArgumentOutOfRangeException(nameof(nu)); - if (x < 0.0) throw new ArgumentOutOfRangeException(nameof(x)); - - if (x == 0.0) { - return Bessel_Zero(nu, +1); - } else if (x < 2.0) { - - // Use the series to determine I and I' - ModifiedBesselI_Series(nu, x, out double I, out double IP); - - // Use the series to determine K and K' at -1/2 <= mu <= 1/2 that is an integer offset from nu - int n = (int) Math.Round(nu); - double mu = nu - n; - Debug.Assert(Math.Abs(mu) <= 0.5); - - // Determine K and K' at order mu - ModifiedBesselK_Series(mu, x, out double K, out double K1); - double KP = (mu / x) * K - K1; - - // Recurr K and K' upward to order nu - ModifiedBesselK_RecurrUpward(mu, x, ref K, ref KP, n); - - // Return the result - return (new SolutionPair(I, IP, K, KP)); - - } else if (x > 32.0 + 0.5 * nu * nu) { - ModifiedBessel_Asymptotic_Scaled(nu, x, out double sI, out double sIP, out double sK, out double sKP); - double e = Math.Exp(x); - return (new SolutionPair(e * sI, e * sIP, sK / e, sKP / e)); - - } else { - - // find 0 <= mu < 1 with same fractional part as nu - // this is necessary because CF2 does not produce K with good accuracy except at very low orders - int n = (int) Math.Floor(nu); - double mu = nu - n; - - // compute K, K' at this point (which is beyond the turning point because mu is small) using CF2 - ModifiedBessel_CF_K(mu, x, out double sK, out double g); - double sKP = g * sK; - - // recurr upward to order nu - ModifiedBesselK_RecurrUpward(mu, x, ref sK, ref sKP, n); - - // determine I'/I at the desired point - double f = ModifiedBessel_CF1(nu, x); - - // Use the Wronskian relationship K I' - I K' = 1/x to determine I and I' separately - double sI = 1.0 / (f * sK - sKP) / x; - double sIP = f * sI; - - double e = Math.Exp(x); - return (new SolutionPair(e * sI, e * sIP, sK / e, sKP / e)); - - } - - } - - // The recurrence for the modified Bessel functions (A&S 9.6.26, DLMF 10.29.2) is - // F_{\nu + 1} = F_{\nu}' - (\nu / x) F - // It is true for F = I or (-1)^{\nu} K, but since I is supressed as \nu increases, it should only be used for K upward, - // and only be used for I downward. - - // Multiplying by e^{\pm x}, we see that the exact same recurrence applies to sI or sK. - - // If you can compute K_0 and K_1, it can be started using K_0' = -K_1. - - private static void ModifiedBesselK_RecurrUpward (double mu, double x, ref double K, ref double KP, int n) { - for (int i = 0; i < n; i++) { - double t = K; - K = (mu / x) * K - KP; - mu += 1.0; - KP = -(t + (mu / x) * K); - } - } - - /// - /// Computes the regular modified cylindrical Bessel function. - /// - /// The order parameter. - /// The argument, which must be non-negative. - /// The value of Iν(x). - /// - /// The modified Bessel functions appear as the solutions of hyperbolic differential equations with - /// cylindrical or circular symmetry, for example the conduction of heat through a cylindrical pipe. - /// The regular modified Bessel functions are related to the Bessel functions with pure imaginary arguments. - /// - /// The regular modified Bessel functions increase monotonically and exponentially from the origin. - /// Because they increase exponentially, this function overflows for even moderately large arguments. In this - /// regeime, you can still obtain the value of the scaled modified e-x Iν(x) - /// by calling . - /// - /// - /// - public static double ModifiedBesselI (double nu, double x) { - return ModifiedBesselI(nu, x, false); - } - - /// - /// Computes the exponentially re-scaled regular modified cylindrical Bessel function. - /// - /// The order parameter. - /// The argument, which must be non-negative. - /// The value of e-x Iν(x). - /// - /// - public static double ScaledModifiedBesselI (double nu, double x) { - return ModifiedBesselI(nu, x, true); - } - - private static double ModifiedBesselI (double nu, double x, bool scaled) { - - if (x < 0.0) throw new ArgumentOutOfRangeException(nameof(x)); - - // Reflect negative nu to positive nu. - if (nu < 0.0) { - if (Math.Round(nu) == nu) { - return (ModifiedBesselI(-nu, x)); - } else { - return (ModifiedBesselI(-nu, x) + 2.0 / Math.PI * MoreMath.SinPi(-nu) * ModifiedBesselK(-nu, x)); - } - } - - if (x == 0.0) { - if (nu == 0.0) { - return (1.0); - } else { - return (0.0); - } - } else if (x <= 2.0 * Math.Sqrt(nu + 1.0)) { - // Close to the origin, use the power series. - double I = ModifiedBesselI_Series(nu, x); - return scaled ? Math.Exp(-x) * I : I; - } else if (x < 32.0 + 0.5 * nu * nu) { - // In the intermediate region, we will use CF1 + CF2 in a slightly different way than for the standard Bessel functions. - - // Find 0 <= mu < 1 with same fractional part as nu. - // This is necessary because CF2 does not produce K with good accuracy except at very low orders. - int n = (int) Math.Floor(nu); - double mu = nu - n; - Debug.Assert(0.0 <= mu && mu < 1.0); - - // Compute K, K' at this point (which is beyond the turning point because mu is small) using CF2. - ModifiedBessel_CF_K(mu, x, out double sK, out double g); - double sKP = g * sK; - - // Recurr upward to the desired order. - // This is okay because K increases with increasing order. - ModifiedBesselK_RecurrUpward(mu, x, ref sK, ref sKP, n); - - // Determine I'/I at the desired point. - double f = ModifiedBessel_CF1(nu, x); - - // Use the Wronskian to determine I from (I'/I), (K'/K), and K. - double sI = 1.0 / (f * sK - sKP) / x; - - return scaled ? sI : Math.Exp(x) * sI; - } else { - // Far from the origin, use the asymptotic expansion. - ModifiedBessel_Asymptotic_Scaled(nu, x, out double sI, out double sIP, out double sK, out double sKP); - return scaled ? sI : Math.Exp(x) * sI; - } - - } - - /// - /// Computes the irregular modified cylindrical Bessel function. - /// - /// The order parameter. - /// The argument, which must be non-negative. - /// The value of Kν(x). - /// - /// The modified Bessel functions are related to the Bessel functions with pure imaginary arguments. - /// The irregular modified Bessel function decreases monotonically and exponentially from the origin. - /// - public static double ModifiedBesselK (double nu, double x) { - return ModifiedBesselK(nu, x, false); - } - - /// - /// Computes the exponentially re-scaled irregular modified cylindrical Bessel function. - /// - /// The order parameter. - /// The argument, which must be non-negative. - /// The value of ex Kν(x). - public static double ScaledModifiedBesselK (double nu, double x) { - return ModifiedBesselK(nu, x, true); - } - - private static double ModifiedBesselK (double nu, double x, bool scaled) { - - if (x < 0.0) throw new ArgumentOutOfRangeException(nameof(x)); - - if (nu < 0.0) return (ModifiedBesselK(-nu, x)); - - if (x == 0.0) { - return (Double.PositiveInfinity); - } else if (x <= 2.0) { - // For small x, determine the value at small mu using the series and recurr up. - int n = (int) Math.Round(nu); - double mu = nu - n; - Debug.Assert(n >= 0); - Debug.Assert(Math.Abs(mu) <= 0.5); - - ModifiedBesselK_Series(mu, x, out double K, out double K1); - - if (n == 0) { - // K is already for order nu - } else if (n == 1) { - K = K1; - } else { - double KP = (mu / x) * K - K1; - ModifiedBesselK_RecurrUpward(mu, x, ref K, ref KP, n); - } - - return scaled ? Math.Exp(x) * K : K; - } else if (x < 32.0 + 0.5 * nu * nu) { - // In the intermediate region, determine the value at small mu using continued fraction and recurr up. - int n = (int) Math.Round(nu); - double mu = nu - n; - Debug.Assert(n >= 0); - Debug.Assert(Math.Abs(mu) <= 0.5); - - ModifiedBessel_CF_K(mu, x, out double sK, out double g); - double sKP = g * sK; - ModifiedBesselK_RecurrUpward(mu, x, ref sK, ref sKP, n); - return scaled ? sK : Math.Exp(-x) * sK; - } else { - // For large x, use the asymptotic series. - ModifiedBessel_Asymptotic_Scaled(nu, x, out double sI, out double sIP, out double sK, out double sKP); - return scaled ? sK : Math.Exp(-x) * sK; - } - - } - - // series near the origin; this is entirely analogous to the Bessel series near the origin - // it has a corresponding radius of rapid convergence, x < 4 + 2 Sqrt(nu) - - // This is exactly the same as BesselJ_Series with xx -> -xx. - // We could even factor this out into a common method with an additional parameter. - - private static void ModifiedBesselI_Series (double nu, double x, out double I, out double IP) { - Debug.Assert(x > 0.0); - - double x2 = 0.5 * x; - double xx = x2 * x2; - double dI = AdvancedMath.PowerOverFactorial(x2, nu); - I = dI; - IP = nu * dI; - for (int k = 1; k < Global.SeriesMax; k++) { - double I_old = I; - double IP_old = IP; - dI *= xx / (k * (nu + k)); - I += dI; - IP += (2 * k + nu) * dI; - if ((I == I_old) && (IP == IP_old)) { - IP = IP / x; - return; - } - } - - throw new NonconvergenceException(); - } - - private static double ModifiedBesselI_Series (double nu, double x) { - - double x2 = 0.5 * x; - double xx = x2 * x2; - double dI = AdvancedMath.PowerOverFactorial(x2, nu); - double I = dI; - for (int k = 1; k < Global.SeriesMax; k++) { - double I_old = I; - dI = dI * xx / (k * (nu + k)); - I += dI; - if (I == I_old) { - return (I); - } - } - throw new NonconvergenceException(); - - } - - // good for x > 32 + nu^2 / 2 - // this returns scaled values; I = e^(x) sI, K = e^(-x) sK, I' = e^(x) sI', K' = e^(-x) K' - // Note I'/K' != d/dx I/K under this scaling convention - - // Asymptotic expansions - // I_{\nu}(x) = \frac{e^x}{\sqrt{2 \pi x}} \left[ 1 - \frac{\mu-1}{8x} + \frac{(\mu-1)(\mu-9)}{2! (8x)^2} + \cdots \right] - // K_{\nu}(x) = \sqrt{\frac{\pi}{2x}}e^{-x} \left[ 1 + \frac{\mu-1}{8x} + \frac{(\mu-1)(\mu-9)}{2! (8x)^2} + \cdots \right] - // where \mu = 4 \nu^2. Derivatives - // I_{\nu}'(x) = \frac{e^x}{\sqrt{2 \pi x}} \left[ 1 - \frac{\mu+3}{8x} + \frac{(\mu-1)(\mu+15)}{2! (8x)^2} + \cdots \right] - // K_{\nu}'(x) = -\sqrt{\frac{\pi}{2x}}e^{-x} \left[ 1 + \frac{\mu+3}{8x} + \frac{(\mu-1)(\mu+15)}{2! (8x)^2} + \cdots \right] - // Note series differ only by alternating vs. same sign of terms. - - private static void ModifiedBessel_Asymptotic_Scaled (double nu, double x, out double sI, out double sIP, out double sK, out double sKP) { - - // precompute some values we will use - double mu = 4.0 * nu * nu; - double xx = 8.0 * x; - - // initialize the series - sI = 1.0; - sIP = 1.0; - sK = 1.0; - sKP = 1.0; - - // intialize the term value - double t = 1.0; - - for (int k = 1; k < Global.SeriesMax; k++) { - - double sI_old = sI; double sK_old = sK; - double sIP_old = sIP; double sKP_old = sKP; - - // determine next term values - int k2 = 2 * k; - t /= k * xx; - double tp = (mu + (k2 * k2 - 1)) * t; - //t *= (mu - (k2-1) * (k2-1)); - t *= (2.0 * (nu - k) + 1.0) * (2.0 * (nu + k) - 1.0); - - // add them, with alternating-sign for I series and same-sign for K series - if (k % 2 == 0) { - sI += t; - sIP += tp; - } else { - sI -= t; - sIP -= tp; - } - sK += t; - sKP += tp; - - // check for convergence - if ((sI == sI_old) && (sK == sK_old) && (sIP == sIP_old) && (sKP == sKP_old)) { - double fI = Math.Sqrt(2.0 * Math.PI * x); - double fK = Math.Sqrt(0.5 * Math.PI / x); - sI /= fI; - sIP /= fI; - sK *= fK; - sKP *= -fK; - return; - } - - } - - throw new NonconvergenceException(); - - } - - // compute I'/I via continued fraction - // this is analogous to CF1 for normal Bessel functions - // it converges quickly for x < nu and slowly for x > nu - - private static double ModifiedBessel_CF1 (double nu, double x) { - double q = 2.0 * (nu + 1.0); - double a = 1.0; - double b = (q / x); - double D = 1.0 / b; - double Df = a / b; - double f = nu / x + Df; - for (int k = 2; k < Global.SeriesMax; k++) { - double f_old = f; - - q += 2.0; - a = 1.0; - b = q / x; - D = 1.0 / (b + a * D); - Df = (b * D - 1.0) * Df; - f += Df; - - if (f == f_old) { - return (f); - } - } - - throw new NonconvergenceException(); - } - - // the continued fraction here is - - private static void ModifiedBessel_CF_K (double nu, double x, out double sK, out double g) { - - double a = 1.0; - double b = 2.0 * (1.0 + x); - - double D = 1.0 / b; - double Df = a / b; - double f = Df; - - // Thompson-Barnett do normalizing sum using these auxiliary variables (see Numerical Recipes) - - double C = (0.5 - nu) * (0.5 + nu); - double q0 = 0.0; - double q1 = 1.0; - double Q = C * q1; - double S = 1.0 + Q * Df; - - for (int k = 2; k < Global.SeriesMax; k++) { - - double f_old = f; - - // Lenz method for z1/z0 - - a = -(k - nu - 0.5) * (k + nu - 0.5); - b += 2.0; - - D = 1.0 / (b + a * D); - Df = (b * D - 1.0) * Df; - f += Df; - - // Thompson-Barnett sum trick - - double S_old = S; - - C = -a / k * C; - double q2 = (q0 - (b - 2.0) * q1) / a; - Q += C * q2; - S += Q * Df; - - if ((S == S_old) && (f == f_old)) { - g = -((x + 0.5) + (nu + 0.5) * (nu - 0.5) * f) / x; - sK = Math.Sqrt(Math.PI / 2.0 / x) / S; - return; - } - - q0 = q1; q1 = q2; - - } - - throw new NonconvergenceException(); - } - - // Temme's series for K0 and K1 for -1/2 \le \ny \le 1/2 and small x is described in NR3 6.6.37-6.6.40. - - // K_{\nu} = \sum_{k=0}^{\infty} c_k f_k, K_{\nu+1} = \frac{2}{x} \sum_{k=0}^{\infty} c_k h_k - - // Here c_k = \frac{(x/2)^{2k}}{k!} and f and h are determined by recurrences - - // p_{k} = \frac{p_{k-1}}{k - \nu} - // q_{k} = \frac{q_{k-1}}{k + \nu} - // f_{k} = \frac{k f_{k-1} + p_{k-1} + q_{k-1}}{k^2 - \nu^2} - // h_{k} = p_{k} - k f_{k} - - // with initial values - - // p_0 = \frac{1}{2} \left( \frac{x}{2} \right)^{-\nu} \Gamma(1 + \nu) - // q_0 = \frac{1}{2} \left( \frac{x}{2} \right)^{\nu} \Gamma(1 - \nu) - // f_0 = \frac{\nu\pi}{\sin(\nu\pi)} \left[ ... \right] - - // I observe it to be accurate to all but last couple digits for x < 2; it converges further out but - // its accuracy deteriorates rapidly. - - // The initial constant determination logic is shared with BesselY_Series; we should factor out the common parts. - - private static void ModifiedBesselK_Series (double nu, double x, out double K0, out double K1) { - - Debug.Assert((-0.5 <= nu) && (nu <= 0.5)); - Debug.Assert(x > 0.0); - - double x2 = 0.5 * x; - double x4 = x2 * x2; - - // determine initial p, q - - double GP = Gamma(1.0 + nu); - double GM = Gamma(1.0 - nu); - double z = Math.Pow(x2, nu); - double p = 0.5 / z * GP; - double q = 0.5 * z * GM; - - // determine initial f; this is rather complicated - - double ln = -Math.Log(x2); - - double s = nu * ln; - double cosh = Math.Cosh(s); - double sinhc = (s == 0.0) ? 1.0 : Math.Sinh(s) / s; - - double G1 = (Math.Abs(nu) < 0.25) ? NewG(1.0 - nu, 2.0 * nu) : (1.0 / GM - 1.0 / GP) / (2.0 * nu); - double G2 = (1.0 / GM + 1.0 / GP) / 2.0; - double F = (nu == 0.0) ? 1.0 : (nu * Math.PI) / MoreMath.SinPi(nu); - - double f = F * (cosh * G1 + sinhc * G2 * ln); - - // run the series - - double c = 1.0; - K0 = f; - K1 = p; - for (int k = 1; k < Global.SeriesMax; k++) { - - double K0_old = K0; - double K1_old = K1; - - c *= x4 / k; - - double kpnu = k + nu; - double kmnu = k - nu; - f = (k * f + p + q) / (kpnu * kmnu); - p = p / kmnu; - q = q / kpnu; - double h = p - k * f; - - K0 += c * f; - K1 += c * h; - - if ((K0 == K0_old) && (K1 == K1_old)) { - K1 = 2.0 / x * K1; - return; - } - } - throw new NonconvergenceException(); - - } - - - /* - /// - /// Computes the modified Bessel function of the first kind. - /// - /// The order. - /// The argument. - /// The value of In(x). - /// - /// The modified Bessel function are essentially the Bessel functions with an imaginary argument. - /// In particular, In(x) = (-1)n Jn(i x). - /// - /// - public static double ScaledBesselI (int n, double x) { - - if (n < 0) { - return (ScaledBesselI(-n, x)); - } - - if (x < 0) { - double I = ScaledBesselI(n, -x); - if (n % 2 != 0) I = -I; - return (I); - } - - if (x <= 2.0 * Math.Sqrt(n + 1)) { - return (Math.Exp(-x) * BesselI_Series(n, x)); - } else if (x > 20.0 + n * n / 4.0) { - return (Math.Exp(x) * BesselI_Reduced_Asymptotic(n, x)); - } else { - return (Math.Exp(x) * BesselI_Reduced_Miller(n, x)); - } - } - - private static double BesselI_Series (double nu, double x) { - - if (x == 0.0) { - if (nu == 0.0) { - return (1.0); - } else { - return (0.0); - } - } else { - double xh = 0.5 * x; - double dI = 1.0; - double I = dI; - double xx = xh * xh; - for (int k = 1; k < Global.SeriesMax; k++) { - double I_old = I; - dI = dI * xx / (k * (nu + k)); - I += dI; - if (I == I_old) { - return (PowerOverFactorial(xh, nu) * I); - } - } - throw new NonconvergenceException(); - } - } - - private static double BesselI_Reduced_Miller (int n, double x) { - - // use I_{k-1} = (2k/x) I_k + I_{k+1} to iterate down from a high k - // use I_0 + 2 I_1 + 2 I_2 + 2 I_3 + 2 I_4 = e^x to normalize - - // this should be stable for all x, since I_k increases as k decreases - - int m = Bessel_Miller_Limit(n, x); - - // assume zero values at a high order - double IP = 0.0; - double I = 1.0; - - // do the renormalization sum as we go - double sum = 0.0; - - // iterate down to 0 - double In = Double.NaN; - for (int k = m; k > 0; k--) { - if (k == n) { - In = I; - } - sum += I; - double IM = (2 * k / x) * I + IP; - IP = I; - I = IM; - } - - // add the I_0 term to sum - sum = 2.0 * sum + I; - - // since we are dividing by e^x, sum should be 1, so divide all terms by the actual sum - - return (In / sum); - - } - - private static double BesselI_Reduced_Asymptotic (double nu, double x) { - - // asmyptotic expansion for modified Bessel I - // development is ~ (nu^2 / x) - - double mu = 4.0 * nu * nu; - double xx = 8.0 * x; - - double dI = 1.0; - double I = dI; - for (int k = 1; k < Global.SeriesMax; k++) { - double I_old = I; - int kk = 2 * k - 1; - dI = - dI * (mu - kk * kk) / xx / k; - I += dI; - if (I == I_old) { - return (I / Math.Sqrt(2.0 * Math.PI * x)); - } - } - - // can rewrite as dI = - dI * a * b / xx / k, where a+=2 and b -= 2 for each iteration - - throw new NonconvergenceException(); - } - */ - - // Neuman series no good for computing K0; it relies on a near-perfect cancelation between the I0 term - // and the higher terms to achieve an exponentially supressed small value - - } - -} +using System; +using System.Diagnostics; + +using Meta.Numerics; + +namespace Meta.Numerics.Functions { + + // mostly direct and straightforward adaptations of the Bessel routines + + public static partial class AdvancedMath { + + /// + /// Computes modified cylindrical Bessel functions. + /// + /// The order, which must be non-negative. + /// The argument, which must be non-negative. + /// The values of I, I', K, and K' for the given order and argument. + /// + /// The modified Bessel functions fulfill a differential equation similar to the Bessel differential equation. + /// + /// + /// + public static SolutionPair ModifiedBessel (double nu, double x) { + if (nu < 0.0) throw new ArgumentOutOfRangeException(nameof(nu)); + if (x < 0.0) throw new ArgumentOutOfRangeException(nameof(x)); + + if (x == 0.0) { + return Bessel_Zero(nu, +1); + } else if (x < 2.0) { + + // Use the series to determine I and I' + ModifiedBesselI_Series(nu, x, out double I, out double IP); + + // Use the series to determine K and K' at -1/2 <= mu <= 1/2 that is an integer offset from nu + int n = (int) Math.Round(nu); + double mu = nu - n; + Debug.Assert(Math.Abs(mu) <= 0.5); + + // Determine K and K' at order mu + ModifiedBesselK_Series(mu, x, out double K, out double K1); + double KP = (mu / x) * K - K1; + + // Recurr K and K' upward to order nu + ModifiedBesselK_RecurrUpward(mu, x, ref K, ref KP, n); + + // Return the result + return (new SolutionPair(I, IP, K, KP)); + + } else if (x > 32.0 + 0.5 * nu * nu) { + ModifiedBessel_Asymptotic_Scaled(nu, x, out double sI, out double sIP, out double sK, out double sKP); + double e = Math.Exp(x); + return (new SolutionPair(e * sI, e * sIP, sK / e, sKP / e)); + + } else { + + // find 0 <= mu < 1 with same fractional part as nu + // this is necessary because CF2 does not produce K with good accuracy except at very low orders + int n = (int) Math.Floor(nu); + double mu = nu - n; + + // compute K, K' at this point (which is beyond the turning point because mu is small) using CF2 + ModifiedBessel_CF_K(mu, x, out double sK, out double g); + double sKP = g * sK; + + // recurr upward to order nu + ModifiedBesselK_RecurrUpward(mu, x, ref sK, ref sKP, n); + + // determine I'/I at the desired point + double f = ModifiedBessel_CF1(nu, x); + + // Use the Wronskian relationship K I' - I K' = 1/x to determine I and I' separately + double sI = 1.0 / (f * sK - sKP) / x; + double sIP = f * sI; + + double e = Math.Exp(x); + return (new SolutionPair(e * sI, e * sIP, sK / e, sKP / e)); + + } + + } + + // The recurrence for the modified Bessel functions (A&S 9.6.26, DLMF 10.29.2) is + // F_{\nu + 1} = F_{\nu}' - (\nu / x) F + // It is true for F = I or (-1)^{\nu} K, but since I is supressed as \nu increases, it should only be used for K upward, + // and only be used for I downward. + + // Multiplying by e^{\pm x}, we see that the exact same recurrence applies to sI or sK. + + // If you can compute K_0 and K_1, it can be started using K_0' = -K_1. + + private static void ModifiedBesselK_RecurrUpward (double mu, double x, ref double K, ref double KP, int n) { + for (int i = 0; i < n; i++) { + double t = K; + K = (mu / x) * K - KP; + mu += 1.0; + KP = -(t + (mu / x) * K); + } + } + + /// + /// Computes the regular modified cylindrical Bessel function. + /// + /// The order parameter. + /// The argument, which must be non-negative. + /// The value of Iν(x). + /// + /// The modified Bessel functions appear as the solutions of hyperbolic differential equations with + /// cylindrical or circular symmetry, for example the conduction of heat through a cylindrical pipe. + /// The regular modified Bessel functions are related to the Bessel functions with pure imaginary arguments. + /// + /// The regular modified Bessel functions increase monotonically and exponentially from the origin. + /// Because they increase exponentially, this function overflows for even moderately large arguments. In this + /// regeime, you can still obtain the value of the scaled modified e-x Iν(x) + /// by calling . + /// + /// is negative. + /// + /// + /// + public static double ModifiedBesselI (double nu, double x) { + return ModifiedBesselI(nu, x, false); + } + + /// + /// Computes the exponentially re-scaled regular modified cylindrical Bessel function. + /// + /// The order parameter. + /// The argument, which must be non-negative. + /// The value of e-x Iν(x). + /// + /// + public static double ScaledModifiedBesselI (double nu, double x) { + return ModifiedBesselI(nu, x, true); + } + + private static double ModifiedBesselI (double nu, double x, bool scaled) { + + if (x < 0.0) throw new ArgumentOutOfRangeException(nameof(x)); + + // Reflect negative nu to positive nu. + if (nu < 0.0) { + if (Math.Round(nu) == nu) { + return (ModifiedBesselI(-nu, x)); + } else { + return (ModifiedBesselI(-nu, x) + 2.0 / Math.PI * MoreMath.SinPi(-nu) * ModifiedBesselK(-nu, x)); + } + } + + if (x == 0.0) { + if (nu == 0.0) { + return (1.0); + } else { + return (0.0); + } + } else if (x <= 2.0 * Math.Sqrt(nu + 1.0)) { + // Close to the origin, use the power series. + double I = ModifiedBesselI_Series(nu, x); + return scaled ? Math.Exp(-x) * I : I; + } else if (x < 32.0 + 0.5 * nu * nu) { + // In the intermediate region, we will use CF1 + CF2 in a slightly different way than for the standard Bessel functions. + + // Find 0 <= mu < 1 with same fractional part as nu. + // This is necessary because CF2 does not produce K with good accuracy except at very low orders. + int n = (int) Math.Floor(nu); + double mu = nu - n; + Debug.Assert(0.0 <= mu && mu < 1.0); + + // Compute K, K' at this point (which is beyond the turning point because mu is small) using CF2. + ModifiedBessel_CF_K(mu, x, out double sK, out double g); + double sKP = g * sK; + + // Recurr upward to the desired order. + // This is okay because K increases with increasing order. + ModifiedBesselK_RecurrUpward(mu, x, ref sK, ref sKP, n); + + // Determine I'/I at the desired point. + double f = ModifiedBessel_CF1(nu, x); + + // Use the Wronskian to determine I from (I'/I), (K'/K), and K. + double sI = 1.0 / (f * sK - sKP) / x; + + return scaled ? sI : Math.Exp(x) * sI; + } else { + // Far from the origin, use the asymptotic expansion. + ModifiedBessel_Asymptotic_Scaled(nu, x, out double sI, out double sIP, out double sK, out double sKP); + return scaled ? sI : Math.Exp(x) * sI; + } + + } + + /// + /// Computes the irregular modified cylindrical Bessel function. + /// + /// The order parameter. + /// The argument, which must be non-negative. + /// The value of Kν(x). + /// + /// The modified Bessel functions are related to the Bessel functions with pure imaginary arguments. + /// The irregular modified Bessel function decreases monotonically and exponentially from the origin. + /// + /// is negative. + /// + public static double ModifiedBesselK (double nu, double x) { + return ModifiedBesselK(nu, x, false); + } + + /// + /// Computes the exponentially re-scaled irregular modified cylindrical Bessel function. + /// + /// The order parameter. + /// The argument, which must be non-negative. + /// The value of ex Kν(x). + public static double ScaledModifiedBesselK (double nu, double x) { + return ModifiedBesselK(nu, x, true); + } + + private static double ModifiedBesselK (double nu, double x, bool scaled) { + + if (x < 0.0) throw new ArgumentOutOfRangeException(nameof(x)); + + if (nu < 0.0) return (ModifiedBesselK(-nu, x)); + + if (x == 0.0) { + return (Double.PositiveInfinity); + } else if (x <= 2.0) { + // For small x, determine the value at small mu using the series and recurr up. + int n = (int) Math.Round(nu); + double mu = nu - n; + Debug.Assert(n >= 0); + Debug.Assert(Math.Abs(mu) <= 0.5); + + ModifiedBesselK_Series(mu, x, out double K, out double K1); + + if (n == 0) { + // K is already for order nu + } else if (n == 1) { + K = K1; + } else { + double KP = (mu / x) * K - K1; + ModifiedBesselK_RecurrUpward(mu, x, ref K, ref KP, n); + } + + return scaled ? Math.Exp(x) * K : K; + } else if (x < 32.0 + 0.5 * nu * nu) { + // In the intermediate region, determine the value at small mu using continued fraction and recurr up. + int n = (int) Math.Round(nu); + double mu = nu - n; + Debug.Assert(n >= 0); + Debug.Assert(Math.Abs(mu) <= 0.5); + + ModifiedBessel_CF_K(mu, x, out double sK, out double g); + double sKP = g * sK; + ModifiedBesselK_RecurrUpward(mu, x, ref sK, ref sKP, n); + return scaled ? sK : Math.Exp(-x) * sK; + } else { + // For large x, use the asymptotic series. + ModifiedBessel_Asymptotic_Scaled(nu, x, out double sI, out double sIP, out double sK, out double sKP); + return scaled ? sK : Math.Exp(-x) * sK; + } + + } + + // series near the origin; this is entirely analogous to the Bessel series near the origin + // it has a corresponding radius of rapid convergence, x < 4 + 2 Sqrt(nu) + + // This is exactly the same as BesselJ_Series with xx -> -xx. + // We could even factor this out into a common method with an additional parameter. + + private static void ModifiedBesselI_Series (double nu, double x, out double I, out double IP) { + Debug.Assert(x > 0.0); + + double x2 = 0.5 * x; + double xx = x2 * x2; + double dI = AdvancedMath.PowerOverFactorial(x2, nu); + I = dI; + IP = nu * dI; + for (int k = 1; k < Global.SeriesMax; k++) { + double I_old = I; + double IP_old = IP; + dI *= xx / (k * (nu + k)); + I += dI; + IP += (2 * k + nu) * dI; + if ((I == I_old) && (IP == IP_old)) { + IP = IP / x; + return; + } + } + + throw new NonconvergenceException(); + } + + private static double ModifiedBesselI_Series (double nu, double x) { + + double x2 = 0.5 * x; + double xx = x2 * x2; + double dI = AdvancedMath.PowerOverFactorial(x2, nu); + double I = dI; + for (int k = 1; k < Global.SeriesMax; k++) { + double I_old = I; + dI = dI * xx / (k * (nu + k)); + I += dI; + if (I == I_old) { + return (I); + } + } + throw new NonconvergenceException(); + + } + + // good for x > 32 + nu^2 / 2 + // this returns scaled values; I = e^(x) sI, K = e^(-x) sK, I' = e^(x) sI', K' = e^(-x) K' + // Note I'/K' != d/dx I/K under this scaling convention + + // Asymptotic expansions + // I_{\nu}(x) = \frac{e^x}{\sqrt{2 \pi x}} \left[ 1 - \frac{\mu-1}{8x} + \frac{(\mu-1)(\mu-9)}{2! (8x)^2} + \cdots \right] + // K_{\nu}(x) = \sqrt{\frac{\pi}{2x}}e^{-x} \left[ 1 + \frac{\mu-1}{8x} + \frac{(\mu-1)(\mu-9)}{2! (8x)^2} + \cdots \right] + // where \mu = 4 \nu^2. Derivatives + // I_{\nu}'(x) = \frac{e^x}{\sqrt{2 \pi x}} \left[ 1 - \frac{\mu+3}{8x} + \frac{(\mu-1)(\mu+15)}{2! (8x)^2} + \cdots \right] + // K_{\nu}'(x) = -\sqrt{\frac{\pi}{2x}}e^{-x} \left[ 1 + \frac{\mu+3}{8x} + \frac{(\mu-1)(\mu+15)}{2! (8x)^2} + \cdots \right] + // Note series differ only by alternating vs. same sign of terms. + + private static void ModifiedBessel_Asymptotic_Scaled (double nu, double x, out double sI, out double sIP, out double sK, out double sKP) { + + // precompute some values we will use + double mu = 4.0 * nu * nu; + double xx = 8.0 * x; + + // initialize the series + sI = 1.0; + sIP = 1.0; + sK = 1.0; + sKP = 1.0; + + // intialize the term value + double t = 1.0; + + for (int k = 1; k < Global.SeriesMax; k++) { + + double sI_old = sI; double sK_old = sK; + double sIP_old = sIP; double sKP_old = sKP; + + // determine next term values + int k2 = 2 * k; + t /= k * xx; + double tp = (mu + (k2 * k2 - 1)) * t; + //t *= (mu - (k2-1) * (k2-1)); + t *= (2.0 * (nu - k) + 1.0) * (2.0 * (nu + k) - 1.0); + + // add them, with alternating-sign for I series and same-sign for K series + if (k % 2 == 0) { + sI += t; + sIP += tp; + } else { + sI -= t; + sIP -= tp; + } + sK += t; + sKP += tp; + + // check for convergence + if ((sI == sI_old) && (sK == sK_old) && (sIP == sIP_old) && (sKP == sKP_old)) { + double fI = Math.Sqrt(2.0 * Math.PI * x); + double fK = Math.Sqrt(0.5 * Math.PI / x); + sI /= fI; + sIP /= fI; + sK *= fK; + sKP *= -fK; + return; + } + + } + + throw new NonconvergenceException(); + + } + + // compute I'/I via continued fraction + // this is analogous to CF1 for normal Bessel functions + // it converges quickly for x < nu and slowly for x > nu + + private static double ModifiedBessel_CF1 (double nu, double x) { + double q = 2.0 * (nu + 1.0); + double a = 1.0; + double b = (q / x); + double D = 1.0 / b; + double Df = a / b; + double f = nu / x + Df; + for (int k = 2; k < Global.SeriesMax; k++) { + double f_old = f; + + q += 2.0; + a = 1.0; + b = q / x; + D = 1.0 / (b + a * D); + Df = (b * D - 1.0) * Df; + f += Df; + + if (f == f_old) { + return (f); + } + } + + throw new NonconvergenceException(); + } + + // the continued fraction here is + + private static void ModifiedBessel_CF_K (double nu, double x, out double sK, out double g) { + + double a = 1.0; + double b = 2.0 * (1.0 + x); + + double D = 1.0 / b; + double Df = a / b; + double f = Df; + + // Thompson-Barnett do normalizing sum using these auxiliary variables (see Numerical Recipes) + + double C = (0.5 - nu) * (0.5 + nu); + double q0 = 0.0; + double q1 = 1.0; + double Q = C * q1; + double S = 1.0 + Q * Df; + + for (int k = 2; k < Global.SeriesMax; k++) { + + double f_old = f; + + // Lenz method for z1/z0 + + a = -(k - nu - 0.5) * (k + nu - 0.5); + b += 2.0; + + D = 1.0 / (b + a * D); + Df = (b * D - 1.0) * Df; + f += Df; + + // Thompson-Barnett sum trick + + double S_old = S; + + C = -a / k * C; + double q2 = (q0 - (b - 2.0) * q1) / a; + Q += C * q2; + S += Q * Df; + + if ((S == S_old) && (f == f_old)) { + g = -((x + 0.5) + (nu + 0.5) * (nu - 0.5) * f) / x; + sK = Math.Sqrt(Math.PI / 2.0 / x) / S; + return; + } + + q0 = q1; q1 = q2; + + } + + throw new NonconvergenceException(); + } + + // Temme's series for K0 and K1 for -1/2 \le \ny \le 1/2 and small x is described in NR3 6.6.37-6.6.40. + + // K_{\nu} = \sum_{k=0}^{\infty} c_k f_k, K_{\nu+1} = \frac{2}{x} \sum_{k=0}^{\infty} c_k h_k + + // Here c_k = \frac{(x/2)^{2k}}{k!} and f and h are determined by recurrences + + // p_{k} = \frac{p_{k-1}}{k - \nu} + // q_{k} = \frac{q_{k-1}}{k + \nu} + // f_{k} = \frac{k f_{k-1} + p_{k-1} + q_{k-1}}{k^2 - \nu^2} + // h_{k} = p_{k} - k f_{k} + + // with initial values + + // p_0 = \frac{1}{2} \left( \frac{x}{2} \right)^{-\nu} \Gamma(1 + \nu) + // q_0 = \frac{1}{2} \left( \frac{x}{2} \right)^{\nu} \Gamma(1 - \nu) + // f_0 = \frac{\nu\pi}{\sin(\nu\pi)} \left[ ... \right] + + // I observe it to be accurate to all but last couple digits for x < 2; it converges further out but + // its accuracy deteriorates rapidly. + + // The initial constant determination logic is shared with BesselY_Series; we should factor out the common parts. + + private static void ModifiedBesselK_Series (double nu, double x, out double K0, out double K1) { + + Debug.Assert((-0.5 <= nu) && (nu <= 0.5)); + Debug.Assert(x > 0.0); + + double x2 = 0.5 * x; + double x4 = x2 * x2; + + // determine initial p, q + + double GP = Gamma(1.0 + nu); + double GM = Gamma(1.0 - nu); + double z = Math.Pow(x2, nu); + double p = 0.5 / z * GP; + double q = 0.5 * z * GM; + + // determine initial f; this is rather complicated + + double ln = -Math.Log(x2); + + double s = nu * ln; + double cosh = Math.Cosh(s); + double sinhc = (s == 0.0) ? 1.0 : Math.Sinh(s) / s; + + double G1 = (Math.Abs(nu) < 0.25) ? NewG(1.0 - nu, 2.0 * nu) : (1.0 / GM - 1.0 / GP) / (2.0 * nu); + double G2 = (1.0 / GM + 1.0 / GP) / 2.0; + double F = (nu == 0.0) ? 1.0 : (nu * Math.PI) / MoreMath.SinPi(nu); + + double f = F * (cosh * G1 + sinhc * G2 * ln); + + // run the series + + double c = 1.0; + K0 = f; + K1 = p; + for (int k = 1; k < Global.SeriesMax; k++) { + + double K0_old = K0; + double K1_old = K1; + + c *= x4 / k; + + double kpnu = k + nu; + double kmnu = k - nu; + f = (k * f + p + q) / (kpnu * kmnu); + p = p / kmnu; + q = q / kpnu; + double h = p - k * f; + + K0 += c * f; + K1 += c * h; + + if ((K0 == K0_old) && (K1 == K1_old)) { + K1 = 2.0 / x * K1; + return; + } + } + throw new NonconvergenceException(); + + } + + + /* + /// + /// Computes the modified Bessel function of the first kind. + /// + /// The order. + /// The argument. + /// The value of In(x). + /// + /// The modified Bessel function are essentially the Bessel functions with an imaginary argument. + /// In particular, In(x) = (-1)n Jn(i x). + /// + /// + public static double ScaledBesselI (int n, double x) { + + if (n < 0) { + return (ScaledBesselI(-n, x)); + } + + if (x < 0) { + double I = ScaledBesselI(n, -x); + if (n % 2 != 0) I = -I; + return (I); + } + + if (x <= 2.0 * Math.Sqrt(n + 1)) { + return (Math.Exp(-x) * BesselI_Series(n, x)); + } else if (x > 20.0 + n * n / 4.0) { + return (Math.Exp(x) * BesselI_Reduced_Asymptotic(n, x)); + } else { + return (Math.Exp(x) * BesselI_Reduced_Miller(n, x)); + } + } + + private static double BesselI_Series (double nu, double x) { + + if (x == 0.0) { + if (nu == 0.0) { + return (1.0); + } else { + return (0.0); + } + } else { + double xh = 0.5 * x; + double dI = 1.0; + double I = dI; + double xx = xh * xh; + for (int k = 1; k < Global.SeriesMax; k++) { + double I_old = I; + dI = dI * xx / (k * (nu + k)); + I += dI; + if (I == I_old) { + return (PowerOverFactorial(xh, nu) * I); + } + } + throw new NonconvergenceException(); + } + } + + private static double BesselI_Reduced_Miller (int n, double x) { + + // use I_{k-1} = (2k/x) I_k + I_{k+1} to iterate down from a high k + // use I_0 + 2 I_1 + 2 I_2 + 2 I_3 + 2 I_4 = e^x to normalize + + // this should be stable for all x, since I_k increases as k decreases + + int m = Bessel_Miller_Limit(n, x); + + // assume zero values at a high order + double IP = 0.0; + double I = 1.0; + + // do the renormalization sum as we go + double sum = 0.0; + + // iterate down to 0 + double In = Double.NaN; + for (int k = m; k > 0; k--) { + if (k == n) { + In = I; + } + sum += I; + double IM = (2 * k / x) * I + IP; + IP = I; + I = IM; + } + + // add the I_0 term to sum + sum = 2.0 * sum + I; + + // since we are dividing by e^x, sum should be 1, so divide all terms by the actual sum + + return (In / sum); + + } + + private static double BesselI_Reduced_Asymptotic (double nu, double x) { + + // asmyptotic expansion for modified Bessel I + // development is ~ (nu^2 / x) + + double mu = 4.0 * nu * nu; + double xx = 8.0 * x; + + double dI = 1.0; + double I = dI; + for (int k = 1; k < Global.SeriesMax; k++) { + double I_old = I; + int kk = 2 * k - 1; + dI = - dI * (mu - kk * kk) / xx / k; + I += dI; + if (I == I_old) { + return (I / Math.Sqrt(2.0 * Math.PI * x)); + } + } + + // can rewrite as dI = - dI * a * b / xx / k, where a+=2 and b -= 2 for each iteration + + throw new NonconvergenceException(); + } + */ + + // Neuman series no good for computing K0; it relies on a near-perfect cancelation between the I0 term + // and the higher terms to achieve an exponentially supressed small value + + } + +} diff --git a/Numerics/Functions/AdvancedMath_Polylog.cs b/Numerics/Functions/AdvancedMath_Polylog.cs index 2143472..6f0ebb5 100644 --- a/Numerics/Functions/AdvancedMath_Polylog.cs +++ b/Numerics/Functions/AdvancedMath_Polylog.cs @@ -100,6 +100,9 @@ private static double DiLog_Series_1 (double e) { /// series. For n = 1 it reduces to -log(1-x). For n = 2 it reduces to the function. /// The polylogarithm function becomes complex for arguments larger than one. /// + /// is negative or is greater than one. + /// + /// public static double PolyLog (int n, double x) { if (x > 1.0) throw new ArgumentOutOfRangeException(nameof(x)); diff --git a/Numerics/Functions/AdvancedMath_Riemann.cs b/Numerics/Functions/AdvancedMath_Riemann.cs index e9495c4..a2d1707 100644 --- a/Numerics/Functions/AdvancedMath_Riemann.cs +++ b/Numerics/Functions/AdvancedMath_Riemann.cs @@ -16,19 +16,21 @@ public partial class AdvancedMath { /// /// /// + /// + /// public static double RiemannZeta (double x) { if (x < 0.0) { // for negative numbers, use the reflection formula double t = 1.0 - x; - return (2.0 * Math.Pow(Global.TwoPI, -t) * MoreMath.Cos(Global.HalfPI * t) * AdvancedMath.Gamma(t) * RiemannZeta(t)); + return 2.0 * Math.Pow(2.0 * Math.PI, -t) * MoreMath.CosPi(0.5 * t) * AdvancedMath.Gamma(t) * RiemannZeta(t); } else { double xm1 = x - 1.0; if (Math.Abs(xm1) < 0.25) { // near the singularity, use the Stjielts expansion - return (RiemannZeta_LaurentSeries(xm1)); + return RiemannZeta_LaurentSeries(xm1); } else { // call Dirichlet function, which converges faster - return (DirichletEta(x) / (1.0 - Math.Pow(2.0, 1.0 - x))); + return DirichletEta(x) / (1.0 - Math.Pow(2.0, 1.0 - x)); } } } @@ -51,7 +53,7 @@ public static double RiemannZeta (double x) { /// public static double DirichletEta (double x) { if (x < 0.0) throw new ArgumentOutOfRangeException(nameof(x)); - return (DirichletEta_Borwein(x)); + return DirichletEta_Borwein(x); } // Borwein's amazing method for computing eta is detailed at http://numbers.computation.free.fr/Constants/Miscellaneous/zetaevaluations.html. diff --git a/Numerics/Numerics.csproj b/Numerics/Numerics.csproj index a14e03e..4fab0b3 100644 --- a/Numerics/Numerics.csproj +++ b/Numerics/Numerics.csproj @@ -1,37 +1,37 @@ - - - - netstandard1.1 - Meta.Numerics - Meta.Numerics - 4.1.0-alpha - Copyright © David Wright 2008-2020 - David Wright - Meta Numerics - A library for numeric computing with support for data manipulation, statistics, matrix algebra, and advanced functions. - http://www.meta-numerics.net - numeric scientific technical math statistics matrix fft data optimization - http://www.meta-numerics.net/Images/ComplexPsiPlot.png - https://github.com/dcwuser/metanumerics - portable - true - https://opensource.org/licenses/MS-PL - false - - - - TRACE;RELEASE;NETSTANDARD1_1 - Numerics1.ruleset - true - - - - Numerics1.ruleset - - - - - - - - + + + + netstandard1.1 + Meta.Numerics + Meta.Numerics + 4.1.2-alpha + Copyright © David Wright 2008-2020 + David Wright + Meta Numerics + A library for numeric computing with support for data manipulation, statistics, matrix algebra, and advanced functions. + http://www.meta-numerics.net + numeric scientific technical math statistics matrix fft data optimization + http://www.meta-numerics.net/Images/ComplexPsiPlot.png + https://github.com/dcwuser/metanumerics + portable + true + https://opensource.org/licenses/MS-PL + false + + + + TRACE;RELEASE;NETSTANDARD1_1 + Numerics1.ruleset + true + + + + Numerics1.ruleset + + + + + + + + diff --git a/README.md b/README.md index 9d52de6..cb912e0 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,11 @@ -# metanumerics -Meta.Numerics is library for advanced numerical computing on the .NET platform. It offers an object-oriented API for statistical analysis, advanced functions, numerical integration and optimization, Fourier transforms, and matrix algebra. +# Meta.Numerics -For more information, visit http://www.meta-numerics.net. +Meta.Numerics is a library for advanced numerical computing for the .NET platform. +It offers an object-oriented API for data manipulation, statistical analysis, advanced functions, +matrix algebra, Fourier transforms, advanced functions, extended precision arithmetic, +and solver functionality such as integration, optimization, and root finding. + +Meta.Numerics is copyright 2008-2020 by David Wright. It is licensed under the +Microsoft Public License, which is a BSD-style open-source license. -Meta.Numerics is copyright 2008-2018 by David Wright +For more information, visit http://www.meta-numerics.net. diff --git a/ReleaseNotes.md b/ReleaseNotes.md new file mode 100644 index 0000000..abf8d3d --- /dev/null +++ b/ReleaseNotes.md @@ -0,0 +1,18 @@ +# Meta.Numerics 4.1 + +## More Extended Precision Types +We have made multiple improvements to the the quadruple-precision (~32 decimal digit) DoubleDouble floating point type. It now handles infinities and NaNs appropriately. Its string rendering and parsing are improved. We have also added more DoubleDouble functions, including trig functions and LogGamma. + +We have added the 128-bit integer types Int128 and UInt128. These behave like the built-in integer types Int64 and UInt64 (long and ulong), but support integer values up to ~1038. Arithmetic using these types is 1-4 times faster than using BigInteger, and unlike BigInteger, they behave like the other fixed-width register types with respect to overflow. + +## More Advanced Functions +We have added a few more advanced functions. These include the complete elliptic integral of the third kind and computation of the elliptic nome (so now only the incomplete elliptic integral of the third kind remains unimplemented). We also added a scaled version of the incomplete Bessel function (allowing you to work with the function for arguments where the function value itself would over- or under-flow), functions that return the zeros of the Airy and Bessel functions, and the hyperbolic integral functions Cin and Shi. + +We have also made many improvements to the internals of long-implemented functions to improve their speed, accuracy, and behavior at extreme arguments including infinities and NaN. + +## Other Improvements +We have added the RegressionResult type with exposes residuals and the sum of squared residuals on all FitResults that have them. We fixed a bug which could cause non-convergence in the multi-dimensional FindLocalMaximum and FindLocalMinimum methods. + +## Future Plans +The next release of the Meta.Numerics will use .NET Standard 2.1. +We will eliminate classes like Sample, BivariateSample, and MultivariateSample in favor of the Univariate, Bivariate, and Multivariate extension method classes which can be used with any type implementing IReadOnlyList, including the built-in collection types and our own DataFrame types. \ No newline at end of file diff --git a/Test/AdvancedComplexMathTest.cs b/Test/AdvancedComplexMathTest.cs index 76e4481..f040a01 100644 --- a/Test/AdvancedComplexMathTest.cs +++ b/Test/AdvancedComplexMathTest.cs @@ -172,6 +172,13 @@ public void ComplexLogGammaRecurrance () { } } + [TestMethod] + public void ComplexLogGammaSpecialValues () { + Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedComplexMath.LogGamma(-0.5), new Complex(Math.Log(2.0 * Math.Sqrt(Math.PI)), -Math.PI))); + Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedComplexMath.LogGamma(-1.5), new Complex(Math.Log(4.0 * Math.Sqrt(Math.PI) / 3.0), -2.0 * Math.PI))); + Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedComplexMath.LogGamma(-2.5), new Complex(Math.Log(8.0 * Math.Sqrt(Math.PI) / 15.0), -3.0 * Math.PI))); + } + // periodicity in imaginary part of ln z means that LogGamma recurrence and LogGamma duplication tests fail [TestMethod] diff --git a/Test/AdvancedIntegerMathTest.cs b/Test/AdvancedIntegerMathTest.cs index 37730a8..3e9beb2 100644 --- a/Test/AdvancedIntegerMathTest.cs +++ b/Test/AdvancedIntegerMathTest.cs @@ -516,6 +516,22 @@ public void BernoulliNumbers () { Assert.IsTrue(AdvancedIntegerMath.BernoulliNumber(1) == -1.0 / 2.0); Assert.IsTrue(AdvancedIntegerMath.BernoulliNumber(2) == 1.0 / 6.0); Assert.IsTrue(AdvancedIntegerMath.BernoulliNumber(3) == 0.0); + Assert.IsTrue(AdvancedIntegerMath.BernoulliNumber(4) == -1.0 / 30.0); + } + + //[TestMethod] + // Cancellation too bad to be useful + public void BernoulliStirlingRelationship () { + // This involves significant cancellation, so don't pick n too high + foreach (int n in TestUtilities.GenerateIntegerValues(2, 16, 4)) { + double S = 0.0; + for (int k = 0; k <= n; k++) { + double dS = AdvancedIntegerMath.Factorial(k) / (k + 1) * AdvancedIntegerMath.StirlingNumber2(n, k); + if (k % 2 != 0) dS = -dS; + S += dS; + } + Assert.IsTrue(TestUtilities.IsNearlyEqual(S, AdvancedIntegerMath.BernoulliNumber(n))); + } } } diff --git a/Test/AdvancedMathTest.cs b/Test/AdvancedMathTest.cs index f77ae08..ca60467 100644 --- a/Test/AdvancedMathTest.cs +++ b/Test/AdvancedMathTest.cs @@ -740,29 +740,28 @@ public void IntegralEInvalidArgumentTest () { [TestMethod] public void RiemannZetaSpecialCaseTest () { Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.RiemannZeta(-3.0), 1.0 / 120.0)); - //Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.RiemannZeta(-2.0), 0.0)); + Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.RiemannZeta(-2.0), 0.0)); Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.RiemannZeta(-1.0), -1.0 / 12.0)); Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.RiemannZeta(0.0), -0.5)); Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.RiemannZeta(2.0), Math.PI * Math.PI / 6.0)); Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.RiemannZeta(4.0), Math.Pow(Math.PI, 4) / 90.0)); } + // DLMF 25.4.1 [TestMethod] public void RiemannZetaReflectionTest () { foreach (double x in TestUtilities.GenerateRealValues(1.0E-2 ,1.0E2, 16)) { - Console.WriteLine("x={0}", x); double zx = AdvancedMath.RiemannZeta(x); double gx = AdvancedMath.Gamma(x); double zr = AdvancedMath.RiemannZeta(1.0 - x); - Console.WriteLine(" {0} vs. {1}", zr, 2.0 * Math.Pow(2.0 * Math.PI, -x) * Math.Cos(Math.PI * x / 2.0) * gx * zx); - Assert.IsTrue(TestUtilities.IsNearlyEqual(zr, 2.0 * Math.Pow(2.0 * Math.PI, -x) * Math.Cos(Math.PI * x / 2.0) * gx * zx)); + Assert.IsTrue(TestUtilities.IsNearlyEqual(zr, 2.0 * Math.Pow(2.0 * Math.PI, -x) * MoreMath.CosPi(x / 2.0) * gx * zx)); } } - + // A&S 23.2.2, DLMF 25.2.11 [TestMethod] public void ReimannZetaPrimesTest () { - // pick high enough values so that p^-x == 1 within double precision before we reach the end of our list of primes + // pick high enough values so that 1 - p^-x == 1 within double precision, i.e. i.e. p^x > ~10^16, before we reach the end of our list of primes foreach (double x in TestUtilities.GenerateRealValues(10.0, 1.0E3, 8)) { double f = 1.0; foreach (int p in primes) { @@ -774,7 +773,43 @@ public void ReimannZetaPrimesTest () { } } - private int[] primes = new int[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 }; + private readonly int[] primes = new int[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 }; + + // A&S 23.2.15, 23.2.16 + [TestMethod] + public void ReimannZetaBernoulliRelationship () { + foreach (int n in TestUtilities.GenerateIntegerValues(2, 32, 8)) { + int m = 2 * n; + Assert.IsTrue(TestUtilities.IsNearlyEqual( + AdvancedMath.RiemannZeta(m), + MoreMath.Pow(2.0 * Math.PI, m) / AdvancedIntegerMath.Factorial(m) / 2.0 * Math.Abs(AdvancedIntegerMath.BernoulliNumber(m)) + )); + Assert.IsTrue(TestUtilities.IsNearlyEqual( + AdvancedMath.RiemannZeta(1.0 - m), + -AdvancedIntegerMath.BernoulliNumber(m) / m + )); + } + } + + + // DLMF 25.8.3 + [TestMethod] + public void ReimannZetaPochammerSum () { + // For s >~ 10, converges to nearly 1, and takes a long time to do so. + foreach (double s in TestUtilities.GenerateRealValues(0.1, 10.0, 4)) { + if (Math.Abs(s - 1.0) < 1.0E-2) continue; + double S = 0.0; + for (int k = 0; true; k++) { + if (k > 100) throw new NonconvergenceException(); + double S_old = S; + double dS = AdvancedMath.Pochhammer(s, k) * AdvancedMath.RiemannZeta(s + k) / AdvancedIntegerMath.Factorial(k) / Math.Pow(2.0, s + k); + S += dS; + if (S == S_old) break; + } + // RHS is just lambda function (Riemann sum over odd integers only), but we don't have a seperate function for that + Assert.IsTrue(TestUtilities.IsNearlyEqual(S, AdvancedMath.RiemannZeta(s) * (1.0 - Math.Pow(2.0, -s)))); + } + } [TestMethod] public void DirichletEtaSpecialCaseTest () { @@ -1166,6 +1201,31 @@ public void PolyLogAsFermiDiracIntegral () { } + // DLMF 19.20.1 + [TestMethod] + public void CarlsonFSpecialCases () { + foreach (double x in TestUtilities.GenerateRealValues(1.0E-2, 1.0E2, 4)) { + Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.CarlsonF(x, x, x), 1.0 / Math.Sqrt(x))); + Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.CarlsonF(x, x, 0.0), Math.PI / 2.0 / Math.Sqrt(x))); + Assert.IsTrue(AdvancedMath.CarlsonF(x, 0.0, 0.0) == Double.PositiveInfinity); + } + } + + [TestMethod] + public void CarlsonFScaling () { + foreach (double x in TestUtilities.GenerateRealValues(1.0E-2, 1.0E2, 3, 1)) { + foreach (double y in TestUtilities.GenerateRealValues(1.0E-3, 1.0E1, 3, 2)) { + foreach (double z in TestUtilities.GenerateRealValues(1.0E-1, 1.0E3, 3, 3)) { + foreach (double lambda in TestUtilities.GenerateRealValues(1.0E-1, 1.0E1, 3, 4)) { + Assert.IsTrue(TestUtilities.IsNearlyEqual( + AdvancedMath.CarlsonF(lambda * x, lambda * y, lambda * z), + AdvancedMath.CarlsonF(x, y, z) / Math.Sqrt(lambda) + )); + } + } + } + } + } [TestMethod] public void CarlsonFDuplication () { diff --git a/Test/AdvancedMathTest_Gamma.cs b/Test/AdvancedMathTest_Gamma.cs index b759a6d..a237add 100644 --- a/Test/AdvancedMathTest_Gamma.cs +++ b/Test/AdvancedMathTest_Gamma.cs @@ -34,7 +34,7 @@ public void GammaSpecialCases () { Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.Gamma(0.5), Math.Sqrt(Math.PI))); Assert.IsTrue(AdvancedMath.Gamma(1.0) == 1.0); Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.Gamma(1.5), Math.Sqrt(Math.PI) / 2.0)); - Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.Gamma(2.0), 1.0)); + Assert.IsTrue(AdvancedMath.Gamma(2.0) == 1.0); Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.Gamma(3.0), 2.0)); Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.Gamma(4.0), 6.0)); Assert.IsTrue(Double.IsPositiveInfinity(AdvancedMath.Gamma(Double.PositiveInfinity))); @@ -187,5 +187,52 @@ public void GammaRatioInequality () { } } + + [TestMethod] + public void PochammerSpecialIntegerArguments () { + Assert.IsTrue(AdvancedMath.Pochhammer(0.0, 0) == 1.0); + foreach (int n in TestUtilities.GenerateIntegerValues(1, 100, 6)) { + Assert.IsTrue(AdvancedMath.Pochhammer(0.0, n) == 0.0); + Assert.IsTrue(TestUtilities.IsNearlyEqual( + AdvancedMath.Pochhammer(0.0, -n), + (n % 2 == 0 ? 1.0 : -1.0) / AdvancedIntegerMath.Factorial(n) + )); + Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.Pochhammer(1.0, n), AdvancedIntegerMath.Factorial(n))); + } + } + + // DLMF 5.2.4-5.2.8 + [TestMethod] + public void PochammerIntegerArguments () { + foreach (double x in TestUtilities.GenerateRealValues(1.0E-2, 1.0E2, 5)) { + Assert.IsTrue(AdvancedMath.Pochhammer(x, 0) == 1.0); + foreach (int n in TestUtilities.GenerateIntegerValues(1, 100, 4)) { + double P = 1.0; + for (int k = 0; k < n; k++) { + P *= (x + k); + } + Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.Pochhammer(x, n), P)); + if (Math.Abs(x - Math.Round(x)) > 0.01) { + double Q = 1.0; + for (int k = 1; k <= n; k++) { + Q /= (x - k); + } + Assert.IsTrue(TestUtilities.IsNearlyEqual(AdvancedMath.Pochhammer(x, -n), Q)); + } + Assert.IsTrue(TestUtilities.IsNearlyEqual( + AdvancedMath.Pochhammer(-x, n), + (n % 2 == 0 ? 1 : -1) * AdvancedMath.Pochhammer(x - n + 1, n) + )); + Assert.IsTrue(TestUtilities.IsNearlyEqual( + AdvancedMath.Pochhammer(x, 2 * n), + MoreMath.Pow(2.0, 2 * n) * AdvancedMath.Pochhammer(x / 2.0, n) * AdvancedMath.Pochhammer((x + 1.0) / 2.0, n) + )); + } + + } + + + } + } } diff --git a/Test/DistributionTest.cs b/Test/DistributionTest.cs index 36a6ad8..cf3ae3d 100644 --- a/Test/DistributionTest.cs +++ b/Test/DistributionTest.cs @@ -47,7 +47,8 @@ private static List CreateDistributions () { new LaplaceDistribution(4.5, 6.0), new ChiDistribution(1), new ChiDistribution(4), new RayleighDistribution(3.0), - new FrechetDistribution(2.9, 4.0) + new FrechetDistribution(2.9, 4.0), + new TestDistribution() }); // Add some distributions that come from tests. @@ -763,10 +764,10 @@ public override Interval Support { public override double ProbabilityDensity (double x) { if (x < 0.0) { return (0.0); - } else if (x > 1.0) { - return (1.0); + } else if (x <= 1.0) { + return (2.0 * x); } else { - return (2.0 * x); + return (0.0); } } diff --git a/Test/DoubleDoubleTest.cs b/Test/DoubleDoubleTest.cs index c1be18e..818ee65 100644 --- a/Test/DoubleDoubleTest.cs +++ b/Test/DoubleDoubleTest.cs @@ -187,6 +187,11 @@ public void DoubleDoubleSerializationRoundtrip () { } + [TestMethod] + public void DoubleDoubleToStringSpecial () { + Assert.IsTrue(DoubleDouble.NaN.ToString() == Double.NaN.ToString()); + } + [TestMethod] public void DoubleDoubleParseSpecialCases () { Assert.IsTrue(DoubleDouble.Parse("0.0") == DoubleDouble.Zero); diff --git a/Test/Int128Test.cs b/Test/Int128Test.cs index c9d7271..5cea290 100644 --- a/Test/Int128Test.cs +++ b/Test/Int128Test.cs @@ -186,15 +186,38 @@ public void Int128RandomArithmetic () { // Subtraction and addition are inverses Assert.IsTrue(y + (x - y) == x); - // Division and multiplication are inverses + // Division methods agree Int128 q = Int128.DivRem(x, y, out Int128 r); Assert.IsTrue(q == x / y); + Assert.IsTrue(r == x % y); + + // Division and multiplication are inverses Assert.IsTrue(q * y + r == x); } } } + [TestMethod] + public void Int128Increment () { + foreach (Int128 u in GetRandomInt128(4)) { + Int128 v = u; + Int128 p = v++; + Assert.IsTrue(v == u + 1); + Assert.IsTrue(p == u); + + v = u; + p = v--; + Assert.IsTrue(v == u - 1); + Assert.IsTrue(p == u); + + v = u; + p = ++v; + Assert.IsTrue(v == u + 1); + Assert.IsTrue(p == v); + } + } + [TestMethod] public void DivisionByZero () { @@ -214,6 +237,19 @@ public void Int128Overflow () { } + [TestMethod] + public void Int128DoubleRoundtrip () { + // For integers exactly representable by doubles (52 or fewer binary digits), + // roundtrip through double should preserve value + Random rng = new Random(8); + foreach (Int128 s in GetRandomInt64(8, rng)) { + if (Int128.Abs(s) > (1UL << 51)) continue; + double d = (double) s; + Int128 s1 = (Int128) d; + Assert.IsTrue(s1 == s); + } + } + [TestMethod] public void Int128SepcialStrings () { @@ -337,17 +373,6 @@ public void Int128BigIntegerAgreement () { } - [TestMethod] - public void Int128DoubleRoundtrip () { - // If value is less than 2^52, we should be able to roundtrip to double. - foreach (long x in GetRandomInt64(8)) { - if (Math.Abs(x) < ((long) (1UL << 52))) { - Int128 y = (Int128) x; - Assert.IsTrue(((Int128) ((double) y)) == y); - } - } - } - [TestMethod] public void Int128Abs () { Assert.IsTrue(Int128.Abs(Int128.Zero) == Int128.Zero); diff --git a/Test/MomentMathTest.cs b/Test/MomentMathTest.cs index 8ae80a1..86128ba 100644 --- a/Test/MomentMathTest.cs +++ b/Test/MomentMathTest.cs @@ -9,6 +9,32 @@ namespace Test { [TestClass] public class MomentMathTest { + [TestMethod] + public void MomentMathOrderOne () { + + double mu = 2.0; + double[] M = new double[] { 1.0, mu }; + double[] C = MomentMath.RawToCentral(M); + Assert.IsTrue(C[0] == 1.0); + Assert.IsTrue(C[1] == 0.0); + double[] K = MomentMath.RawToCumulant(M); + Assert.IsTrue(K[1] == mu); + + M = MomentMath.CentralToRaw(mu, C); + Assert.IsTrue(M[0] == 1.0); + Assert.IsTrue(M[1] == mu); + K = MomentMath.CentralToCumulant(mu, C); + Assert.IsTrue(K[1] == mu); + + M = MomentMath.CumulantToRaw(K); + Assert.IsTrue(M[0] == 1.0); + Assert.IsTrue(M[1] == mu); + C = MomentMath.CumulantToCentral(K); + Assert.IsTrue(C[0] == 1.0); + Assert.IsTrue(C[1] == 0.0); + + } + [TestMethod] public void MomentMathConsistency () { diff --git a/Test/UInt128Test.cs b/Test/UInt128Test.cs index 658e089..b3c58c3 100644 --- a/Test/UInt128Test.cs +++ b/Test/UInt128Test.cs @@ -143,16 +143,41 @@ public void UInt128Arithmetic () { } [TestMethod] - public void UInt128Overflow () { + public void UInt128Increment () { + foreach (UInt128 u in GetRandomUInt128(new Random(1), 4)) { + UInt128 v = u; + UInt128 p = v++; + Assert.IsTrue(v == u + 1); + Assert.IsTrue(p == u); + + v = u; + p = v--; + Assert.IsTrue(v == u - 1); + Assert.IsTrue(p == u); + + v = u; + p = ++v; + Assert.IsTrue(v == u + 1); + Assert.IsTrue(p == v); + } + } + - UInt128 two = 2; + [TestMethod] + public void UInt128Overflow () { unchecked { Assert.IsTrue(UInt128.MaxValue + UInt128.One == UInt128.Zero); // 99 + 1 = [1]00 Assert.IsTrue(UInt128.MaxValue + UInt128.MaxValue == UInt128.MaxValue - UInt128.One); // 99 + 99 = [1]98 Assert.IsTrue(UInt128.Zero - UInt128.One == UInt128.MaxValue); // [1]00 - 1 = 99 Assert.IsTrue(UInt128.Zero - UInt128.MaxValue == UInt128.One); // [1]00 - 99 = 1 - Assert.IsTrue(UInt128.MaxValue * two == UInt128.MaxValue - UInt128.One); // 99 * 2 = [1]98 + Assert.IsTrue(UInt128.MaxValue * 2 == UInt128.MaxValue - UInt128.One); // 99 * 2 = [1]98 Assert.IsTrue(UInt128.MaxValue * UInt128.MaxValue == UInt128.One); // 99 * 99 = [98]01 + + UInt128 m = UInt128.MaxValue; + m++; + Assert.IsTrue(m == UInt128.Zero); + m--; + Assert.IsTrue(m == UInt128.MaxValue); } }