NBody benchmark program

Benchmark.cs

using System;
using System.Collections.Generic;
using System.Text;
using System.Linq;
using System.Runtime.InteropServices;

namespace Benchmark {

  public static class Extensions {
    public static string ToStringReset(this StringBuilder sb) {
      var result = sb.ToString();
      sb.Clear();
      return result;
    }
    public static double Sorted_QMedian(this List<double> list, int quarter) {
      int len = list.Count;
      if (len == 0) { return 0; }

      int idx = (len / 4) * quarter;
      double med = list[idx];
      return len % 2 != 0 ? med : ((list[idx - 1] + med) / 2);
    }
  }

  public class BenchmarkRunner {
    private static StringBuilder SB = new StringBuilder(2048);

    public int Iterations { get; set; } = 3;

    public int RepeatTimes { get; set; } = 2;

    public string[] CaseArgs { get; set; } = new string[] { };

    public Statistics BaseLine { get; set; }

    public List<Statistics> CaseList { get; private set; } = new List<Statistics>(16);

    public Statistics AddCase(string name, Action<string[]> action) { 
      Statistics stats = new Statistics(this, name, action);
      CaseList.Add(stats);
      if (this.BaseLine == null) { this.BaseLine = stats; }
      return stats;
    }

    public void Reset() {
      foreach(var stat in CaseList) { stat.Reset(); }
    }


    public void Run() {
      for (int iter = this.Iterations; iter > 0; --iter) {
        for (int i = 0; i < CaseList.Count; ++i) {
          GC.Collect(10, GCCollectionMode.Forced, true, true);
          Statistics stats = CaseList[i];
          Console.Out.WriteLine("-------------------------------\n Beginning: " + stats.Name);
          for (int j = 0; j < RepeatTimes; ++j) {
            var timer = Timer.Start();
            stats.Action(CaseArgs);
            var seconds = timer.Seconds;
            stats.AddTime(seconds);
            Console.Out.WriteLine(" Elapsed: " + seconds.ToString());
          }
        }
      }
      foreach (var stat in CaseList) { stat.Update(); }
    }

    public override string ToString() {
      foreach (var stat in CaseList) { stat.Update(); }

      var columns = StatColumn.AllColumns;
      foreach (var col in columns) { col.UpdateSpacing(CaseList); }
      foreach (var col in columns) { SB.Append(col.Name.PadLeft(col.Spacing)).Append(" | "); }

      var headers = SB.ToStringReset();
      var separator = SB.Append('-', headers.Length - 1).ToStringReset();
      SB.AppendLine(separator).Append(headers);

      foreach (var stat in CaseList) {
        SB.AppendLine();
        foreach (var col in columns) { SB.Append(col.Format(stat)).Append(" | "); }
      }
      SB.AppendLine().Append(separator);
      return SB.ToStringReset();
    }
  }


  public class StatColumn {
    public string Name { get; }
    public Func<Statistics, string> GetValue { get; }
    public virtual int Spacing { get; set; }


    public static readonly StatColumn CaseName = new StatColumn("Name", x => x.Name);
    public static readonly StatColumn Mean = new DecimalColumn("Mean", x => x.Mean);
    public static readonly StatColumn Min = new DecimalColumn("Min", x => x.Min);
    public static readonly StatColumn Max = new DecimalColumn("Max", x => x.Max);
    public static readonly StatColumn Q1 = new DecimalColumn("Q1", x => x.Q1);
    public static readonly StatColumn Median = new DecimalColumn("Median", x => x.Median);
    public static readonly StatColumn Q3 = new DecimalColumn("Q3", x => x.Q3);
    public static readonly StatColumn First = new DecimalColumn("First", x => x.First);
    public static readonly StatColumn Error = new DecimalColumn("Error", x => x.Error);
    public static readonly StatColumn StdErr = new DecimalColumn("StdErr", x => x.StdErr);
    public static readonly StatColumn StdDev = new DecimalColumn("StdDev", x => x.StdDev);
    public static readonly StatColumn TotalTime = new DecimalColumn("Elapsed", x => x.Elapsed);
    public static readonly StatColumn Scale = new StatColumn("Scale", x => x.Scale.ToString("P"));
    public static readonly StatColumn Runs = new IntegerColumn("Runs", x => x.N);
    public static readonly StatColumn[] AllColumns = new[] {
      CaseName, Mean, Median, Min, Max, Q1, Q3, Scale, Runs, TotalTime, First, Error, StdErr, StdDev
    };


