Benchmarking Different Hypotenuse Algorithms

HypothenuseBenchmarks.cs

/*
    Starting Speed Test for 100000000 iterations
    Distance_Normal:                        17187ms    AvgErr: 0.00%    MaxErr: 0.00%
    Distance_Sqrt2:                         585  ms    AvgErr: 5.28%    MaxErr: 8.24%
    Distance_4142:                          623  ms    AvgErr: 5.28%    MaxErr: 8.24%
    Distance_4142_Cast:                     3824 ms    AvgErr: 5.28%    MaxErr: 8.24%
    Distance_4142_Int_Only:                 584  ms    AvgErr: 5.28%    MaxErr: 8.24%
    Distance_4142_Int_Only_Manual_Inlined:  666  ms    AvgErr: 5.28%    MaxErr: 8.24%
    Distance_AShift1:                       606  ms    AvgErr: 8.84%    MaxErr:11.80%
    Distance_337:                           613  ms    AvgErr: 3.38%    MaxErr: 5.52%
    Distance_918:                           1555 ms    AvgErr: 1.37%    MaxErr: 2.63%
    Distance_428:                           1049 ms    AvgErr: 0.55%    MaxErr: 1.40%
    Distance_FlipCode:                      608  ms    AvgErr: 2.80%    MaxErr: 3.79%
    Distance_FlipCode_ShiftCode:            592  ms    AvgErr: 2.39%    MaxErr: 4.30%
    Distance_Pearls:                        6998 ms    AvgErr: 2.37%    MaxErr:94.60%
    Distance_Pearls_2:                      4589 ms    AvgErr: 7.53%    MaxErr:98.20%
*/

using System;
using System.Runtime.CompilerServices;

namespace HypothenuseCalculations {
    class Program {

        static void ABSAndOrder(ref int a, ref int b) {
            if (a < 0) a = -a;
            if (b < 0) b = -b;

            if (a > b) {
                int c = a;
                a = b;
                b = c;
            }
        }

        static void ABSAndOrder(ref double a, ref double b) {
            if (a < 0) a = -a;
            if (b < 0) b = -b;

            if (a > b) {
                double c = a;
                a = b;
                b = c;
            }
        }

