From 8736508a2fd50bd5690b5416e1d46b660e17bfaf Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Mon, 28 Dec 2015 10:19:27 +0100 Subject: [PATCH] Precision: single-precision EpsilonOf, simplify PCL variation --- .../Differentiation/NumericalDerivative.cs | 2 +- src/Numerics/Precision.cs | 169 ++++++++++-------- 2 files changed, 100 insertions(+), 71 deletions(-) diff --git a/src/Numerics/Differentiation/NumericalDerivative.cs b/src/Numerics/Differentiation/NumericalDerivative.cs index 58534dea..7f042f7f 100644 --- a/src/Numerics/Differentiation/NumericalDerivative.cs +++ b/src/Numerics/Differentiation/NumericalDerivative.cs @@ -333,7 +333,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 handle. diff --git a/src/Numerics/Precision.cs b/src/Numerics/Precision.cs index dfe75eb2..9d4db21c 100644 --- a/src/Numerics/Precision.cs +++ b/src/Numerics/Precision.cs @@ -29,10 +29,8 @@ // using System; - -#if PORTABLE +using System.Runtime; using System.Runtime.InteropServices; -#endif namespace MathNet.Numerics { @@ -115,14 +113,14 @@ namespace MathNet.Numerics 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. + /// Actual double precision 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. + /// Actual 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. /// On a standard machine this is equivalent to `PositiveDoublePrecision`. /// @@ -164,12 +162,7 @@ namespace MathNet.Numerics // Note that we need the absolute value of the input because Log10 doesn't // work for negative numbers (obviously). double magnitude = Math.Log10(Math.Abs(value)); - -#if PORTABLE var truncated = (int)Truncate(magnitude); -#else - var truncated = (int) Math.Truncate(magnitude); -#endif // 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 @@ -196,12 +189,7 @@ namespace MathNet.Numerics // Note that we need the absolute value of the input because Log10 doesn't // work for negative numbers (obviously). var magnitude = Convert.ToSingle(Math.Log10(Math.Abs(value))); - -#if PORTABLE var truncated = (int)Truncate(magnitude); -#else - var truncated = (int) Math.Truncate(magnitude); -#endif // 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 @@ -227,22 +215,6 @@ namespace MathNet.Numerics return value*Math.Pow(10, -magnitude); } - /// - /// Gets the equivalent long value for the given double value. - /// - /// The double value which should be turned into a long value. - /// - /// The resulting long value. - /// - static long AsInt64(double value) - { -#if PORTABLE - return DoubleToInt64Bits(value); -#else - return BitConverter.DoubleToInt64Bits(value); -#endif - } - /// /// Returns a 'directional' long value. This is a long value which acts the same as a double, /// e.g. a negative double value will return a negative double value starting at 0 and going @@ -253,7 +225,7 @@ namespace MathNet.Numerics static long AsDirectionalInt64(double value) { // Convert in the normal way. - long result = AsInt64(value); + long result = DoubleToInt64Bits(value); // Now find out where we're at in the range // If the value is larger/equal to zero then we can just return the value @@ -271,7 +243,7 @@ namespace MathNet.Numerics static int AsDirectionalInt32(float value) { // Convert in the normal way. - int result = FloatToInt32Bits(value); + int result = SingleToInt32Bits(value); // Now find out where we're at in the range // If the value is larger/equal to zero then we can just return the value @@ -304,10 +276,10 @@ namespace MathNet.Numerics // Translate the bit pattern of the double to an integer. // Note that this leads to: // double > 0 --> long > 0, growing as the double value grows - // double < 0 --> long < 0, increasing in absolute magnitude as the double + // double < 0 --> long < 0, increasing in absolute magnitude as the double // gets closer to zero! // i.e. 0 - double.epsilon will give the largest long value! - long intValue = AsInt64(value); + long intValue = DoubleToInt64Bits(value); if (intValue < 0) { intValue -= count; @@ -323,13 +295,9 @@ namespace MathNet.Numerics return 0; } - // Note that not all long values can be translated into double values. There's a whole bunch of them + // Note that not all long values can be translated into double values. There's a whole bunch of them // which return weird values like infinity and NaN -#if PORTABLE return Int64BitsToDouble(intValue); -#else - return BitConverter.Int64BitsToDouble(intValue); -#endif } /// @@ -357,12 +325,12 @@ namespace MathNet.Numerics // Translate the bit pattern of the double to an integer. // Note that this leads to: // double > 0 --> long > 0, growing as the double value grows - // double < 0 --> long < 0, increasing in absolute magnitude as the double + // double < 0 --> long < 0, increasing in absolute magnitude as the double // gets closer to zero! // i.e. 0 - double.epsilon will give the largest long value! - long intValue = AsInt64(value); + long intValue = DoubleToInt64Bits(value); - // If the value is zero then we'd really like the value to be -0. So we'll make it -0 + // If the value is zero then we'd really like the value to be -0. So we'll make it -0 // and then everything else should work out. if (intValue == 0) { @@ -379,13 +347,9 @@ namespace MathNet.Numerics intValue -= count; } - // Note that not all long values can be translated into double values. There's a whole bunch of them + // Note that not all long values can be translated into double values. There's a whole bunch of them // which return weird values like infinity and NaN -#if PORTABLE return Int64BitsToDouble(intValue); -#else - return BitConverter.Int64BitsToDouble(intValue); -#endif } /// @@ -506,10 +470,10 @@ namespace MathNet.Numerics // Translate the bit pattern of the double to an integer. // Note that this leads to: // double > 0 --> long > 0, growing as the double value grows - // double < 0 --> long < 0, increasing in absolute magnitude as the double + // double < 0 --> long < 0, increasing in absolute magnitude as the double // gets closer to zero! // i.e. 0 - double.epsilon will give the largest long value! - long intValue = AsInt64(value); + long intValue = DoubleToInt64Bits(value); #if PORTABLE // We need to protect against over- and under-flow of the intValue when @@ -542,7 +506,7 @@ namespace MathNet.Numerics // IntValue is positive var topRangeEnd = long.MaxValue - intValue < maxNumbersBetween // Overflow, 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 + // Also we couldn't translate it back to a double, so we'll return Double.MaxValue ? double.MaxValue // No troubles here : Int64BitsToDouble(intValue + maxNumbersBetween); @@ -652,7 +616,7 @@ namespace MathNet.Numerics /// public static Tuple RangeOfMatchingNumbers(this double value, double relativeDifference) { - // Make sure the relative is non-negative + // Make sure the relative is non-negative if (relativeDifference < 0) { throw new ArgumentOutOfRangeException("relativeDifference"); @@ -675,7 +639,7 @@ namespace MathNet.Numerics // so return the ulps counts for the difference. if (value.Equals(0)) { - var v = AsInt64(relativeDifference); + var v = DoubleToInt64Bits(relativeDifference); return new Tuple(v, v); } @@ -749,7 +713,6 @@ namespace MathNet.Numerics return double.NaN; } -#if PORTABLE long signed64 = DoubleToInt64Bits(value); if (signed64 == 0) { @@ -761,19 +724,35 @@ namespace MathNet.Numerics return Int64BitsToDouble(signed64) - value; } return value - Int64BitsToDouble(signed64); -#else - long signed64 = BitConverter.DoubleToInt64Bits(value); - if (signed64 == 0) + } + + /// + /// Evaluates the minimum distance to the next distinguishable number near the argument value. + /// + /// The value used to determine the minimum distance. + /// + /// Relative Epsilon (positive float or NaN). + /// + /// Evaluates the negative epsilon. The more common positive epsilon is equal to two times this negative epsilon. + /// + public static float EpsilonOf(this float value) + { + if (float.IsInfinity(value) || float.IsNaN(value)) { - signed64++; - return BitConverter.Int64BitsToDouble(signed64) - value; + return float.NaN; } - if (signed64-- < 0) + + int signed32 = SingleToInt32Bits(value); + if (signed32 == 0) { - return BitConverter.Int64BitsToDouble(signed64) - value; + signed32++; + return Int32BitsToSingle(signed32) - value; } - return value - BitConverter.Int64BitsToDouble(signed64); -#endif + if (signed32-- < 0) + { + return Int32BitsToSingle(signed32) - value; + } + return value - Int32BitsToSingle(signed32); } /// @@ -781,13 +760,25 @@ namespace MathNet.Numerics /// /// The value used to determine the minimum distance. /// Relative Epsilon (positive double or NaN) - /// Evaluates the positive epsilon. See also + /// Evaluates the positive epsilon. See also /// public static double PositiveEpsilonOf(this double value) { return 2*EpsilonOf(value); } + /// + /// Evaluates the minimum distance to the next distinguishable number near the argument value. + /// + /// The value used to determine the minimum distance. + /// Relative Epsilon (positive float or NaN) + /// Evaluates the positive epsilon. See also + /// + public static float PositiveEpsilonOf(this float value) + { + return 2 * EpsilonOf(value); + } + /// /// Converts a float value to a bit array stored in an int. /// @@ -795,11 +786,12 @@ namespace MathNet.Numerics /// The bit array. static int FloatToInt32Bits(float value) { - return BitConverter.ToInt32(BitConverter.GetBytes(value), 0); + return SingleToInt32Bits(value); } + /// - /// Calculates the actual positive double precision machine epsilon - the smallest number that can be added to 1, yielding a results different than 1. + /// Calculates the actual (negative) double precision 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. /// /// Positive Machine epsilon @@ -828,24 +820,51 @@ namespace MathNet.Numerics return eps; } + [TargetedPatchingOptOut("Performance critical to inline this type of method across NGen image boundaries")] + static double Truncate(double value) + { #if PORTABLE + return value >= 0.0 ? Math.Floor(value) : Math.Ceiling(value); +#else + return Math.Truncate(value); +#endif + } + + [TargetedPatchingOptOut("Performance critical to inline this type of method across NGen image boundaries")] static long DoubleToInt64Bits(double value) { - var union = new DoubleLongUnion {Double = value}; +#if PORTABLE + var union = new DoubleLongUnion { Double = value }; return union.Int64; +#else + return BitConverter.DoubleToInt64Bits(value); +#endif } + [TargetedPatchingOptOut("Performance critical to inline this type of method across NGen image boundaries")] static double Int64BitsToDouble(long value) { - var union = new DoubleLongUnion {Int64 = value}; +#if PORTABLE + var union = new DoubleLongUnion { Int64 = value }; return union.Double; +#else + return BitConverter.Int64BitsToDouble(value); +#endif } - static double Truncate(double value) + static int SingleToInt32Bits(float value) { - return value >= 0.0 ? Math.Floor(value) : Math.Ceiling(value); + var union = new SingleIntUnion { Single = value }; + return union.Int32; } + static float Int32BitsToSingle(int value) + { + var union = new SingleIntUnion { Int32 = value }; + return union.Single; + } + +#if PORTABLE [StructLayout(LayoutKind.Explicit)] struct DoubleLongUnion { @@ -856,5 +875,15 @@ namespace MathNet.Numerics public long Int64; } #endif + + [StructLayout(LayoutKind.Explicit)] + struct SingleIntUnion + { + [FieldOffset(0)] + public float Single; + + [FieldOffset(0)] + public int Int32; + } } }