    protected StatColumn(string name, Func<Statistics, string> getter) {
      this.Spacing = name.Length;
      this.Name = name;
      this.GetValue = getter;
    }
    public virtual void UpdateSpacing(IEnumerable<Statistics> stats) {
      int longest = this.Name.Length;
      foreach (var stat in stats) { longest = Math.Max(GetValue(stat).Length, longest); }
      this.Spacing = longest;
    }
    public virtual string Format(Statistics stat) {
      return GetValue(stat).PadLeft(this.Spacing);
    }


    public class DecimalColumn : StatColumn {
      public new Func<Statistics, double> GetValue { get; }
      public int LeftSpace { get; set; }
      public int RightSpace { get; set; }
      public override int Spacing => Math.Max(base.Spacing, 1 + LeftSpace + RightSpace);

      public DecimalColumn(string name, Func<Statistics, double> getter) : base(name, null) {
        GetValue = getter;
      }
      public override void UpdateSpacing(IEnumerable<Statistics> stats) {
        int ideal = Math.Max(5, base.Spacing);
        int whole = 1;
        int precision = 0;

        foreach (var stat in stats) {
          double value = this.GetValue(stat);
          string[] parts = value.ToString().Split(".");
          whole = Math.Max(whole, parts[0].Length);
          precision = Math.Max(precision, parts.Length > 1 ? parts[1].Length : 0);
        }
        this.LeftSpace = whole;
        this.RightSpace = Math.Max(0, ideal - whole - 1);
      }

      public override string Format(Statistics stat) {
        double value = GetValue(stat);
        string[] parts = value.ToString().Split(".");

        string whole = parts[0].PadLeft(this.LeftSpace);
        string precision = parts.Length > 1 ? parts[1] : "";
        if (precision.Length > this.RightSpace) {
          precision = precision.Substring(0, this.RightSpace);
        } else if (precision.Length < this.RightSpace) {
          precision = precision.PadRight(this.RightSpace);
        }
        return whole + (precision.Length > 0 ? "." : "") + precision;
      }
    }

    public class IntegerColumn : StatColumn {
      public new Func<Statistics, int> GetValue { get; }
      public IntegerColumn(string name, Func<Statistics, int> getter) : base(name, null) {
        GetValue = getter;
      }
      public override void UpdateSpacing(IEnumerable<Statistics> stats) {
        int largest = 0;
        foreach (var stat in stats) { largest = Math.Max(GetValue(stat), largest); }
        this.Spacing = Math.Max(this.Name.Length, largest.ToString().Length);
      }
      public override string Format(Statistics stat) {
        return GetValue(stat).ToString().PadLeft(this.Spacing);
      }
    }
  }


  public class Statistics {
    private double _elapsed;
    private double _variance;
    private double _mean;
    private double _stdDev;
    private double _stdErr;
    private double _skewness;
    private double _kurtosis;
    private double _error;
    private BenchmarkRunner _owner;
    private List<double> _timings = new List<double>(32);
    private List<double> _sorted = new List<double>(32);
    private List<double> _lowerOutliers = new List<double>();
    private List<double> _upperOutliers = new List<double>();

    public string Name { get; }
    public Action<string[]> Action { get; }

    public List<double> SortedValues => EnsureLatest()._sorted;
    public List<double> LowerOutliers => EnsureLatest()._lowerOutliers;
    public List<double> UpperOutliers => EnsureLatest()._upperOutliers;

    public int N => _timings.Count;
    public double First => _timings.Count > 0 ? _timings[0] : 0;

    public double Elapsed => EnsureLatest()._elapsed;
    public double Mean => EnsureLatest()._mean;

    public double Min => SortedValues.FirstOrDefault();
    public double Max => SortedValues.LastOrDefault();
    public double Median => SortedValues.Sorted_QMedian(2);
    public double Q1 => SortedValues.Sorted_QMedian(1);
    public double Q3 => SortedValues.Sorted_QMedian(3);
    public double InterquartileRange => Q3 - Q1;
    public double LowerFence => Q1 - 1.5 * InterquartileRange;
    public double UpperFence => Q3 + 1.5 * InterquartileRange;

    public double Variance => EnsureLatest()._variance;
    public double StdDev => EnsureLatest()._stdDev;
    public double StdErr => EnsureLatest()._stdErr;
    public double Skewness => EnsureLatest()._skewness;
    public double Kurtosis => EnsureLatest()._kurtosis;
    public double Error => EnsureLatest()._error;

    public double Scale => this.Elapsed / _owner.BaseLine.Elapsed;

    private Statistics EnsureLatest() {
      if (this._timings.Count != this._sorted.Count) { Update(); }
      return this;
    }

    public Statistics(BenchmarkRunner owner, string method, Action<string[]> action) {
      this._owner = owner;
      Name = method;
      Action = action;
    }

