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;