diff --git a/src/Managed.UnitTests/AssertHelpers.cs b/src/Managed.UnitTests/AssertHelpers.cs index ae6b4cf0..bef96af9 100644 --- a/src/Managed.UnitTests/AssertHelpers.cs +++ b/src/Managed.UnitTests/AssertHelpers.cs @@ -26,8 +26,10 @@ // OTHER DEALINGS IN THE SOFTWARE. // + namespace MathNet.Numerics.UnitTests { + using System.Collections.Generic; using MbUnit.Framework; /// @@ -71,5 +73,44 @@ namespace MathNet.Numerics.UnitTests Assert.Fail("Imaginary components are not equal within {0} places. Expected:{1}; Actual:{2}", decimalPlaces, expected.Imaginary, actual.Imaginary); } } + + /// + /// Asserts that the expected value and the actual value are equal up to a certain + /// maximum error. + /// + /// The type of the structures. Must implement + /// . + /// The expected value. + /// The actual value. + /// The accuracy required for being almost equal. + public static void AlmostEqual(T expected, T actual, double maximumError) + where T : IPrecisionSupport + { + if (!actual.AlmostEqualWithError(expected, maximumError)) + { + Assert.Fail("Not equal within a maximum error {0}. Expected:{1}; Actual:{2}", maximumError, expected, actual); + } + } + + /// + /// Asserts that the expected value and the actual value are equal up to a certain + /// maximum error. + /// + /// The type of the structures. Must implement + /// . + /// The expected value list. + /// The actual value list. + /// The accuracy required for being almost equal. + public static void AlmostEqualList(IList expected, IList actual, double maximumError) + where T : IPrecisionSupport + { + for (int i = 0; i < expected.Count; i++) + { + if (!actual[i].AlmostEqualWithError(expected[i], maximumError)) + { + Assert.Fail("Not equal within a maximum error {0}. Expected:{1}; Actual:{2}", maximumError, expected[i], actual[i]); + } + } + } } } diff --git a/src/Managed.UnitTests/PrecisionTest.cs b/src/Managed.UnitTests/PrecisionTest.cs index 4b038ad1..188496ea 100644 --- a/src/Managed.UnitTests/PrecisionTest.cs +++ b/src/Managed.UnitTests/PrecisionTest.cs @@ -496,53 +496,53 @@ namespace MathNet.Numerics.UnitTests public void AlmostEqualWithRelativeError() { // compare zero and negative zero - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(0.0, -0.0, 1e-5)); - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(0.0, -0.0, 1e-15)); + Assert.IsTrue(Precision.AlmostEqualWithError(0.0, -0.0, 1e-5)); + Assert.IsTrue(Precision.AlmostEqualWithError(0.0, -0.0, 1e-15)); // compare two nearby numbers - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(1.0, 1.0 + 3 * _doublePrecision, 1e-15)); - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(1.0, 1.0 + _doublePrecision, 1e-15)); - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(1.0, 1.0 + 1e-16, 1e-15)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(1.0, 1.0 + 1e-15, 1e-15)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(1.0, 1.0 + 1e-14, 1e-15)); + Assert.IsTrue(Precision.AlmostEqualWithError(1.0, 1.0 + 3 * _doublePrecision, 1e-15)); + Assert.IsTrue(Precision.AlmostEqualWithError(1.0, 1.0 + _doublePrecision, 1e-15)); + Assert.IsTrue(Precision.AlmostEqualWithError(1.0, 1.0 + 1e-16, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(1.0, 1.0 + 1e-15, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(1.0, 1.0 + 1e-14, 1e-15)); // compare with the two numbers reversed in compare order - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(1.0 + 3 * _doublePrecision, 1.0, 1e-15)); - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(1.0 + _doublePrecision, 1.0, 1e-15)); - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(1.0 + 1e-16, 1.0, 1e-15)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(1.0 + 1e-15, 1.0, 1e-15)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(1.0 + 1e-14, 1.0, 1e-15)); + Assert.IsTrue(Precision.AlmostEqualWithError(1.0 + 3 * _doublePrecision, 1.0, 1e-15)); + Assert.IsTrue(Precision.AlmostEqualWithError(1.0 + _doublePrecision, 1.0, 1e-15)); + Assert.IsTrue(Precision.AlmostEqualWithError(1.0 + 1e-16, 1.0, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(1.0 + 1e-15, 1.0, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(1.0 + 1e-14, 1.0, 1e-15)); // compare different numbers - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(2.0, 1.0, 1e-15)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(1.0, 2.0, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(2.0, 1.0, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(1.0, 2.0, 1e-15)); // compare different numbers with large tolerance - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(2.0, 1.0, 1e-5)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(1.0, 2.0, 1e-5)); - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(2.0, 1.0, 1e+1)); - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(1.0, 2.0, 1e+1)); + Assert.IsFalse(Precision.AlmostEqualWithError(2.0, 1.0, 1e-5)); + Assert.IsFalse(Precision.AlmostEqualWithError(1.0, 2.0, 1e-5)); + Assert.IsTrue(Precision.AlmostEqualWithError(2.0, 1.0, 1e+1)); + Assert.IsTrue(Precision.AlmostEqualWithError(1.0, 2.0, 1e+1)); // compare inf & inf - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(double.PositiveInfinity, double.PositiveInfinity, 1e-15)); - Assert.IsTrue(Precision.AlmostEqualWithRelativeError(double.NegativeInfinity, double.NegativeInfinity, 1e-15)); + Assert.IsTrue(Precision.AlmostEqualWithError(double.PositiveInfinity, double.PositiveInfinity, 1e-15)); + Assert.IsTrue(Precision.AlmostEqualWithError(double.NegativeInfinity, double.NegativeInfinity, 1e-15)); // compare -inf and inf - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(double.PositiveInfinity, double.NegativeInfinity, 1e-15)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(double.NegativeInfinity, double.PositiveInfinity, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(double.PositiveInfinity, double.NegativeInfinity, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(double.NegativeInfinity, double.PositiveInfinity, 1e-15)); // compare inf and non-inf - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(double.PositiveInfinity, 1.0, 1e-15)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(1.0, double.PositiveInfinity, 1e-15)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(double.NegativeInfinity, 1.0, 1e-15)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(1.0, double.NegativeInfinity, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(double.PositiveInfinity, 1.0, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(1.0, double.PositiveInfinity, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(double.NegativeInfinity, 1.0, 1e-15)); + Assert.IsFalse(Precision.AlmostEqualWithError(1.0, double.NegativeInfinity, 1e-15)); // compare tiny numbers with opposite signs - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(1e-12, -1e-12, 1e-14)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(-1e-12, 1e-12, 1e-14)); + Assert.IsFalse(Precision.AlmostEqualWithError(1e-12, -1e-12, 1e-14)); + Assert.IsFalse(Precision.AlmostEqualWithError(-1e-12, 1e-12, 1e-14)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(-2.0, 2.0, 1e-14)); - Assert.IsFalse(Precision.AlmostEqualWithRelativeError(2.0, -2.0, 1e-14)); + Assert.IsFalse(Precision.AlmostEqualWithError(-2.0, 2.0, 1e-14)); + Assert.IsFalse(Precision.AlmostEqualWithError(2.0, -2.0, 1e-14)); } [Test] diff --git a/src/Managed/Complex.cs b/src/Managed/Complex.cs index 994304a4..2ff5fc02 100644 --- a/src/Managed/Complex.cs +++ b/src/Managed/Complex.cs @@ -66,7 +66,7 @@ namespace MathNet.Numerics /// [Serializable] [StructLayout(LayoutKind.Sequential)] - public struct Complex : IFormattable, IEquatable + public struct Complex : IFormattable, IEquatable, IPrecisionSupport { #region fields @@ -1061,5 +1061,51 @@ namespace MathNet.Numerics } #endregion + + #region Trigonometric Functions + + /// + /// Trigonometric Sine (sin, Sinus) of this Complex. + /// + /// + /// The sine of the complex number. + /// + public Complex Sine() + { + if (this.IsReal) + { + return new Complex(Math.Sin(this._real), 0.0); + } + + return new Complex( + Math.Sin(this._real) * Math.Cosh(this._imag), Math.Cos(this._real) * Math.Sinh(this._imag)); + } + + #endregion + + #region IPrecisionSupport + + /// + /// Returns a Norm of a value of this type, which is appropriate for measuring how + /// close this value is to zero. + /// + /// A norm of this value. + double IPrecisionSupport.Norm() + { + return ModulusSquared; + } + + /// + /// Returns a Norm of the difference of two values of this type, which is + /// appropriate for measuring how close together these two values are. + /// + /// The value to compare with. + /// A norm of the difference between this and the other value. + double IPrecisionSupport.NormOfDifference(Complex otherValue) + { + return (this - otherValue).ModulusSquared; + } + + #endregion } } \ No newline at end of file diff --git a/src/Managed/IPrecisionSupport.cs b/src/Managed/IPrecisionSupport.cs new file mode 100644 index 00000000..2cdfc188 --- /dev/null +++ b/src/Managed/IPrecisionSupport.cs @@ -0,0 +1,52 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics +{ + /// + /// Support Interface for Precision Operations (like AlmostEquals). + /// + /// Type of the implementing class. + public interface IPrecisionSupport + { + /// + /// Returns a Norm of a value of this type, which is appropriate for measuring how + /// close this value is to zero. + /// + /// A norm of this value. + double Norm(); + + /// + /// Returns a Norm of the difference of two values of this type, which is + /// appropriate for measuring how close together these two values are. + /// + /// The value to compare with. + /// A norm of the difference between this and the other value. + double NormOfDifference(T otherValue); + } +} diff --git a/src/Managed/Integration/Algorithms/NewtonCotesTrapeziumRule.cs b/src/Managed/Integration/Algorithms/NewtonCotesTrapeziumRule.cs index d2b848c9..a3b5cc91 100644 --- a/src/Managed/Integration/Algorithms/NewtonCotesTrapeziumRule.cs +++ b/src/Managed/Integration/Algorithms/NewtonCotesTrapeziumRule.cs @@ -94,13 +94,13 @@ namespace MathNet.Numerics.Integration.Algorithms /// The analytic smooth function to integrate. /// Where the interval starts, inclusive and finite. /// Where the interval stops, inclusive and finite. - /// The expected relative accuracy of the approximation. + /// The expected accuracy of the approximation. /// Approximation of the finite integral in the given interval. public static double IntegrateAdaptive( Func f, double intervalBegin, double intervalEnd, - double targetRelativeError) + double targetError) { int numberOfPartitions = 1; double step = intervalEnd - intervalBegin; @@ -118,7 +118,7 @@ namespace MathNet.Numerics.Integration.Algorithms step *= 0.5; numberOfPartitions *= 2; - if (sum.AlmostEqualWithRelativeError(midpointsum, targetRelativeError)) + if (sum.AlmostEqualWithError(midpointsum, targetError)) { break; } diff --git a/src/Managed/Managed.csproj b/src/Managed/Managed.csproj index a8b5c64d..7c53c4f8 100644 --- a/src/Managed/Managed.csproj +++ b/src/Managed/Managed.csproj @@ -54,6 +54,7 @@ + diff --git a/src/Managed/Precision.cs b/src/Managed/Precision.cs index a4c93101..718f4b3c 100644 --- a/src/Managed/Precision.cs +++ b/src/Managed/Precision.cs @@ -26,10 +26,11 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System; - namespace MathNet.Numerics { + using System; + using System.Collections.Generic; + /// /// Utilities for working with floating point numbers. /// @@ -102,12 +103,12 @@ namespace MathNet.Numerics /// /// The maximum relative precision of a double /// - private static readonly double _doubleMachinePrecision = System.Math.Pow(BinaryBaseNumber, -DoublePrecision); + private static readonly double _doubleMachinePrecision = Math.Pow(BinaryBaseNumber, -DoublePrecision); /// /// The maximum relative precision of a single /// - private static readonly double _singleMachinePrecision = System.Math.Pow(BinaryBaseNumber, -SinglePrecision); + private static readonly double _singleMachinePrecision = Math.Pow(BinaryBaseNumber, -SinglePrecision); /// /// The number of significant figures that a double-precision floating point has. @@ -129,8 +130,8 @@ namespace MathNet.Numerics /// static Precision() { - _numberOfDecimalPlacesForFloats = (int)System.Math.Ceiling(System.Math.Abs(System.Math.Log10(_singleMachinePrecision))); - _numberOfDecimalPlacesForDoubles = (int)System.Math.Ceiling(System.Math.Abs(System.Math.Log10(_doubleMachinePrecision))); + _numberOfDecimalPlacesForFloats = (int)Math.Ceiling(Math.Abs(Math.Log10(_singleMachinePrecision))); + _numberOfDecimalPlacesForDoubles = (int)Math.Ceiling(Math.Abs(Math.Log10(_doubleMachinePrecision))); } /// @@ -161,7 +162,7 @@ namespace MathNet.Numerics /// Returns the magnitude of the number. /// /// The value. - /// The magnitute of the number. + /// The magnitude of the number. public static int Magnitude(this double value) { // Can't do this with zero because the 10-log of zero doesn't exist. @@ -172,17 +173,17 @@ namespace MathNet.Numerics // Note that we need the absolute value of the input because Log10 doesn't // work for negative numbers (obviously). - double magnitude = System.Math.Log10(System.Math.Abs(value)); + double magnitude = Math.Log10(Math.Abs(value)); // To get the right number we need to know if the value is negative or positive // truncating a positive number will always give use the correct magnitude // truncating a negative number will give us a magnitude that is off by 1 if (magnitude < 0) { - return (int)System.Math.Truncate(magnitude - 1); + return (int)Math.Truncate(magnitude - 1); } - return (int)System.Math.Truncate(magnitude); + return (int)Math.Truncate(magnitude); } /// @@ -198,7 +199,7 @@ namespace MathNet.Numerics } int magnitude = Magnitude(value); - return value * System.Math.Pow(10, -magnitude); + return value * Math.Pow(10, -magnitude); } /// @@ -409,7 +410,7 @@ namespace MathNet.Numerics /// Determines the range of floating point numbers that will match the specified value with the given tolerance. /// /// The value. - /// The ulps difference. + /// The ulps difference. /// The bottom range end. /// The top range end. /// @@ -455,7 +456,7 @@ namespace MathNet.Numerics // Note that long.MinValue has the same bit pattern as // -0.0. Therefore we're working in opposite direction (i.e. add if we want to // go more negative and subtract if we want to go less negative) - if (System.Math.Abs(long.MinValue - intValue) < maxNumbersBetween) + if (Math.Abs(long.MinValue - intValue) < maxNumbersBetween) { // Got underflow, which can be fixed by splitting the calculation into two bits // first get the remainder of the intValue after subtracting it from the long.MinValue @@ -468,7 +469,7 @@ namespace MathNet.Numerics topRangeEnd = BitConverter.Int64BitsToDouble(intValue - maxNumbersBetween); } - if (System.Math.Abs(intValue) < maxNumbersBetween) + if (Math.Abs(intValue) < maxNumbersBetween) { // Underflow, which means we'd have to go further than a long would allow us. // Also we couldn't translate it back to a double, so we'll return -Double.MaxValue @@ -517,7 +518,7 @@ namespace MathNet.Numerics /// always bigger than the value) /// /// The value. - /// The ulps difference. + /// The ulps difference. /// The maximum floating point number which is larger than the given . public static double MaximumMatchingFloatingPointNumber(this double value, long maxNumbersBetween) { @@ -531,7 +532,7 @@ namespace MathNet.Numerics /// always smaller than the value) /// /// The value. - /// The ulps difference. + /// The ulps difference. /// The minimum floating point number which is smaller than the given . public static double MinimumMatchingFloatingPointNumber(this double value, long maxNumbersBetween) { @@ -541,7 +542,7 @@ namespace MathNet.Numerics } /// - /// Determines the range of ulps that will match the specified value with the given tolerance. + /// Determines the range of ulps that will match the specified value with the given tolerance. /// /// The value. /// The relative difference. @@ -588,15 +589,15 @@ namespace MathNet.Numerics // Calculate the ulps for the maximum and minimum values // Note that these can overflow - long max = GetDirectionalLongFromDouble(value + (relativeDifference * System.Math.Abs(value))); - long min = GetDirectionalLongFromDouble(value - (relativeDifference * System.Math.Abs(value))); + long max = GetDirectionalLongFromDouble(value + (relativeDifference * Math.Abs(value))); + long min = GetDirectionalLongFromDouble(value - (relativeDifference * Math.Abs(value))); // Calculate the ulps from the value long intValue = GetDirectionalLongFromDouble(value); // Determine the ranges - topRangeEnd = System.Math.Abs(max - intValue); - bottomRangeEnd = System.Math.Abs(intValue - min); + topRangeEnd = Math.Abs(max - intValue); + bottomRangeEnd = Math.Abs(intValue - min); } /// @@ -708,38 +709,153 @@ 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 AlmostEqualWithRelativeError(a, b, a - b, _defaultRelativeAccuracy); + return AlmostEqualWithError(a, b, a - b, _defaultRelativeAccuracy); + } + + /// + /// Checks whether two structures with precision support are almost equal. + /// + /// The type of the structures. Must implement . + /// The first structure + /// The second structure + /// true if the two values differ by no more than 10 * 2^(-52); false otherwise. + public static bool AlmostEqual(this T a, T b) + where T : IPrecisionSupport + { + return AlmostEqualWithError(a.Norm(), b.Norm(), a.NormOfDifference(b), _defaultRelativeAccuracy); } /// /// Compares two doubles and determines if they are equal within - /// the specified maximum relative error. + /// the specified maximum error. /// /// The first value. /// The second value. - /// The relative accuracy required for being almost equal. + /// The accuracy required for being almost equal. /// /// if both doubles are almost equal up to the - /// specified maximum relative error, otherwise. + /// specified maximum error, otherwise. /// - public static bool AlmostEqualWithRelativeError(this double a, double b, double maximumRelativeError) + public static bool AlmostEqualWithError(this double a, double b, double maximumError) { - return AlmostEqualWithRelativeError(a, b, a - b, maximumRelativeError); + return AlmostEqualWithError(a, b, a - b, maximumError); } /// - /// Compares two doubles and determines if they are equal within - /// the specified maximum relative error. + /// Compares two lists of doubles and determines if they are equal within the + /// specified maximum error. + /// + /// The first value list. + /// The second value list. + /// + /// The accuracy required for being almost equal. + /// + /// + /// if both doubles are almost equal up to the specified + /// maximum error, otherwise. + /// + public static bool AlmostEqualListWithError(this IList a, IList b, double maximumError) + { + if (a == null && b == null) + { + return true; + } + + if (a == null || b == null || a.Count != b.Count) + { + return false; + } + + for (int i = 0; i < a.Count; i++) + { + if (!AlmostEqualWithError(a[i], b[i], a[i] - b[i], maximumError)) + { + return false; + } + } + + return true; + } + + /// + /// Compares two structure with precision support and determines if they are equal + /// within the specified maximum relative error. + /// + /// + /// The type of the structures. Must implement . + /// + /// The first structure. + /// The second structure. + /// + /// The accuracy required for being almost equal. + /// + /// + /// if both doubles are almost equal up to the specified + /// maximum relative error, otherwise. + /// + public static bool AlmostEqualWithError(this T a, T b, double maximumError) + where T : IPrecisionSupport + { + return AlmostEqualWithError(a.Norm(), b.Norm(), a.NormOfDifference(b), maximumError); + } + + /// + /// Compares two lists of structures with precision support and determines if they + /// are equal within the specified maximum error. + /// + /// + /// The type of the structures. Must implement . + /// + /// The first structure list. + /// The second structure list. + /// + /// The accuracy required for being almost equal. + /// + /// + /// if both doubles are almost equal up to the specified + /// maximum error, otherwise. + /// + public static bool AlmostEqualListWithError(this IList a, IList b, double maximumError) + where T : IPrecisionSupport + { + if (a == null && b == null) + { + return true; + } + + if (a == null || b == null || a.Count != b.Count) + { + return false; + } + + for (int i = 0; i < a.Count; i++) + { + if (!AlmostEqualWithError(a[i].Norm(), b[i].Norm(), a[i].NormOfDifference(b[i]), maximumError)) + { + return false; + } + } + + return true; + } + + /// + /// Compares two doubles and determines if they are equal within the specified + /// maximum error. /// /// The first value. /// The second value. - /// The difference of the two values (according to some norm). - /// The relative accuracy required for being almost equal. + /// + /// The difference of the two values (according to some norm). + /// + /// + /// The accuracy required for being almost equal. + /// /// - /// if both doubles are almost equal up to the - /// specified maximum relative error, otherwise. + /// if both doubles are almost equal up to the specified + /// maximum error, otherwise. /// - public static bool AlmostEqualWithRelativeError(this double a, double b, double diff, double maximumRelativeError) + public static bool AlmostEqualWithError(this double a, double b, double diff, double maximumError) { // If A or B are infinity (positive or negative) then // only return true if they are exactly equal to each other - @@ -756,47 +872,66 @@ namespace MathNet.Numerics return false; } - if ((a == 0 && Math.Abs(b) < maximumRelativeError) - || (b == 0 && Math.Abs(a) < maximumRelativeError)) + if (AlmostZero(a) || AlmostZero(b)) { - return true; + return AlmostEqualWithAbsoluteError(a, b, diff, maximumError); } - return Math.Abs(diff) < maximumRelativeError * Math.Max(Math.Abs(a), Math.Abs(b)); + return AlmostEqualWithRelativeError(a, b, diff, maximumError); } /// - /// Compares two doubles and determines if they are equal to within the tolerance or not. Equality comparison is based on the binary representation. + /// Compares two doubles and determines if they are equal within the specified + /// maximum absolute error. /// - /// - /// - /// Determines the 'number' of floating point numbers between two values (i.e. the number of discrete steps - /// between the two numbers) and then checks if that is within the specified tolerance. So if a tolerance - /// of 1 is passed then the result will be true only if the two numbers have the same binary representation - /// OR if they are two adjacent numbers that only differ by one step. - /// - /// - /// The comparison method used is explained in http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm . The article - /// at http://www.extremeoptimization.com/resources/Articles/FPDotNetConceptsAndFormats.aspx explains how to transform the C code to - /// .NET enabled code without using pointers and unsafe code. - /// - /// /// The first value. /// The second value. - /// The maximum number of floating point values between the two values. Must be 1 or larger. - /// if both doubles are equal to each other within the specified tolerance; otherwise . - /// - /// Thrown if is smaller than one. - /// - public static bool AlmostEqual(this double a, double b, long maxNumbersBetween) + /// + /// The difference of the two values (according to some norm). + /// + /// + /// The absolute accuracy required for being almost equal. + /// + /// + /// if both doubles are almost equal up to the specified + /// maximum absolute error, otherwise. + /// + public static bool AlmostEqualWithAbsoluteError(this double a, double b, double diff, double maximumAbsoluteError) { - // Make sure maxNumbersBetween is non-negative and small enough that the - // default NAN won't compare as equal to anything. - if (maxNumbersBetween < 1) + // If A or B are infinity (positive or negative) then + // only return true if they are exactly equal to each other - + // that is, if they are both infinities of the same sign. + if (double.IsInfinity(a) || double.IsInfinity(b)) { - throw new ArgumentOutOfRangeException("maxNumbersBetween"); + return a == b; + } + + // If A or B are a NAN, return false. NANs are equal to nothing, + // not even themselves. + if (double.IsNaN(a) || double.IsNaN(b)) + { + return false; } + return Math.Abs(diff) < maximumAbsoluteError; + } + + /// + /// Compares two doubles and determines if they are equal within the specified + /// maximum relative error. + /// + /// The first value. + /// The second value. + /// The difference of the two values (according to some norm). + /// + /// The relative accuracy required for being + /// almost equal. + /// + /// if both doubles are almost equal up to the specified + /// maximum relative error, otherwise. + /// + public static bool AlmostEqualWithRelativeError(this double a, double b, double diff, double maximumRelativeError) + { // If A or B are infinity (positive or negative) then // only return true if they are exactly equal to each other - // that is, if they are both infinities of the same sign. @@ -812,16 +947,13 @@ namespace MathNet.Numerics return false; } - // Get the first double and convert it to an integer value (by using the binary representation) - long firstUlong = GetDirectionalLongFromDouble(a); - - // Get the second double and convert it to an integer value (by using the binary representation) - long secondUlong = GetDirectionalLongFromDouble(b); + if ((a == 0 && Math.Abs(b) < maximumRelativeError) + || (b == 0 && Math.Abs(a) < maximumRelativeError)) + { + return true; + } - // Now compare the values. - // Note that this comparison can overflow so we'll approach this differently - // Do note that we could overflow this way too. We should probably check that we don't. - return (a > b) ? (secondUlong + maxNumbersBetween >= firstUlong) : (firstUlong + maxNumbersBetween >= secondUlong); + return Math.Abs(diff) < maximumRelativeError * Math.Max(Math.Abs(a), Math.Abs(b)); } /// @@ -831,7 +963,7 @@ namespace MathNet.Numerics /// /// /// The values are equal if the difference between the two numbers is smaller than 10^(-numberOfDecimalPlaces). We divide by - /// two so that we have half the range on each side of the numbers, e.g. if decimalPlaces == 2, then 0.01 will equal between + /// two so that we have half the range on each side of the numbers, e.g. if == 2, then 0.01 will equal between /// 0.005 and 0.015, but not 0.02 and not 0.00 /// /// @@ -886,7 +1018,7 @@ namespace MathNet.Numerics /// /// /// The values are equal if the difference between the two numbers is smaller than 10^(-numberOfDecimalPlaces). We divide by - /// two so that we have half the range on each side of the numbers, e.g. if decimalPlaces == 2, then 0.01 will equal between + /// two so that we have half the range on each side of the numbers, e.g. if == 2, then 0.01 will equal between /// 0.005 and 0.015, but not 0.02 and not 0.00 /// /// @@ -899,13 +1031,13 @@ namespace MathNet.Numerics // If the magnitudes of the two numbers are equal to within one magnitude the numbers could potentially be equal int magnitudeOfFirst = Magnitude(a); int magnitudeOfSecond = Magnitude(b); - if (System.Math.Max(magnitudeOfFirst, magnitudeOfSecond) > (System.Math.Min(magnitudeOfFirst, magnitudeOfSecond) + 1)) + if (Math.Max(magnitudeOfFirst, magnitudeOfSecond) > (Math.Min(magnitudeOfFirst, magnitudeOfSecond) + 1)) { return false; } // Get the power of the number of decimalPlaces - double decimalPlaceMagnitude = System.Math.Pow(10, -(decimalPlaces - 1)); + double decimalPlaceMagnitude = Math.Pow(10, -(decimalPlaces - 1)); // The values are equal if the difference between the two numbers is smaller than // 10^(-numberOfDecimalPlaces). We divide by two so that we have half the range @@ -914,11 +1046,11 @@ namespace MathNet.Numerics double maxDifference = decimalPlaceMagnitude / 2.0; if (a > b) { - return (a * System.Math.Pow(10, -magnitudeOfFirst)) - maxDifference < (b * System.Math.Pow(10, -magnitudeOfFirst)); + return (a * Math.Pow(10, -magnitudeOfFirst)) - maxDifference < (b * Math.Pow(10, -magnitudeOfFirst)); } else { - return (b * System.Math.Pow(10, -magnitudeOfSecond)) - maxDifference < (a * System.Math.Pow(10, -magnitudeOfSecond)); + return (b * Math.Pow(10, -magnitudeOfSecond)) - maxDifference < (a * Math.Pow(10, -magnitudeOfSecond)); } } @@ -929,7 +1061,7 @@ namespace MathNet.Numerics /// /// /// The values are equal if the difference between the two numbers is smaller than 10^(-numberOfDecimalPlaces). We divide by - /// two so that we have half the range on each side of the numbers, e.g. if decimalPlaces == 2, then 0.01 will equal between + /// two so that we have half the range on each side of the numbers, e.g. if == 2, then 0.01 will equal between /// 0.005 and 0.015, but not 0.02 and not 0.00 /// /// @@ -939,13 +1071,72 @@ namespace MathNet.Numerics /// if both doubles are equal to each other within the specified number of decimal places; otherwise . private static bool AlmostEqualWithAbsoluteDecimalPlaces(this double a, double b, int decimalPlaces) { - double decimalPlaceMagnitude = System.Math.Pow(10, -(decimalPlaces - 1)); + double decimalPlaceMagnitude = Math.Pow(10, -(decimalPlaces - 1)); // The values are equal if the difference between the two numbers is smaller than // 10^(-numberOfDecimalPlaces). We divide by two so that we have half the range // on each side of the numbers, e.g. if decimalPlaces == 2, // then 0.01 will equal between 0.005 and 0.015, but not 0.02 and not 0.00 - return System.Math.Abs((a - b)) < decimalPlaceMagnitude / 2.0; + return Math.Abs((a - b)) < decimalPlaceMagnitude / 2.0; + } + + /// + /// Compares two doubles and determines if they are equal to within the tolerance or not. Equality comparison is based on the binary representation. + /// + /// + /// + /// Determines the 'number' of floating point numbers between two values (i.e. the number of discrete steps + /// between the two numbers) and then checks if that is within the specified tolerance. So if a tolerance + /// of 1 is passed then the result will be true only if the two numbers have the same binary representation + /// OR if they are two adjacent numbers that only differ by one step. + /// + /// + /// The comparison method used is explained in http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm . The article + /// at http://www.extremeoptimization.com/resources/Articles/FPDotNetConceptsAndFormats.aspx explains how to transform the C code to + /// .NET enabled code without using pointers and unsafe code. + /// + /// + /// The first value. + /// The second value. + /// The maximum number of floating point values between the two values. Must be 1 or larger. + /// if both doubles are equal to each other within the specified tolerance; otherwise . + /// + /// Thrown if is smaller than one. + /// + public static bool AlmostEqual(this double a, double b, long maxNumbersBetween) + { + // Make sure maxNumbersBetween is non-negative and small enough that the + // default NAN won't compare as equal to anything. + if (maxNumbersBetween < 1) + { + throw new ArgumentOutOfRangeException("maxNumbersBetween"); + } + + // If A or B are infinity (positive or negative) then + // only return true if they are exactly equal to each other - + // that is, if they are both infinities of the same sign. + if (double.IsInfinity(a) || double.IsInfinity(b)) + { + return a == b; + } + + // If A or B are a NAN, return false. NANs are equal to nothing, + // not even themselves. + if (double.IsNaN(a) || double.IsNaN(b)) + { + return false; + } + + // Get the first double and convert it to an integer value (by using the binary representation) + long firstUlong = GetDirectionalLongFromDouble(a); + + // Get the second double and convert it to an integer value (by using the binary representation) + long secondUlong = GetDirectionalLongFromDouble(b); + + // Now compare the values. + // Note that this comparison can overflow so we'll approach this differently + // Do note that we could overflow this way too. We should probably check that we don't. + return (a > b) ? (secondUlong + maxNumbersBetween >= firstUlong) : (firstUlong + maxNumbersBetween >= secondUlong); } /// @@ -976,7 +1167,7 @@ namespace MathNet.Numerics /// /// /// The values are equal if the difference between the two numbers is smaller than 10^(-numberOfDecimalPlaces). We divide by - /// two so that we have half the range on each side of the numbers, e.g. if decimalPlaces == 2, then 0.01 will equal between + /// two so that we have half the range on each side of the numbers, e.g. if == 2, then 0.01 will equal between /// 0.005 and 0.015, but not 0.02 and not 0.00 /// /// @@ -1025,7 +1216,7 @@ namespace MathNet.Numerics /// /// /// The values are equal if the difference between the two numbers is smaller than 10^(-numberOfDecimalPlaces). We divide by - /// two so that we have half the range on each side of the numbers, e.g. if decimalPlaces == 2, then 0.01 will equal between + /// two so that we have half the range on each side of thg. if == 2, then 0.01 will equal between /// 0.005 and 0.015, but not 0.02 and not 0.00 /// /// @@ -1051,7 +1242,7 @@ namespace MathNet.Numerics /// /// The first value. /// The second value. - /// The maximum error in terms of Units in Last Place (ulps), i.e. the maximum number of decimals that may be different. Must be 1 or larger. + /// The maximum error in terms of Units in Last Place (ulps), i.e. the maximum number of decimals that may be different. Must be 1 or larger. /// /// /// diff --git a/src/Native/Native.csproj b/src/Native/Native.csproj index dffba52b..6064b589 100644 --- a/src/Native/Native.csproj +++ b/src/Native/Native.csproj @@ -119,6 +119,9 @@ Interpolation\SplineBoundaryCondition.cs + + IPrecisionSupport.cs + NumberTheory\IntegerTheory.cs