    public void Update() {
      _lowerOutliers.Clear();
      _upperOutliers.Clear();
      _sorted.Clear();
      _sorted.AddRange(_timings);
      _sorted.Sort();

      int len = _sorted.Count;
      _elapsed = _sorted.Sum();
      _mean = _elapsed / Math.Max(len, 1);

      double lowerFence = this.LowerFence;
      double upperFence = this.UpperFence;

      _variance = 0;
      _skewness = 0;
      _kurtosis = 0;
      for (int i = 0; i < len; ++i) {
        double val = _sorted[i];
        double offset = val - _mean;
        _variance += Math.Pow(offset, 2);
        _skewness += Math.Pow(offset, 3);
        _kurtosis += Math.Pow(offset, 4);

        if (val < lowerFence) {
          _lowerOutliers.Add(val);
        } else if (val > upperFence) {
          _upperOutliers.Add(val);
        }
      }
      _variance = len <= 1 ? 0 : _variance / (len - 1);
      _stdDev = Math.Sqrt(_variance);
      _stdErr = _stdDev / (len > 0 ? Math.Sqrt(len) : 1);

      _skewness = (len == 0 ? 0 : (_skewness / len)) / Math.Pow(_stdDev, 3);
      _kurtosis = (len == 0 ? 0 : (_kurtosis / len)) / Math.Pow(_stdDev, 4);
      _error = len <= 2 ? double.NaN : StdErr * Helpers.InverseStudent(.1, N - 1);
    }

    public void Reset() {
      this._timings.Clear();
      EnsureLatest();
    }

    public Statistics AddTime(double time) {
      if (!double.IsNaN(time)) {
        this._timings.Add(time);
      }
      return this;
    }

  }


  public struct Timer {
    [DllImport("kernel32.dll")] private static extern bool QueryPerformanceCounter(out long ticks);
    [DllImport("kernel32.dll")] private static extern bool QueryPerformanceFrequency(out long frequency);

    private static long GetTicks() { QueryPerformanceCounter(out long ticks); return ticks; }

    private static readonly long FREQUENCY;
    static Timer() {
      QueryPerformanceFrequency(out FREQUENCY);
    }


    public long StartTime;

    public double Seconds => (double)(GetTicks() - StartTime) / FREQUENCY;

    public double MicroSeconds => Seconds * 1e6;

    public static Timer Start() {
      return new Timer { StartTime = GetTicks() };
    }
  }


  static class Helpers {
    public static double InverseStudent(double p, double n) {
      double lower = 0.0;
      double upper = 1000.0;
      while (upper - lower > 1e-9) {
        double t = (lower + upper) / 2;
        double p2 = StudentTwoTail(t, n);
        if (p2 < p) { upper = t; } else { lower = t; }
      }
      return (lower + upper) / 2;
    }

    static double StudentTwoTail(double t, double n) {
      t = t * t;
      double y = t / n;
      double b = y + 1.0;
      double z = 1;
      double a;
      int nn = (int)Math.Round(n);
      if (Math.Abs(n - nn) > 1e-9 || n >= 20 || t < n && n > 200) {
        if (y > 1.0e-6) { y = Math.Log(b); }
        a = n - 0.5;
        b = 48.0 * (a * a);
        y = a * y;
        y = (((((-0.4 * y - 3.3) * y - 24.0) * y - 85.5) / (0.8 * (y * y) + 100.0 + b) + y + 3.0) / b + 1.0) * Math.Sqrt(y);
        return 2 * Gauss(-y);
      } else if (n < 20 && t < 4.0) {
        y = Math.Sqrt(y);
        a = y;
        if (nn == 1) { a = 0; }
      } else {
        a = Math.Sqrt(b);
        y = a * nn;
        int j = 0;
        while (Math.Abs(a - z) > 0) {
          j += 2;
          z = a;
          y *= (j - 1) / (b * j);
          a += y / (nn + j);
        }
        nn += 2;
        z = 0;
        y = 0;
        a = -a;
      }
      while (true) {
        nn -= 2;
        if (nn > 1) { a = (nn - 1) / (b * nn) * a + y; } else break;
      }
      a = nn == 0 ? a / Math.Sqrt(b) : (Math.Atan(y) + a / b) * 2 / Math.PI;
      return z - a;
    }

