diff --git a/build/build.proj b/build/build.proj index 9dca1773..9aa578a2 100644 --- a/build/build.proj +++ b/build/build.proj @@ -47,6 +47,11 @@ + + + + + diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebra.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebra.cs index dab836da..e87fdb33 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebra.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebra.cs @@ -55,12 +55,12 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLength); } - if (alpha.AlmostZero()) + if (alpha == 0.0) { return; } - if (alpha.AlmostEqual(1.0)) + if (alpha == 1.0) { Parallel.For(0, y.Length, i => y[i] += x[i]); } @@ -70,12 +70,19 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } } + /// + /// Scales an array. Can be used to scale a vector and a matrix. + /// + /// The scalar. + /// The values to scale. + /// This is equivalent to the SCAL BLAS routine. public void ScaleArray(double alpha, double[] x) { - if (alpha.AlmostEqual(1.0)) + if (alpha == 1.0) { return; } + Parallel.For(0, x.Length, i => x[i] = alpha * x[i]); } } diff --git a/src/Numerics/Complex.cs b/src/Numerics/Complex.cs index 8430127c..895d7420 100644 --- a/src/Numerics/Complex.cs +++ b/src/Numerics/Complex.cs @@ -206,12 +206,12 @@ namespace MathNet.Numerics } /// - /// Gets a value indicating whether whether the Complex is zero. + /// Gets a value indicating whether the Complex is zero. /// /// true if this instance is zero; otherwise, false. public bool IsZero { - get { return _real.AlmostZero() && _imag.AlmostZero(); } + get { return _real == 0.0 && _imag == 0.0; } } /// @@ -220,7 +220,7 @@ namespace MathNet.Numerics /// true if this instance is one; otherwise, false. public bool IsOne { - get { return _real.AlmostEqual(1.0) && _imag.AlmostZero(); } + get { return _real == 1.0 && _imag == 0.0; } } /// @@ -229,7 +229,7 @@ namespace MathNet.Numerics /// true if this instance is I; otherwise, false. public bool IsI { - get { return _real.AlmostZero() && _imag.AlmostEqual(1.0); } + get { return _real == 0.0 && _imag == 1.0; } } /// @@ -267,7 +267,7 @@ namespace MathNet.Numerics /// true if this instance is a real number; otherwise, false. public bool IsReal { - get { return _imag.AlmostZero(); } + get { return _imag == 0.0; } } /// @@ -278,7 +278,7 @@ namespace MathNet.Numerics /// public bool IsRealNonNegative { - get { return _imag.AlmostZero() && _real >= 0; } + get { return _imag == 0.0 && _real >= 0; } } /// @@ -369,7 +369,7 @@ namespace MathNet.Numerics // don't replace this with "Modulus"! var mod = SpecialFunctions.Hypotenuse(_real, _imag); - if (mod.AlmostZero()) + if (mod == 0.0) { return Zero; } @@ -438,7 +438,7 @@ namespace MathNet.Numerics if (exponent.Real < 0) { - if (exponent.Imaginary.AlmostZero()) + if (exponent.Imaginary == 0.0) { return new Complex(double.PositiveInfinity, 0.0); } @@ -652,14 +652,14 @@ namespace MathNet.Numerics var ret = new StringBuilder(); - if (!_real.AlmostZero()) + if (_real != 0.0) { ret.Append(_real.ToString(format, formatProvider)); } - if (!_imag.AlmostZero()) + if (_imag != 0.0) { - if (!_real.AlmostZero()) + if (_real != 0.0) { if (_imag < 0) { @@ -676,7 +676,7 @@ namespace MathNet.Numerics if (ret.Length == 0) { - ret.Append((0.0).ToString(format, formatProvider)); + ret.Append(0.0.ToString(format, formatProvider)); } return ret.ToString(); @@ -910,7 +910,7 @@ namespace MathNet.Numerics /// The divisor. public static Complex operator /(Complex dividend, double divisor) { - if (divisor.AlmostZero()) + if (divisor == 0.0) { return Infinity; } diff --git a/src/Numerics/Interpolation/Algorithms/AkimaSplineInterpolation.cs b/src/Numerics/Interpolation/Algorithms/AkimaSplineInterpolation.cs index 6d32cd5f..238239bd 100644 --- a/src/Numerics/Interpolation/Algorithms/AkimaSplineInterpolation.cs +++ b/src/Numerics/Interpolation/Algorithms/AkimaSplineInterpolation.cs @@ -156,7 +156,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms for (int i = 2; i < derivatives.Length - 2; i++) { derivatives[i] = - weights[i - 1].AlmostZero() && weights[i + 1].AlmostZero() + weights[i - 1].AlmostEqual(0.0) && weights[i + 1].AlmostEqual(0.0) ? (((samplePoints[i + 1] - samplePoints[i]) * differences[i - 1]) + ((samplePoints[i] - samplePoints[i - 1]) * differences[i])) / (samplePoints[i + 1] - samplePoints[i - 1]) diff --git a/src/Numerics/Interpolation/Algorithms/BulirschStoerRationalInterpolation.cs b/src/Numerics/Interpolation/Algorithms/BulirschStoerRationalInterpolation.cs index 51f5552c..308723cf 100644 --- a/src/Numerics/Interpolation/Algorithms/BulirschStoerRationalInterpolation.cs +++ b/src/Numerics/Interpolation/Algorithms/BulirschStoerRationalInterpolation.cs @@ -141,7 +141,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms for (int i = 0; i < n; i++) { double distance = Math.Abs(t - _points[i]); - if (distance.AlmostZero()) + if (distance.AlmostEqual(0.0)) { return _values[i]; } @@ -166,7 +166,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms double ho = (_points[i] - t) * d[i] / hp; double den = ho - c[i + 1]; - if (den.AlmostZero()) + if (den.AlmostEqual(0.0)) { return double.NaN; // zero-div, singularity } diff --git a/src/Numerics/LinearAlgebra/Double/DenseVector.cs b/src/Numerics/LinearAlgebra/Double/DenseVector.cs index ffa5e643..ebf0c501 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseVector.cs @@ -226,7 +226,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The scalar to add. public override void Add(double scalar) { - if (scalar.AlmostZero()) + if (scalar == 0.0) { return; } @@ -379,7 +379,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The scalar to subtract. public override void Subtract(double scalar) { - if (scalar.AlmostZero()) + if (scalar == 0.0) { return; } diff --git a/src/Numerics/LinearAlgebra/Double/Vector.cs b/src/Numerics/LinearAlgebra/Double/Vector.cs index 43ad3f46..a90cfa26 100644 --- a/src/Numerics/LinearAlgebra/Double/Vector.cs +++ b/src/Numerics/LinearAlgebra/Double/Vector.cs @@ -31,7 +31,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double using System; using System.Collections; using System.Collections.Generic; - using System.Diagnostics.CodeAnalysis; using System.Text; using Properties; @@ -117,7 +116,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The scalar to add. public virtual void Add(double scalar) { - if (scalar.AlmostZero()) + if (scalar == 0.0) { return; } @@ -223,7 +222,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The scalar to subtract. public virtual void Subtract(double scalar) { - if (scalar.AlmostZero()) + if (scalar == 0.0) { return; } @@ -667,7 +666,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double { var norm = Norm(); var clone = Clone(); - if (norm.AlmostZero()) + if (norm == 0.0) { return clone; } diff --git a/src/Numerics/Precision.cs b/src/Numerics/Precision.cs index b4c1adf7..3e2c6726 100644 --- a/src/Numerics/Precision.cs +++ b/src/Numerics/Precision.cs @@ -71,33 +71,7 @@ namespace MathNet.Numerics private const int SinglePrecision = 24; #endregion - - #region Properties - /// - /// Gets the maximum relative precision of a double. - /// - /// The maximum relative precision of a double. - public static double DoubleMachinePrecision - { - get - { - return _doubleMachinePrecision; - } - } - - /// - /// Gets the maximum relative precision of a single. - /// - /// The maximum relative precision of a single. - public static double SingleMachinePrecision - { - get - { - return _singleMachinePrecision; - } - } - #endregion - + #region Fields /// @@ -125,6 +99,32 @@ namespace MathNet.Numerics #endregion + #region Properties + /// + /// Gets the maximum relative precision of a double. + /// + /// The maximum relative precision of a double. + public static double DoubleMachinePrecision + { + get + { + return _doubleMachinePrecision; + } + } + + /// + /// Gets the maximum relative precision of a single. + /// + /// The maximum relative precision of a single. + public static double SingleMachinePrecision + { + get + { + return _singleMachinePrecision; + } + } + #endregion + /// /// Initializes static members of the Precision class. /// @@ -440,7 +440,7 @@ namespace MathNet.Numerics /// Forces small numbers near zero to zero. /// /// The real number to coerce to zero, if it is almost zero. - /// Zero if || is smaller than 10*2^(-52) = 0.22e-14, otherwise. + /// Zero if || is smaller than 2^(-53) = 1.11e-16, otherwise. public static double CoerceZero(this double a) { return CoerceZero(a, _doubleMachinePrecision); @@ -682,65 +682,6 @@ namespace MathNet.Numerics return (a >= b) ? (ulong)(intA - intB) : (ulong)(intB - intA); } - /// - /// True if the given number is almost equal to zero, according to the specified absolute accuracy. - /// - /// The real number to check for being almost zero. - /// The maximum number of floating point values between the two values. Must be 1 or larger. - /// - /// True if || is less than steps from zero, False otherwise. - /// - public static bool AlmostZero(this double a, long maxNumbersBetween) - { - if (maxNumbersBetween < 0) - { - throw new ArgumentOutOfRangeException("maxNumbersBetween"); - } - - if (double.IsNaN(a) || double.IsInfinity(a)) - { - return false; - } - - ulong realNumbersBetween = NumbersBetween(0.0, a); - return realNumbersBetween <= (ulong)maxNumbersBetween; - } - - /// - /// True if the given number is almost equal to zero, according to the specified absolute accuracy. - /// - /// The real number to check for being almost zero. - /// The absolute threshold for to consider it as zero. - /// - /// True if || is smaller than , False otherwise. - /// - public static bool AlmostZero(this double a, double maximumAbsoluteError) - { - if (maximumAbsoluteError < 0) - { - throw new ArgumentOutOfRangeException("maximumAbsoluteError"); - } - - if (double.IsNaN(a) || double.IsInfinity(a)) - { - return false; - } - - return Math.Abs(a) < maximumAbsoluteError; - } - - /// - /// True if the given number is almost equal to zero. - /// - /// The real number to check for being almost zero. - /// - /// True if || is smaller than 10*2^(-52) = 0.22e-14, False otherwise. - /// - public static bool AlmostZero(this double a) - { - return AlmostZero(a, _doubleMachinePrecision); - } - /// /// Checks whether two real numbers are almost equal. /// @@ -749,7 +690,8 @@ namespace MathNet.Numerics /// true if the two values differ by no more than 10 * 2^(-52); false otherwise. public static bool AlmostEqual(this double a, double b) { - return AlmostEqualWithError(a, b, a - b, _defaultRelativeAccuracy); + double diff = a - b; + return AlmostEqualWithError(a, b, diff, _defaultRelativeAccuracy); } /// @@ -762,7 +704,8 @@ namespace MathNet.Numerics public static bool AlmostEqual(this T a, T b) where T : IPrecisionSupport { - return AlmostEqualWithError(a.Norm(), b.Norm(), a.NormOfDifference(b), _defaultRelativeAccuracy); + double diff = a.NormOfDifference(b); + return AlmostEqualWithError(a.Norm(), b.Norm(), diff, _defaultRelativeAccuracy); } /// @@ -912,7 +855,7 @@ namespace MathNet.Numerics return false; } - if (AlmostZero(a) || AlmostZero(b)) + if (Math.Abs(a) < _doubleMachinePrecision || Math.Abs(b) < _doubleMachinePrecision) { return AlmostEqualWithAbsoluteError(a, b, diff, maximumError); } @@ -1036,6 +979,11 @@ namespace MathNet.Numerics { return a == b; } + + if (Math.Abs(a) < _doubleMachinePrecision || Math.Abs(b) < _doubleMachinePrecision) + { + return AlmostEqualWithAbsoluteDecimalPlaces(a, b, decimalPlaces); + } // If both numbers are equal, get out now. This should remove the possibility of both numbers being zero // and any problems associated with that. @@ -1044,11 +992,6 @@ namespace MathNet.Numerics return true; } - if (AlmostZero(a) || AlmostZero(b)) - { - return AlmostEqualWithAbsoluteDecimalPlaces(a, b, decimalPlaces); - } - return AlmostEqualWithRelativeDecimalPlaces(a, b, decimalPlaces); } @@ -1385,5 +1328,48 @@ namespace MathNet.Numerics // larger than the second. return a.CompareTo(b); } + + /// + /// Evaluates the minimum distance to the next distinguishable number near the argument value. + /// + /// The value used to determine the minimum distance. + /// + /// Relative Epsilon (positive double or NaN). + /// + /// Evaluates the negative epsilon. The more common positive epsilon is equal to two times this negative epsilon. + /// + public static double EpsilonOf(this double value) + { + if (double.IsInfinity(value) || double.IsNaN(value)) + { + return double.NaN; + } + + long signed64 = BitConverter.DoubleToInt64Bits(value); + if (signed64 == 0) + { + signed64++; + return BitConverter.Int64BitsToDouble(signed64) - value; + } + + if (signed64-- < 0) + { + return BitConverter.Int64BitsToDouble(signed64) - value; + } + + return value - BitConverter.Int64BitsToDouble(signed64); + } + + /// + /// Evaluates the minimum distance to the next distinguishable number near the argument value. + /// + /// The value used to determine the minimum distance. + /// Relative Epsilon (positive double or NaN) + /// Evaluates the positive epsilon. See also + /// + public static double PositiveEpsilonOf(this double value) + { + return 2 * EpsilonOf(value); + } } } \ No newline at end of file diff --git a/src/Numerics/SpecialFunctions/Stability.cs b/src/Numerics/SpecialFunctions/Stability.cs index 7bd0acaa..f4179e08 100644 --- a/src/Numerics/SpecialFunctions/Stability.cs +++ b/src/Numerics/SpecialFunctions/Stability.cs @@ -45,7 +45,7 @@ namespace MathNet.Numerics return Math.Exp(power) - 1.0; } - if (x < Precision.DoubleMachinePrecision) + if (x < x.PositiveEpsilonOf()) { return x; } diff --git a/src/UnitTests/PrecisionTest.cs b/src/UnitTests/PrecisionTest.cs index 188496ea..3dc5cac2 100644 --- a/src/UnitTests/PrecisionTest.cs +++ b/src/UnitTests/PrecisionTest.cs @@ -184,6 +184,8 @@ namespace MathNet.Numerics.UnitTests public void CoerceZero() { Assert.AreEqual(0.0, Precision.CoerceZero(0d)); + Console.WriteLine(0.0.EpsilonOf()); + Console.WriteLine(Precision.Increment(0.0)); Assert.AreEqual(0.0, Precision.CoerceZero(Precision.Increment(0.0))); Assert.AreEqual(0.0, Precision.CoerceZero(Precision.Decrement(0.0))); @@ -450,31 +452,6 @@ namespace MathNet.Numerics.UnitTests Assert.AreEqual(3, Precision.NumbersBetween(double.Epsilon, -2 * double.Epsilon)); } - // AlmostZero - [Test] - public void AlmostZeroWithMaxNumbersBetween() - { - Assert.IsTrue(Precision.AlmostZero(0, 1)); - Assert.IsTrue(Precision.AlmostZero(1 * double.Epsilon, 1)); - Assert.IsTrue(Precision.AlmostZero(10 * double.Epsilon, 10)); - - Assert.IsFalse(Precision.AlmostZero(double.NegativeInfinity, 1)); - Assert.IsFalse(Precision.AlmostZero(double.PositiveInfinity, 1)); - Assert.IsFalse(Precision.AlmostZero(double.NaN, 1)); - } - - [Test] - public void AlmostZeroWithTolerance() - { - Assert.IsTrue(Precision.AlmostZero(0)); - Assert.IsTrue(Precision.AlmostZero(1 * double.Epsilon, 2 * double.Epsilon)); - Assert.IsTrue(Precision.AlmostZero(10 * double.Epsilon, 20 * double.Epsilon)); - - Assert.IsFalse(Precision.AlmostZero(double.NegativeInfinity, 1.0)); - Assert.IsFalse(Precision.AlmostZero(double.PositiveInfinity, 1.0)); - Assert.IsFalse(Precision.AlmostZero(double.NaN, 1.0)); - } - [Test] [ExpectedException(typeof(ArgumentOutOfRangeException))] public void AlmostEqualWithMaxNumbersBetweenWithLessThanOneNumber() diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 82cdf056..8534bd45 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -121,6 +121,234 @@ + + data\NIST\AtmWtAgt.dat + PreserveNewest + + + data\NIST\Bennett5.dat + PreserveNewest + + + data\NIST\BoxBOD.dat + PreserveNewest + + + data\NIST\Chwirut1.dat + PreserveNewest + + + data\NIST\Chwirut2.dat + PreserveNewest + + + data\NIST\DanWood.dat + PreserveNewest + + + data\NIST\Eckerle4.dat + PreserveNewest + + + data\NIST\ENSO.dat + PreserveNewest + + + data\NIST\Filip.dat + PreserveNewest + + + data\NIST\Gauss1.dat + PreserveNewest + + + data\NIST\Gauss2.dat + PreserveNewest + + + data\NIST\Gauss3.dat + PreserveNewest + + + data\NIST\Hahn1.dat + PreserveNewest + + + data\NIST\Kirby2.dat + PreserveNewest + + + data\NIST\Lanczos1.dat + PreserveNewest + + + data\NIST\Lanczos2.dat + PreserveNewest + + + data\NIST\Lanczos3.dat + PreserveNewest + + + data\NIST\Lew.dat + PreserveNewest + + + data\NIST\Longley.dat + PreserveNewest + + + data\NIST\Lottery.dat + PreserveNewest + + + data\NIST\Mavro.dat + PreserveNewest + + + data\NIST\MGH09.dat + PreserveNewest + + + data\NIST\MGH10.dat + PreserveNewest + + + data\NIST\MGH17.dat + PreserveNewest + + + data\NIST\Michelso.dat + PreserveNewest + + + data\NIST\Misra1a.dat + PreserveNewest + + + data\NIST\Misra1b.dat + PreserveNewest + + + data\NIST\Misra1c.dat + PreserveNewest + + + data\NIST\Misra1d.dat + PreserveNewest + + + data\NIST\Nelson.dat + PreserveNewest + + + data\NIST\NoInt1.dat + PreserveNewest + + + data\NIST\NoInt2.dat + PreserveNewest + + + data\NIST\Norris.dat + PreserveNewest + + + data\NIST\NumAcc1.dat + PreserveNewest + + + data\NIST\NumAcc2.dat + PreserveNewest + + + data\NIST\NumAcc3.dat + PreserveNewest + + + data\NIST\NumAcc4.dat + PreserveNewest + + + data\NIST\Pontius.dat + PreserveNewest + + + data\NIST\Rat42.dat + PreserveNewest + + + data\NIST\Rat43.dat + PreserveNewest + + + data\NIST\Roszman1.dat + PreserveNewest + + + data\NIST\SiRstvt.dat + PreserveNewest + + + data\NIST\SmLs01t.dat + PreserveNewest + + + data\NIST\SmLs02t.dat + PreserveNewest + + + data\NIST\SmLs03t.dat + PreserveNewest + + + data\NIST\SmLs04t.dat + PreserveNewest + + + data\NIST\SmLs05t.dat + PreserveNewest + + + data\NIST\SmLs06t.dat + PreserveNewest + + + data\NIST\SmLs07t.dat + PreserveNewest + + + data\NIST\SmLs08t.dat + PreserveNewest + + + data\NIST\SmLs09t.dat + PreserveNewest + + + data\NIST\Thurber.dat + PreserveNewest + + + data\NIST\Wampler1.dat + PreserveNewest + + + data\NIST\Wampler2.dat + PreserveNewest + + + data\NIST\Wampler3.dat + PreserveNewest + + + data\NIST\Wampler4.dat + PreserveNewest + + + data\NIST\Wampler5.dat + PreserveNewest + MathNet.Numerics.snk