From c72bb3662ea04ffadb9a8b925247eeb8b37aacbf Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 10 Jan 2015 10:51:18 +0100 Subject: [PATCH] Precision: migrate epsilon logic from NumericalDerivative to Precision class --- .../Differentiation/NumericalDerivative.cs | 51 +++++--------- src/Numerics/Precision.cs | 66 +++++++++++++++++-- src/Numerics/RootFinding/Brent.cs | 2 +- 3 files changed, 81 insertions(+), 38 deletions(-) diff --git a/src/Numerics/Differentiation/NumericalDerivative.cs b/src/Numerics/Differentiation/NumericalDerivative.cs index a0bfbf9a..a0f88e12 100644 --- a/src/Numerics/Differentiation/NumericalDerivative.cs +++ b/src/Numerics/Differentiation/NumericalDerivative.cs @@ -45,23 +45,23 @@ namespace MathNet.Numerics.Differentiation Absolute, /// - /// A base step size value, h, will be scaled according to the function input parameter. A common example is hx = h*(1+abs(x)), however + /// A base step size value, h, will be scaled according to the function input parameter. A common example is hx = h*(1+abs(x)), however /// this may vary depending on implementation. This definition only guarantees that the only scaling will be relative to the /// function input parameter and not the order of the finite difference derivative. /// RelativeX, /// - /// A base step size value, eps (typically machine precision), is scaled according to the finite difference coefficient order - /// and function input parameter. The initial scaling according to finite different coefficient order can be thought of as producing a - /// base step size, h, that is equivalent to scaling. This stepsize is then scaled according to the function + /// A base step size value, eps (typically machine precision), is scaled according to the finite difference coefficient order + /// and function input parameter. The initial scaling according to finite different coefficient order can be thought of as producing a + /// base step size, h, that is equivalent to scaling. This stepsize is then scaled according to the function /// input parameter. Although implementation may vary, an example of second order accurate scaling may be (eps)^(1/3)*(1+abs(x)). /// Relative }; /// - /// Class to evaluate the numerical derivative of a function using finite difference approximations. + /// Class to evaluate the numerical derivative of a function using finite difference approximations. /// Variable point and center methods can be initialized . /// This class can also be used to return function handles (delagates) for a fixed derivative order and variable. /// It is possible to evaluate the derivative and partial derivative of univariate and multivariate functions respectively. @@ -69,12 +69,12 @@ namespace MathNet.Numerics.Differentiation public class NumericalDerivative { /// - /// Sets and gets the finite difference step size. This value is for each function evaluation if relative stepsize types are used. - /// If the base step size used in scaling is desired, see . + /// Sets and gets the finite difference step size. This value is for each function evaluation if relative stepsize types are used. + /// If the base step size used in scaling is desired, see . /// /// - /// Setting then getting the StepSize may return a different value. This is not unusual since a user-defined step size is converted to a - /// base-2 representable number to improve finite difference accuracy. + /// Setting then getting the StepSize may return a different value. This is not unusual since a user-defined step size is converted to a + /// base-2 representable number to improve finite difference accuracy. /// public double StepSize { @@ -138,7 +138,7 @@ namespace MathNet.Numerics.Differentiation /// /// Type of step size for computing finite differences. If set to absolute, dx = h. - /// If set to relative, dx = (1+abs(x))*h^(2/(order+1)). This provides accurate results when + /// If set to relative, dx = (1+abs(x))*h^(2/(order+1)). This provides accurate results when /// h is approximately equal to the square-root of machine accuracy, epsilon. /// public StepType StepType @@ -174,7 +174,7 @@ namespace MathNet.Numerics.Differentiation throw new ArgumentOutOfRangeException("points", "Points must be two or greater."); _points = points; Center = center; - _epsilon = CalculateMachineEpsilon(); + _epsilon = Precision.PositiveMachineEpsilon; _coefficients = new FiniteDifferenceCoefficients(points); } @@ -224,7 +224,7 @@ namespace MathNet.Numerics.Differentiation else if(c[i] != 0) // Only evaluate function if it will actually be used. { points[i] = f(x + (i - Center) * h); - Evaluations++; + Evaluations++; } } @@ -266,7 +266,7 @@ namespace MathNet.Numerics.Differentiation { x[parameterIndex] = xi + (i - Center) * h; points[i] = f(x); - Evaluations++; + Evaluations++; } } @@ -351,7 +351,7 @@ namespace MathNet.Numerics.Differentiation if (order == 1) return EvaluatePartialDerivative(f, x, parameterIndex[0], order, currentValue); - + int reducedOrder = order - 1; var reducedParameterIndex = new int[reducedOrder]; Array.Copy(parameterIndex, 0, reducedParameterIndex, 0, reducedOrder); @@ -367,7 +367,7 @@ namespace MathNet.Numerics.Differentiation points[i] = EvaluateMixedPartialDerivative(f, x, reducedParameterIndex, reducedOrder); } - // restore original value + // restore original value x[currentParameterIndex] = xi; // This will always be to the first order @@ -379,7 +379,7 @@ namespace MathNet.Numerics.Differentiation /// /// /// This function recursively uses to evaluate mixed partial derivative. - /// Therefore, it is more efficient to call for higher order derivatives of + /// Therefore, it is more efficient to call for higher order derivatives of /// a single independent variable. /// /// Multivariate function array handle. @@ -437,21 +437,6 @@ namespace MathNet.Numerics.Differentiation Evaluations = 0; } - /// - /// Calculates machine epsilon - the smallest number that can be added to 1, yeilding a results different than 1. - /// This is also known as roundoff error. - /// - /// Machine epislon - public static double CalculateMachineEpsilon() - { - double eps = 1; - - while ((1.0d + (eps / 2.0d)) > 1.0d) - eps /= 2.0d; - - return eps; - } - private double[] CalculateStepSize(int points, double[] x, double order) { var h = new double[x.Length]; @@ -463,12 +448,12 @@ namespace MathNet.Numerics.Differentiation private double CalculateStepSize(int points, double x, double order) { - // Step size relative to function input parameter + // Step size relative to function input parameter if (StepType == StepType.RelativeX) { StepSize = BaseStepSize*(1 + Math.Abs(x)); } - // Step size relative to function input parameter and order + // Step size relative to function input parameter and order else if (StepType == StepType.Relative) { var accuracy = points - order; diff --git a/src/Numerics/Precision.cs b/src/Numerics/Precision.cs index 2b42b104..dfe75eb2 100644 --- a/src/Numerics/Precision.cs +++ b/src/Numerics/Precision.cs @@ -4,7 +4,7 @@ // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com // -// Copyright (c) 2009-2013 Math.NET +// Copyright (c) 2009-2015 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation @@ -91,15 +91,43 @@ namespace MathNet.Numerics const int SingleWidth = 24; /// - /// The maximum relative precision of of double-precision floating numbers (64 bit) + /// Standard epsilon, the maximum relative precision of IEEE 754 double-precision floating numbers (64 bit). + /// According to the definition of Prof. Demmel and used in LAPACK and Scilab. /// public static readonly double DoublePrecision = Math.Pow(2, -DoubleWidth); /// - /// The maximum relative precision of of single-precision floating numbers (32 bit) + /// Standard epsilon, the maximum relative precision of IEEE 754 double-precision floating numbers (64 bit). + /// According to the definition of Prof. Higham and used in the ISO C standard and MATLAB. + /// + public static readonly double PositiveDoublePrecision = 2*DoublePrecision; + + /// + /// Standard epsilon, the maximum relative precision of IEEE 754 single-precision floating numbers (32 bit). + /// According to the definition of Prof. Demmel and used in LAPACK and Scilab. /// public static readonly double SinglePrecision = Math.Pow(2, -SingleWidth); + /// + /// Standard epsilon, the maximum relative precision of IEEE 754 single-precision floating numbers (32 bit). + /// According to the definition of Prof. Higham and used in the ISO C standard and MATLAB. + /// + public static readonly double PositiveSinglePrecision = 2*SinglePrecision; + + /// + /// Actual machine epsilon, the smallest number that can be subtracted from 1, yielding a results different than 1. + /// This is also known as unit roundoff error. According to the definition of Prof. Demmel. + /// On a standard machine this is equivalent to `DoublePrecision`. + /// + public static readonly double MachineEpsilon = MeasureMachineEpsilon(); + + /// + /// Actual machine epsilon, the smallest number that can be added to 1, yielding a results different than 1. + /// This is also known as unit roundoff error. According to the definition of Prof. Higham. + /// On a standard machine this is equivalent to `PositiveDoublePrecision`. + /// + public static readonly double PositiveMachineEpsilon = MeasurePositiveMachineEpsilon(); + /// /// The number of significant decimal places of double-precision floating numbers (64 bit). /// @@ -761,7 +789,7 @@ namespace MathNet.Numerics } /// - /// Converts a float valut to a bit array stored in an int. + /// Converts a float value to a bit array stored in an int. /// /// The value to convert. /// The bit array. @@ -770,6 +798,36 @@ namespace MathNet.Numerics return BitConverter.ToInt32(BitConverter.GetBytes(value), 0); } + /// + /// Calculates the actual positive double precision machine epsilon - the smallest number that can be added to 1, yielding a results different than 1. + /// This is also known as unit roundoff error. According to the definition of Prof. Demmel. + /// + /// Positive Machine epsilon + static double MeasureMachineEpsilon() + { + double eps = 1.0d; + + while ((1.0d - (eps / 2.0d)) < 1.0d) + eps /= 2.0d; + + return eps; + } + + /// + /// Calculates the actual positive double precision machine epsilon - the smallest number that can be added to 1, yielding a results different than 1. + /// This is also known as unit roundoff error. According to the definition of Prof. Higham. + /// + /// Machine epsilon + static double MeasurePositiveMachineEpsilon() + { + double eps = 1.0d; + + while ((1.0d + (eps / 2.0d)) > 1.0d) + eps /= 2.0d; + + return eps; + } + #if PORTABLE static long DoubleToInt64Bits(double value) { diff --git a/src/Numerics/RootFinding/Brent.cs b/src/Numerics/RootFinding/Brent.cs index f94b748f..833e535b 100644 --- a/src/Numerics/RootFinding/Brent.cs +++ b/src/Numerics/RootFinding/Brent.cs @@ -119,7 +119,7 @@ namespace MathNet.Numerics.RootFinding } // convergence check - double xAcc1 = 2.0*Precision.DoublePrecision*Math.Abs(root) + 0.5*accuracy; + double xAcc1 = Precision.PositiveDoublePrecision*Math.Abs(root) + 0.5*accuracy; double xMidOld = xMid; xMid = (upperBound - root)/2.0;