168 lines
8.6 KiB
C#
168 lines
8.6 KiB
C#
using Just.PreciseMath.PerformanceComparison.Reference;
|
|
|
|
namespace Just.PreciseMath.PerformanceComparison;
|
|
|
|
/// <summary>
|
|
/// Deterministic operand sets for the implementation comparison: <see cref="Count"/>
|
|
/// values per operation, each with a nonzero high and low component, built from exact
|
|
/// dyadic multiples and held inside a moderate exponent range so that no benchmark
|
|
/// path takes an overflow, underflow, or boundary branch.
|
|
/// </summary>
|
|
/// <remarks>
|
|
/// Both double-double implementations receive bit-identical operands, and the plain
|
|
/// doubles equal the high components, so the binary64 reference measures the leading
|
|
/// component of the same values. The mixed-sign set feeds +, -, *, / and Abs; the
|
|
/// nonnegative set feeds Square and Sqrt.
|
|
/// </remarks>
|
|
internal sealed class PerformanceInputs
|
|
{
|
|
/// <summary>Operands per benchmark invocation. One double-double operation takes only
|
|
/// a few nanoseconds, so an invocation covers many values and the measurement is
|
|
/// normalized back to a single operation.</summary>
|
|
internal const int Count = 200;
|
|
|
|
/// <summary>Mixed-sign operands for +, -, *, / and Abs.</summary>
|
|
internal DoubleDouble[] Left { get; } = new DoubleDouble[Count];
|
|
|
|
internal DoubleDouble[] Right { get; } = new DoubleDouble[Count];
|
|
|
|
/// <summary>Nonnegative operands for Square and Sqrt, with both components positive.</summary>
|
|
internal DoubleDouble[] Positive { get; } = new DoubleDouble[Count];
|
|
|
|
/// <summary>The same operands for the blog reference implementation.</summary>
|
|
internal DoubleDoubleSample[] SampleLeft { get; } = new DoubleDoubleSample[Count];
|
|
|
|
internal DoubleDoubleSample[] SampleRight { get; } = new DoubleDoubleSample[Count];
|
|
|
|
internal DoubleDoubleSample[] SamplePositive { get; } = new DoubleDoubleSample[Count];
|
|
|
|
/// <summary>The high components as plain doubles, for the binary64 reference.</summary>
|
|
internal double[] DoubleLeft { get; } = new double[Count];
|
|
|
|
internal double[] DoubleRight { get; } = new double[Count];
|
|
|
|
internal double[] DoublePositive { get; } = new double[Count];
|
|
|
|
internal static PerformanceInputs Create()
|
|
{
|
|
PerformanceInputs inputs = new PerformanceInputs();
|
|
for (int i = 0; i < Count; i++)
|
|
{
|
|
// Highs are small dyadic multiples of a power of two: exact in binary64,
|
|
// with the significand needing at most nine bits. Lows are multiples of
|
|
// 2^(exponent-57) or 2^(exponent-60), at least two bits below half an ulp
|
|
// of their high, so every pair is normalized and both components survive
|
|
// every operation. The exponents keep highs within 2^+-48, which leaves
|
|
// products, quotients, and squares far from the binary64 exponent limits.
|
|
int leftExponent = (((i / 4) % 5) - 2) * 24;
|
|
int rightExponent = (((i / 4) % 3) - 1) * 16;
|
|
int positiveExponent = ((i % 5) - 2) * 24;
|
|
double leftSign = (i % 2 == 0) ? 1.0 : -1.0;
|
|
double rightSign = (i % 3 == 0) ? 1.0 : -1.0;
|
|
double leftHigh = leftSign * Math.ScaleB(1.25 + ((i % 7) / 32.0), leftExponent);
|
|
double leftLow = leftSign * Math.ScaleB(1.0 + (i % 4), leftExponent - 57);
|
|
double rightHigh = rightSign * Math.ScaleB(1.0 + (((i % 5) - 2) / 128.0), rightExponent);
|
|
// Opposite sign, and a three-bit finer scale than the left low: no
|
|
// addition or multiplication of the two lows can cancel to a value that
|
|
// binary64 would represent without a residual.
|
|
double rightLow = -rightSign * Math.ScaleB(1.0 + ((i / 2) % 3), rightExponent - 60);
|
|
double positiveHigh = Math.ScaleB(1.25 + ((i % 7) / 32.0), positiveExponent);
|
|
double positiveLow = Math.ScaleB(1.0 + (i % 4), positiveExponent - 57);
|
|
|
|
inputs.Left[i] = DoubleDouble.FromComponents(leftHigh, leftLow);
|
|
inputs.Right[i] = DoubleDouble.FromComponents(rightHigh, rightLow);
|
|
inputs.Positive[i] = DoubleDouble.FromComponents(positiveHigh, positiveLow);
|
|
inputs.SampleLeft[i] = DoubleDoubleSample.FromTwoDoubles(leftHigh, leftLow);
|
|
inputs.SampleRight[i] = DoubleDoubleSample.FromTwoDoubles(rightHigh, rightLow);
|
|
inputs.SamplePositive[i] = DoubleDoubleSample.FromTwoDoubles(positiveHigh, positiveLow);
|
|
inputs.DoubleLeft[i] = leftHigh;
|
|
inputs.DoubleRight[i] = rightHigh;
|
|
inputs.DoublePositive[i] = positiveHigh;
|
|
}
|
|
|
|
inputs.Validate();
|
|
return inputs;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Setup-only fixture guard. It checks the operand shape and that every benchmarked
|
|
/// operation consumes both components; it is not an accuracy oracle and does not
|
|
/// compare the two implementations against each other.
|
|
/// </summary>
|
|
private void Validate()
|
|
{
|
|
for (int i = 0; i < Count; i++)
|
|
{
|
|
ValidateOperand(Left[i], nameof(Left), i);
|
|
ValidateOperand(Right[i], nameof(Right), i);
|
|
ValidateOperand(Positive[i], nameof(Positive), i);
|
|
ValidateMatch(Left[i], SampleLeft[i], nameof(Left), i);
|
|
ValidateMatch(Right[i], SampleRight[i], nameof(Right), i);
|
|
ValidateMatch(Positive[i], SamplePositive[i], nameof(Positive), i);
|
|
if (DoubleLeft[i] != Left[i].High || DoubleRight[i] != Right[i].High || DoublePositive[i] != Positive[i].High)
|
|
{
|
|
throw new InvalidOperationException($"[{i}] the binary64 operands do not match the leading components.");
|
|
}
|
|
if (!(Positive[i].High > 0.0) || !(Positive[i].Low > 0.0))
|
|
{
|
|
throw new InvalidOperationException($"[{i}] Square and Sqrt operands must be positive in both components.");
|
|
}
|
|
}
|
|
|
|
for (int i = 0; i < Count; i++)
|
|
{
|
|
CheckResidual("Left + Right", Left[i] + Right[i], i);
|
|
CheckResidual("Left - Right", Left[i] - Right[i], i);
|
|
CheckResidual("Left * Right", Left[i] * Right[i], i);
|
|
CheckResidual("Left / Right", Left[i] / Right[i], i);
|
|
CheckResidual("Positive * Positive", DoubleDouble.Square(Positive[i]), i);
|
|
CheckResidual("Sqrt(Positive)", DoubleDouble.Sqrt(Positive[i]), i);
|
|
CheckResidual("Abs(Left)", DoubleDouble.Abs(Left[i]), i);
|
|
CheckResidual("SampleLeft + SampleRight", SampleLeft[i] + SampleRight[i], i);
|
|
CheckResidual("SampleLeft - SampleRight", SampleLeft[i] - SampleRight[i], i);
|
|
CheckResidual("SampleLeft * SampleRight", SampleLeft[i] * SampleRight[i], i);
|
|
CheckResidual("SampleLeft / SampleRight", SampleLeft[i] / SampleRight[i], i);
|
|
CheckResidual("SamplePositive.Square()", SamplePositive[i].Square(), i);
|
|
CheckResidual("SamplePositive.Sqrt()", SamplePositive[i].Sqrt(), i);
|
|
CheckResidual("SampleLeft.Abs()", SampleLeft[i].Abs(), i);
|
|
}
|
|
}
|
|
|
|
private static void ValidateOperand(DoubleDouble value, string name, int index)
|
|
{
|
|
// FromComponents renormalizes an arbitrary pair, so reproducing the value
|
|
// bit-for-bit is the normalization check for a purely dyadic construction.
|
|
bool normalized = DoubleDouble.FromComponents(value.High, value.Low).Equals(value);
|
|
if (!DoubleDouble.IsFinite(value) || value.High == 0.0 || value.Low == 0.0 || !normalized
|
|
|| Math.Abs(Math.ILogB(value.High)) > 100)
|
|
{
|
|
throw new InvalidOperationException($"{name}[{index}] is not an ordinary normalized value with both components used.");
|
|
}
|
|
}
|
|
|
|
private static void ValidateMatch(DoubleDouble value, DoubleDoubleSample sample, string name, int index)
|
|
{
|
|
bool normalized = DoubleDoubleSample.FromTwoDoubles(sample.Hi, sample.Lo).Equals(sample);
|
|
if (sample.Hi != value.High || sample.Lo != value.Low || !normalized)
|
|
{
|
|
throw new InvalidOperationException($"{name}[{index}] differs between the two double-double implementations.");
|
|
}
|
|
}
|
|
|
|
private static void CheckResidual(string operation, DoubleDouble result, int index)
|
|
{
|
|
if (!DoubleDouble.IsFinite(result) || result.Low == 0.0)
|
|
{
|
|
throw new InvalidOperationException($"{operation}[{index}] left the ordinary range or lost its low component.");
|
|
}
|
|
}
|
|
|
|
private static void CheckResidual(string operation, DoubleDoubleSample result, int index)
|
|
{
|
|
if (!double.IsFinite(result.Hi) || result.Lo == 0.0)
|
|
{
|
|
throw new InvalidOperationException($"{operation}[{index}] left the ordinary range or lost its low component.");
|
|
}
|
|
}
|
|
}
|