From cbc1ef22801a7ea89581ec26fbe9b4f9fb610d2d Mon Sep 17 00:00:00 2001
From: Febin <25066330+febkor@users.noreply.github.com>
Date: Fri, 5 Mar 2021 09:45:08 +0200
Subject: [PATCH 1/2] Precision: pre-compute negative powers of 10

---
 src/Numerics/Interpolation/CubicSpline.cs |  1 -
 src/Numerics/Precision.Equality.cs        | 22 +++++++++++++++++-----
 2 files changed, 17 insertions(+), 6 deletions(-)

diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs
index d2344901..0500c37e 100644
--- a/src/Numerics/Interpolation/CubicSpline.cs
+++ b/src/Numerics/Interpolation/CubicSpline.cs
@@ -237,7 +237,6 @@ namespace MathNet.Numerics.Interpolation
 
             var dd = new double[x.Length];
             var hPrev = x[1] - x[0];
-            // This check is quite costly as it usually involves a Math.Pow().
             var mPrevIs0 = m[0].AlmostEqual(0.0);
 
             for (var i = 1; i < x.Length - 1; ++i)
diff --git a/src/Numerics/Precision.Equality.cs b/src/Numerics/Precision.Equality.cs
index 8f4af6cc..4d64bd41 100644
--- a/src/Numerics/Precision.Equality.cs
+++ b/src/Numerics/Precision.Equality.cs
@@ -34,8 +34,6 @@ using Complex = System.Numerics.Complex;
 
 namespace MathNet.Numerics
 {
-    // TODO PERF: Cache/Precompute 10^x terms
-
     public static partial class Precision
     {
         /// <summary>
@@ -355,7 +353,7 @@ namespace MathNet.Numerics
             // 10^(-numberOfDecimalPlaces). We divide by two so that we have half the range
             // on each side of the numbers, e.g. if decimalPlaces == 2,
             // then 0.01 will equal between 0.005 and 0.015, but not 0.02 and not 0.00
-            return Math.Abs(diff) < Math.Pow(10, -decimalPlaces) / 2d;
+            return Math.Abs(diff) < Pow10(-decimalPlaces) * 0.5;
         }
 
         /// <summary>
@@ -431,7 +429,7 @@ namespace MathNet.Numerics
                 // 10^(-numberOfDecimalPlaces). We divide by two so that we have half the range
                 // on each side of the numbers, e.g. if decimalPlaces == 2,
                 // then 0.01 will equal between 0.005 and 0.015, but not 0.02 and not 0.00
-                return Math.Abs(diff) < Math.Pow(10, -decimalPlaces) / 2d;
+                return Math.Abs(diff) < Pow10(-decimalPlaces) * 0.5;
             }
 
             // If the magnitudes of the two numbers are equal to within one magnitude the numbers could potentially be equal
@@ -447,7 +445,7 @@ namespace MathNet.Numerics
             // 10^(-numberOfDecimalPlaces). We divide by two so that we have half the range
             // on each side of the numbers, e.g. if decimalPlaces == 2,
             // then 0.01 will equal between 0.00995 and 0.01005, but not 0.0015 and not 0.0095
-            return Math.Abs(diff) < Math.Pow(10, magnitudeOfMax - decimalPlaces) / 2d;
+            return Math.Abs(diff) < Pow10(magnitudeOfMax - decimalPlaces) * 0.5;
         }
 
         /// <summary>
@@ -1041,5 +1039,19 @@ namespace MathNet.Numerics
         {
             return AlmostEqualNormRelative(a.L2Norm(), b.L2Norm(), (a - b).L2Norm(), decimalPlaces);
         }
+
+        private static readonly double[] NegativePowersOf10 = new double[]
+        {
+            1, 0.1, 0.01, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9,
+            1e-10, 1e-11, 1e-12, 1e-13, 1e-14, 1e-15, 1e-16,
+            1e-17, 1e-18, 1e-19, 1e-20
+        };
+
+        private static double Pow10(int y)
+        {
+            return -NegativePowersOf10.Length < y && y <= 0
+               ? NegativePowersOf10[-y]
+               : Math.Pow(10.0, y);
+        }
     }
 }

From 79f29c450a68ed2d8e1fd52e48c346e17cb35087 Mon Sep 17 00:00:00 2001
From: Febin <25066330+febkor@users.noreply.github.com>
Date: Fri, 5 Mar 2021 09:46:48 +0200
Subject: [PATCH 2/2] Distributions: minor simplification to Cauchy for
 performance

---
 src/Numerics/Distributions/Cauchy.cs | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/src/Numerics/Distributions/Cauchy.cs b/src/Numerics/Distributions/Cauchy.cs
index c8af2860..f6d9b1bd 100644
--- a/src/Numerics/Distributions/Cauchy.cs
+++ b/src/Numerics/Distributions/Cauchy.cs
@@ -179,7 +179,8 @@ namespace MathNet.Numerics.Distributions
         /// <seealso cref="PDF"/>
         public double Density(double x)
         {
-            return 1.0/(Constants.Pi*_scale*(1.0 + (((x - _location)/_scale)*((x - _location)/_scale))));
+            var z = (x - _location)/_scale;
+            return 1.0/(Constants.Pi*_scale*(1.0 + z * z));
         }
 
         /// <summary>
@@ -190,7 +191,8 @@ namespace MathNet.Numerics.Distributions
         /// <seealso cref="PDFLn"/>
         public double DensityLn(double x)
         {
-            return -Math.Log(Constants.Pi*_scale*(1.0 + (((x - _location)/_scale)*((x - _location)/_scale))));
+            var z = (x - _location)/_scale;
+            return -Math.Log(Constants.Pi*_scale*(1.0 +  z * z));
         }
 
         /// <summary>
@@ -283,9 +285,9 @@ namespace MathNet.Numerics.Distributions
                 throw new ArgumentException("Invalid parametrization for the distribution.");
             }
 
-            return 1.0/(Constants.Pi*scale*(1.0 + (((x - location)/scale)*((x - location)/scale))));
+            var z = (x - location)/scale;
+            return 1.0/(Constants.Pi*scale*(1.0 + z * z));
         }
-
         /// <summary>
         /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
         /// </summary>
@@ -301,7 +303,8 @@ namespace MathNet.Numerics.Distributions
                 throw new ArgumentException("Invalid parametrization for the distribution.");
             }
 
-            return -Math.Log(Constants.Pi*scale*(1.0 + (((x - location)/scale)*((x - location)/scale))));
+            var z = (x - location)/scale;
+            return -Math.Log(Constants.Pi*scale*(1.0 + z * z ));
         }
 
         /// <summary>