    private static double Gauss(double x) {
      double z;
      if (Math.Abs(x) < 1e-9) { z = 0.0; } else {
        double y = Math.Abs(x) / 2;
        if (y >= 3.0) { z = 1.0; } else if (y < 1.0) {
          double w = y * y;
          z = ((((((((0.000124818987 * w - 0.001075204047) * w + 0.005198775019) * w - 0.019198292004) * w + 0.059054035642) *
            w - 0.151968751364) * w + 0.319152932694) * w - 0.531923007300) * w + 0.797884560593) * y * 2.0;
        } else {
          y = y - 2.0;
          z = (((((((((((((-0.000045255659 * y + 0.000152529290) * y - 0.000019538132) * y - 0.000676904986) * y + 0.001390604284) *
            y - 0.000794620820) * y - 0.002034254874) * y + 0.006549791214) * y - 0.010557625006) * y + 0.011630447319) *
            y - 0.009279453341) * y + 0.005353579108) * y - 0.002141268741) * y + 0.000535310849) * y + 0.999936657524;
        }
      }
      return x > 0.0 ? (z + 1.0) / 2 : (1.0 - z) / 2;
    }
  }
}

Program.cs

using System;

using Scalar = NBodyCore.Scalar;
using AVX = NBodyCore.AVX;
using Official = NBodyCore.Official;
using System.Threading;

class Program {

  static void Main(string[] args) {
    Thread.CurrentThread.Priority = ThreadPriority.Highest;

    Benchmark.BenchmarkRunner bench = new Benchmark.BenchmarkRunner();
    bench.Iterations = 2;
    bench.RepeatTimes = 4;

    //	bench.RunCase("Official.Number3", () => Official.Number3.NBody.Main(args));
    //	bench.RunCase("Official.Number8", () => Official.Number7.NBody.Main(args));
    //	bench.BaseLine = bench.RunCase("Official.OriginalNum7", () => Official.OriginalNum7.NBody.Main(args));
    //	bench.RunCase("Derek_Scalar", () => Official.OriginalNum7.NBody.Main(args));
    //	bench.RunCase("Derek.CoreNum7_2", () => Exp.CoreNum7_2.NBody.Main(args));
    //	bench.RunCase("Derek.VectorAVXCorrect2", () => Exp.VectorAVXCorrect2.NBody.Main(args));
    //	//bench.RunCase("Derek.VectorAVXCorrect", () => Derek.VectorAVXCorrect.NBody.Main(args));

    //	//bench.RunCase("Derek.VectorAVX", () => Derek.VectorAVX.NBody.Main(args));
    //	bench.RunCase("Official.CoreNum3", () => Official.CoreNum3.NBody.Main(args));

    //	//bench.RunCase("Derek.CoreNum7", () => Derek.CoreNum7.NBody.Main(args));

    //	//bench.RunCase("Derek.PassObjToMethods", () => Derek.PassObjToMethods.NBody.Main(args));
    //	//bench.RunCase("Derek.PassStructToMethods", () => Derek.PassStructToMethods.NBody.Main(args));
    //	//bench.RunCase("Derek.PassRefStructToMethods", () => Derek.PassRefStructToMethods.NBody.Main(args));


    //bench.AddCase("AVX.MaxTheoretical", AVX.MaxTheoretical.NBody.Main);
    //bench.BaseLine = bench.AddCase("Official.Number8", Official.Number8.NBody.Main);
    //bench.AddCase("Avx-5.1", AVX.V5_1.NBody.Main);
    //bench.AddCase("Official.Number3", Official.Number3.NBody.Main);
    //bench.AddCase("Classic.Number7", Official.ClassicNumber7.NBody.Main);

    bench.AddCase("Official.Number8", Official.Number8.NBody.Main);
    bench.AddCase("Official.AvxNumber1", Official.AvxNumber1.NBody.Main);
    bench.AddCase("Official.AvxNumber2", AVX.AvxNumber1_4.NBody.Main);
    bench.AddCase("AVX-4 Bodies at a Time", AVX.MaxTheoretical2.NBody.Main);
    //bench.AddCase("MaxTheoreticalC", AVX.MaxTheoretical3.NBody.Main);


    //bench.AddCase("AvxNumber1_4", AVX.AvxNumber1_4.NBody.Main);
    //bench.AddCase("AvxNumber1_4ref", AVX.AvxNumber1_4ref.NBody.Main);


    bench.CaseArgs = new[] { "100000" };
    bench.Run();
    Console.Out.WriteLine(bench);
    bench.Reset();

    bench.CaseArgs = new[] { "1000000" };
    bench.Run();
    Console.Out.WriteLine(bench);
    bench.Reset();

    bench.CaseArgs = new[] { "2000000" };
    bench.Run();
    Console.Out.WriteLine(bench);
    bench.Reset();

    //Console.Clear();

    bench.CaseArgs = new[] { "50000000" };
    while (true) {
      bench.Run();
      Console.Out.WriteLine(bench.ToString());
      //Console.In.Read();
    }
  }
}