tsai
/
mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Precision: single-precision EpsilonOf, simplify PCL variation

pull/362/head
Christoph Ruegg 11 years ago
parent
3a72370069
commit
8736508a2f
2 changed files with 100 additions and 71 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 2
      src/Numerics/Differentiation/NumericalDerivative.cs
  2. 169
      src/Numerics/Precision.cs

2
src/Numerics/Differentiation/NumericalDerivative.cs
View File

@ -333,7 +333,7 @@ namespace MathNet.Numerics.Differentiation
/// </summary>
/// <remarks>
/// This function recursively uses <see cref="EvaluatePartialDerivative(Func&lt;double[], double&gt;, double[], int, int, double?)"/> to evaluate mixed partial derivative.
/// Therefore, it is more efficient to call <see cref="EvaluatePartialDerivative(Func&lt;double[], double&gt;, double[], int, int, double?)"/> for higher order derivatives of
/// Therefore, it is more efficient to call <see cref="EvaluatePartialDerivative(Func&lt;double[], double&gt;, double[], int, int, double?)"/> for higher order derivatives of
/// a single independent variable.
/// </remarks>
/// <param name="f">Multivariate function handle.</param>

169
src/Numerics/Precision.cs
View File

@ -29,10 +29,8 @@
// </copyright>
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;
/// <summary>
/// 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`.
/// </summary>
public static readonly double MachineEpsilon = MeasureMachineEpsilon();
/// <summary>
/// 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`.
/// </summary>
@ -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);
}
/// <summary>
/// Gets the equivalent <c>long</c> value for the given <c>double</c> value.
/// </summary>
/// <param name="value">The <c>double</c> value which should be turned into a <c>long</c> value.</param>
/// <returns>
/// The resulting <c>long</c> value.
/// </returns>
static long AsInt64(double value)
{
#if PORTABLE
return DoubleToInt64Bits(value);
#else
return BitConverter.DoubleToInt64Bits(value);
#endif
}
/// <summary>
/// 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
}
/// <summary>
@ -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
}
/// <summary>
@ -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
/// </returns>
public static Tuple<long, long> 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<long, long>(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)
}
/// <summary>
/// Evaluates the minimum distance to the next distinguishable number near the argument value.
/// </summary>
/// <param name="value">The value used to determine the minimum distance.</param>
/// <returns>
/// Relative Epsilon (positive float or NaN).
/// </returns>
/// <remarks>Evaluates the <b>negative</b> epsilon. The more common positive epsilon is equal to two times this negative epsilon.</remarks>
/// <seealso cref="PositiveEpsilonOf(float)"/>
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);
}
/// <summary>
@ -781,13 +760,25 @@ namespace MathNet.Numerics
/// </summary>
/// <param name="value">The value used to determine the minimum distance.</param>
/// <returns>Relative Epsilon (positive double or NaN)</returns>
/// <remarks>Evaluates the <b>positive</b> epsilon. See also <see cref="EpsilonOf"/></remarks>
/// <remarks>Evaluates the <b>positive</b> epsilon. See also <see cref="EpsilonOf(double)"/></remarks>
/// <seealso cref="EpsilonOf(double)"/>
public static double PositiveEpsilonOf(this double value)
{
return 2*EpsilonOf(value);
}
/// <summary>
/// Evaluates the minimum distance to the next distinguishable number near the argument value.
/// </summary>
/// <param name="value">The value used to determine the minimum distance.</param>
/// <returns>Relative Epsilon (positive float or NaN)</returns>
/// <remarks>Evaluates the <b>positive</b> epsilon. See also <see cref="EpsilonOf(float)"/></remarks>
/// <seealso cref="EpsilonOf(float)"/>
public static float PositiveEpsilonOf(this float value)
{
return 2 * EpsilonOf(value);
}
/// <summary>
/// Converts a float value to a bit array stored in an int.
/// </summary>
@ -795,11 +786,12 @@ namespace MathNet.Numerics
/// <returns>The bit array.</returns>
static int FloatToInt32Bits(float value)
{
return BitConverter.ToInt32(BitConverter.GetBytes(value), 0);
return SingleToInt32Bits(value);
}
/// <summary>
/// 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.
/// </summary>
/// <returns>Positive Machine epsilon</returns>
@ -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;
}
}
}

Write Preview
Loading…
Cancel
Save