        static int Distance_Normal(int a, int b) {
            return (int)Math.Sqrt(Math.Pow(a, 2) + Math.Pow(b, 2));
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static int Distance_Sqrt2(int a, int b) {
            ABSAndOrder(ref a, ref b);
            return (int)((Math.Sqrt(2) - 1) * a) + b;
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static uint Distance_4142_Cast(int a, int b) {
            ABSAndOrder(ref a, ref b);
            return (uint)(0.414213562 * (uint)a + (uint)b);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static int Distance_4142(int a, int b) {
            ABSAndOrder(ref a, ref b);
            return (int)(0.414213562 * a + b);
        }
        
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static int Distance_4142_Int_Only(int a, int b) {
            ABSAndOrder(ref a, ref b);
            return 4142 * a / 10000 + b;
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static int Distance_4142_Int_Only_Manual_Inlined(int a, int b) {
            if (a < 0) a = -a;
            if (b < 0) b = -b;

            if (a > b) {
                int c = a;
                a = b;
                b = c;
            }
            return 4142 * a / 10000 + b;
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static int Distance_AShift1(int a, int b) {
            ABSAndOrder(ref a, ref b);
            return (a >> 1) + b;
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static int Distance_337(int a, int b) {
            ABSAndOrder(ref a, ref b);
            return (int)(b + 0.337 * a);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static int Distance_918(int a, int b) {
            ABSAndOrder(ref a, ref b);
            return (int)Math.Max(b, 0.918 * (b + (a >> 1)));
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static int Distance_428(int a, int b) {
            if (a == 0 && b == 0)
                return 0; //Otherwise we get some crazy number
            ABSAndOrder(ref a, ref b);
            return (int)(b + 0.428 * a * a / b);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static int Distance_FlipCode(int a, int b) {
            int approx;

            ABSAndOrder(ref a, ref b);

            approx = (b * 1007) + (a * 441);
            if (b < (a << 4))
                approx -= (b * 40);
            
            return ((approx + 512) >> 10);
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static int Distance_FlipCode_ShiftCode(int a, int b) {
            ABSAndOrder(ref a, ref b);

            return (((b << 8) + (b << 3) - (b << 4) - (b << 1) + (a << 7) - (a << 5) + (a << 3) - (a << 1)) >> 8);
        }

        /* Iterations   Accuracy
         *  2          6.5 digits
         *  3           20 digits
         *  4           62 digits
         * assuming a numeric type able to maintain that degree of accuracy in
         * the individual operations.

         * Transliterated from _More Programming Pearls, Confessions of a Coder_
         * by Jon Bentley, pg. 156.
         */
const int ITER = 3;

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static double Distance_Pearls(double a, double b) {
            ABSAndOrder(ref a, ref b);
            double R;
            for (int i = 0; i < ITER; i++) {
                R = b / a;
                R = R * R;
                R = R / (4.0 + R);
                a = a + 2.0 * R * a;
                b = b * R;
            }
            return a;
        }

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static double Distance_Pearls_2(double a, double b) {
            ABSAndOrder(ref a, ref b);
            double R;
            for (int i = 0; i < 2; i++) {
                R = b / a;
                R = R * R;
                R = R / (4.0 + R);
                a = a + 2.0 * R * a;
                b = b * R;
            }
            return a;
        }

        static void Main(string[] args) {

            int iterations = 100000000;
            
            Console.WriteLine("Starting Speed Test for " + iterations + " iterations");
            var watch = System.Diagnostics.Stopwatch.StartNew();
            for (int i = 0; i < iterations; i += 1) {
                Distance_Normal(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_Normal: " + watch.ElapsedMilliseconds + "ms");*/

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_Sqrt2(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_Sqrt2: " + watch.ElapsedMilliseconds + "ms");

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_4142(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_4142: " + watch.ElapsedMilliseconds + "ms");

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_4142_Cast(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_4142_Cast: " + watch.ElapsedMilliseconds + "ms");

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_4142_Int_Only(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_4142_Int_Only: " + watch.ElapsedMilliseconds + "ms");

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_4142_Int_Only_Manual_Inlined(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_4142_Int_Only_Manual_Inlined: " + watch.ElapsedMilliseconds + "ms");

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_AShift1(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_AShift1: " + watch.ElapsedMilliseconds + "ms");

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_337(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_337: " + watch.ElapsedMilliseconds + "ms");

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_918(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_918: " + watch.ElapsedMilliseconds + "ms");

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_428(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_428: " + watch.ElapsedMilliseconds + "ms"); 

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_FlipCode(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_FlipCode: " + watch.ElapsedMilliseconds + "ms");

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_FlipCode_ShiftCode(777, 999);
            }
            watch.Stop();
            Console.WriteLine("Distance_FlipCode_ShiftCode: " + watch.ElapsedMilliseconds + "ms");

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_Pearls(777.0, 999.0);
            }
            watch.Stop();
            Console.WriteLine("Distance_Pearls: " + watch.ElapsedMilliseconds + "ms");

            watch.Restart();
            for (int i = 0; i < iterations; i += 1) {
                Distance_Pearls_2(777.0, 999.0);
            }
            watch.Stop();
            Console.WriteLine("Distance_Pearls_2: " + watch.ElapsedMilliseconds + "ms");
            
            SetUpControl();
            TestAccuracy("Distance_Sqrt2", CreateWrapper(Distance_Sqrt2));
            TestAccuracy("Distance_4142", CreateWrapper(Distance_4142));
            TestAccuracy("Distance_4142_Int_Only", CreateWrapper(Distance_4142_Int_Only));
            TestAccuracy("Distance_AShift1", CreateWrapper(Distance_AShift1));
            TestAccuracy("Distance_337", CreateWrapper(Distance_337));
            TestAccuracy("Distance_918", CreateWrapper(Distance_918));
            TestAccuracy("Distance_428", CreateWrapper(Distance_428));
            TestAccuracy("Distance_FlipCode", CreateWrapper(Distance_FlipCode));
            TestAccuracy("Distance_FlipCode_ShiftCode", CreateWrapper(Distance_FlipCode_ShiftCode));
            TestAccuracy("Distance_Pearls", CreateWrapperD(Distance_Pearls));
            TestAccuracy("Distance_Pearls_2", CreateWrapperD(Distance_Pearls_2));

            Console.ReadLine();
        }

        static Func<int, int, double> CreateWrapper(Func<int, int, int> func) {
            return (a, b) => { return func(a, b); };
        }

        static Func<int, int, double> CreateWrapperD(Func<double, double, double> func) {
            return (a, b) => { return func(a, b); };
        }

        static Func<int, int, double> CreateWrapperU(Func<int, int, uint> func) {
            return (a, b) => { return func(a, b); };
        }

        const int cWidth = 1000;
        const int cHeight = 1000;
        static double[,] sAccuracyTableControl = new double[cWidth, cHeight];
        static void SetUpControl() {
            Console.WriteLine("\nSetting Up Control For accuracy testing");
            for (int x = 0; x < cWidth; x+=1) {
                for (int y = 0; y < cHeight; y+=1) {
                    int xx = x - cWidth / 2;
                    int yy = y - cHeight / 2;
                    sAccuracyTableControl[x, y] = Distance_Normal(xx*1000, yy*1000);
                }
            }
        }

        static void TestAccuracy (string name, Func<int, int, double> distance) {
            int width = 1000;
            int height = 1000;
            double[,] accuracyTable = new double[width, height];
            double MaxError = 0.0;
            UInt64 AvgError = 0;
            for (int x = 0; x < cWidth; x += 1) {
                for (int y = 0; y < cHeight; y += 1) {
                    int xx = x - cWidth / 2;
                    int yy = y - cHeight / 2;
                    if (xx == 0 && yy == 0)
                        continue;
                    double distanceControl = sAccuracyTableControl[x, y];
                    double distanceTest = distance(xx*1000, yy*1000);
                    double error = Math.Abs((distanceTest - distanceControl) / Math.Abs(distanceControl));
                    if (error > MaxError)
                        MaxError = error;
                    AvgError += (UInt64)(error * 10000);
                    accuracyTable[x, y] = error;
                }
            }
            AvgError = AvgError / (cWidth * cHeight);
            int iMaxError = (int)(MaxError * 10000);
            Console.WriteLine(name + " - AvgError: " + AvgError/100 + "." + AvgError%100 + "% MaxError: " + iMaxError/100 + "." + iMaxError%100 + "%");
        }
    }
}

https://atoms.alife.co.uk/sqrt/index.html - has more algorithms I would like to add

Also, add
i = 0x2A555555 + ( i / 3 ); from - https://timmmm.github.io/fast-inverse-square-root/

https://en.wikipedia.org/wiki/Alpha_max_plus_beta_min_algorithm