Skip to content

Commit

Permalink
fsum just compress option
Browse files Browse the repository at this point in the history
  • Loading branch information
AnthonyLloyd committed Nov 23, 2023
1 parent 2bb539b commit 2e43ca5
Show file tree
Hide file tree
Showing 2 changed files with 1 addition and 133 deletions.
83 changes: 1 addition & 82 deletions Tests/MathX.cs
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ public static double NSum(this IEnumerable<double> values)
}

/// <summary>Shewchuk summation</summary>
public static double FSum(this double[] values, bool compress = false, bool compresstwosum = false, bool renormalise = false)
public static double FSum(this double[] values, bool compress = false)
{
if (values.Length < 3) return values.Length == 2 ? values[0] + values[1] : values.Length == 1 ? values[0] : 0.0;
Span<double> partials = stackalloc double[16];
Expand Down Expand Up @@ -114,14 +114,6 @@ public static double FSum(this double[] values, bool compress = false, bool comp
partials = partials[..count];
if (compress)
Compress(ref lo, ref partials, ref hi);
if(compresstwosum)
CompressTwoSum(ref lo, ref partials, ref hi);
if (renormalise)
{
partials = [lo, .. partials, hi];
lo = 0; hi = 0;
Renormalise(ref partials);
}
foreach (var p in partials)
lo += p;
}
Expand Down Expand Up @@ -161,79 +153,6 @@ static void Compress(ref double lo, ref Span<double> partials, ref double hi)
partials = partials[..top];
}

static void CompressTwoSum(ref double lo, ref Span<double> partials, ref double hi)
{
double q;
hi = TwoSum(hi, partials[^1], out var Q);
var bottom = partials.Length;
for (int i = partials.Length - 2; i >= 0; i--)
{
Q = TwoSum(Q, partials[i], out q);
if (q != 0.0)
{
partials[--bottom] = Q;
Q = q;
}
}
lo = TwoSum(Q, lo, out q);
if (q != 0.0)
{
partials[--bottom] = lo;
lo = q;
}
if (bottom == partials.Length) { partials = []; return; }
Q = TwoSum(partials[bottom], lo, out lo);
var top = 0;
for (int i = bottom + 1; i < partials.Length; i++)
{
Q = TwoSum(partials[i], Q, out q);
if (q != 0.0) partials[top++] = q;
}
hi = TwoSum(hi, Q, out q);
if (q != 0.0) partials[top++] = q;
partials = partials[..top];
}

static void Renormalise(ref Span<double> e)
{
var Q = e[^1];
var bottom = e.Length - 1;
for (int i = e.Length - 2; i >= 0; i--)
{
Q = FastTwoSum(Q, e[i], out var q);
if (q != 0.0)
{
e[bottom--] = Q;
Q = q;
}
}
e[bottom] = Q;
var top = 0;
e[0] = Q;
for (int i = bottom + 1; i < e.Length; i++)
{
Q = TwoSum(e[i], e[top], out var q);
e[top] = Q;
if (q != 0.0)
{
var l = top++ - 1;
while (l >= 0)
{
var c2 = TwoSum(e[l + 1], e[l], out var d2);
if (d2 == 0.0)
{
e[l--] = c2;
top--;
}
else
break;
}
e[top] = q;
}
}
e = e[..(top + 1)];
}

public static double SSum(this double[] values)
{
if (values.Length == 0)
Expand Down
51 changes: 0 additions & 51 deletions Tests/MathXTests.cs
Original file line number Diff line number Diff line change
Expand Up @@ -182,57 +182,6 @@ public void FSum_Shuffle_Error_Distribution_Compress()
}, output.WriteLine/*, time: 10*/);
}

[Fact]
public void FSum_Shuffle_Error_Distribution_Renormalise()
{
genDouble.Array[3, 100]
.SelectMany(a => Gen.Shuffle(a).Select(s => (a, s)))
.Sample((original, shuffled) =>
{
var originalSum = MathX.FSum(original, renormalise: true);
var shuffledSum = MathX.FSum(shuffled, renormalise: true);
return Check.UlpsBetween(originalSum, shuffledSum).ToString().PadLeft(5);
}, output.WriteLine/*, time: 10*/);
}

[Fact]
public void FSum_Shuffle_Error_Distribution_Compress_Renormalise()
{
genDouble.Array[3, 100]
.SelectMany(a => Gen.Shuffle(a).Select(s => (a, s)))
.Sample((original, shuffled) =>
{
var originalSum = MathX.FSum(original, true, true);
var shuffledSum = MathX.FSum(shuffled, true, true);
return Check.UlpsBetween(originalSum, shuffledSum).ToString().PadLeft(5);
}, output.WriteLine/*, time: 10*/);
}

[Fact]
public void FSum_Shuffle_Error_Distribution_Comparison()
{
genDouble.Array[3, 100]
.SelectMany(a => Gen.Shuffle(a).Select(s => (a, s)))
.Sample((original, shuffled) =>
{
var v1 = Check.UlpsBetween(MathX.FSum(original, false, false, false), MathX.FSum(shuffled, false, false, false));
var v2 = Check.UlpsBetween(MathX.FSum(original, true, false, false), MathX.FSum(shuffled, true, false, false));
var v3 = Check.UlpsBetween(MathX.FSum(original, false, true, false), MathX.FSum(shuffled, false, true, false));
var v4 = Check.UlpsBetween(MathX.FSum(original, false, false, true), MathX.FSum(shuffled, false, false, true));
var v5 = Check.UlpsBetween(MathX.FSum(original, true, false, true), MathX.FSum(shuffled, true, false, true));
var v6 = Check.UlpsBetween(MathX.FSum(original, false, true, true), MathX.FSum(shuffled, false, true, true));
return $"{v1}_{v2}_{v3}_{v4}_{v5}_{v6}";
}, output.WriteLine/*, time: 10*/);
}
// Passed Tests.MathXTests.FSum_Shuffle_Error_Distribution_Comparison [14 m 29 s]
// Standard Output Messages:
//| | Count | % | Median | Lower Q | Upper Q | Minimum | Maximum |
//|-------------|------------:|--------:|----------:|----------:|----------:|----------:|--------------:|
//| 0_0_0_0_0_0 | 399,999,616 | 100.00% | 30.0864μs | 12.9196μs | 57.2142μs | 0.2000μs | 16,827.5000μs |
//| 1_0_0_0_0_0 | 299 | 0.00% | 29.4548μs | 16.3683μs | 46.2734μs | 2.5000μs | 91.6000μs |
//| 1_1_1_1_1_1 | 85 | 0.00% | 40.5411μs | 22.8298μs | 60.4956μs | 2.1000μs | 102.9000μs |


[Fact]
public void FSum_Compress_Needed_Example()
{
Expand Down

0 comments on commit 2e43ca5

Please sign in to comment.