diff --git a/src/Numerics.Tests/RootFindingTests/FindRootsTest.cs b/src/Numerics.Tests/RootFindingTests/FindRootsTest.cs index 724269435..da794c670 100644 --- a/src/Numerics.Tests/RootFindingTests/FindRootsTest.cs +++ b/src/Numerics.Tests/RootFindingTests/FindRootsTest.cs @@ -30,6 +30,7 @@ using System; using Complex = System.Numerics.Complex; using NUnit.Framework; +using System.Diagnostics; namespace MathNet.Numerics.Tests.RootFindingTests { @@ -40,7 +41,7 @@ internal class FindRootsTest public void MultipleRoots() { // Roots at -2, 2 - Func f1 = x => x*x - 4; + Func f1 = x => x * x - 4; Assert.AreEqual(0, f1(FindRoots.OfFunction(f1, 0, 5, 1e-14))); Assert.AreEqual(-2, FindRoots.OfFunction(f1, -5, -1, 1e-14)); Assert.AreEqual(2, FindRoots.OfFunction(f1, 1, 4, 1e-14)); @@ -49,7 +50,7 @@ public void MultipleRoots() Assert.AreEqual(2, FindRoots.OfFunction(x => -f1(x), 1, 4, 1e-14)); // Roots at 3, 4 - Func f2 = x => (x - 3)*(x - 4); + Func f2 = x => (x - 3) * (x - 4); Assert.AreEqual(0, f2(FindRoots.OfFunction(f2, 3.5, 5, 1e-14)), 1e-14); Assert.AreEqual(3, FindRoots.OfFunction(f2, -5, 3.5, 1e-14)); Assert.AreEqual(4, FindRoots.OfFunction(f2, 3.2, 5, 1e-14)); @@ -60,7 +61,7 @@ public void MultipleRoots() [Test] public void LocalMinima() { - Func f1 = x => x*x*x - 2*x + 2; + Func f1 = x => x * x * x - 2 * x + 2; Assert.AreEqual(0, f1(FindRoots.OfFunction(f1, -5, 5, 1e-14)), 1e-14); Assert.AreEqual(0, f1(FindRoots.OfFunction(f1, -2, 4, 1e-14)), 1e-14); } @@ -68,7 +69,7 @@ public void LocalMinima() [Test] public void NoRoot() { - Func f1 = x => x*x + 4; + Func f1 = x => x * x + 4; Assert.That(() => FindRoots.OfFunction(f1, -5, 5, 1e-14), Throws.TypeOf()); } @@ -76,7 +77,7 @@ public void NoRoot() public void Oneeq3() { // Test case from http://www.polymath-software.com/library/nle/Oneeq3.htm - Func f1 = T => Math.Exp(21000/T)/(T*T) - 1.11e11; + Func f1 = T => Math.Exp(21000 / T) / (T * T) - 1.11e11; double x = FindRoots.OfFunction(f1, 550, 560, 1e-12); Assert.AreEqual(551.773822885233, x, 1e-5); Assert.AreEqual(0, f1(x), 1e-2); @@ -92,7 +93,7 @@ public void Oneeq5() const double BETA = 2.298E-3; const double GAMMA = 0.283E-6; const double DH = -57798.0; - return DH + TR*(ALPHA + TR*(BETA/2 + TR*GAMMA/3)) - 298.0*(ALPHA + 298.0*(BETA/2 + 298.0*GAMMA/3)); + return DH + TR * (ALPHA + TR * (BETA / 2 + TR * GAMMA / 3)) - 298.0 * (ALPHA + 298.0 * (BETA / 2 + 298.0 * GAMMA / 3)); }; double x = FindRoots.OfFunction(f1, 3000, 5000, 1e-10); @@ -114,11 +115,11 @@ public void Oneeq6a() const double A = 0.01855; const double B = -0.01587; const double P = 100.0; - const double Beta = R*T*B0 - A0 - R*C/(T*T); - const double Gama = -R*T*B0*B + A0*A - R*C*B0/(T*T); - const double Delta = R*B0*B*C/(T*T); + const double Beta = R * T * B0 - A0 - R * C / (T * T); + const double Gama = -R * T * B0 * B + A0 * A - R * C * B0 / (T * T); + const double Delta = R * B0 * B * C / (T * T); - return R*T/V + Beta/(V*V) + Gama/(V*V*V) + Delta/(V*V*V*V) - P; + return R * T / V + Beta / (V * V) + Gama / (V * V * V) + Delta / (V * V * V * V) - P; }; double x = FindRoots.OfFunction(f1, 0.1, 1); @@ -130,7 +131,7 @@ public void Oneeq6a() public void Oneeq7() { // Test case from http://www.polymath-software.com/library/nle/Oneeq7.htm - Func f1 = x => x/(1 - x) - 5*Math.Log(0.4*(1 - x)/(0.4 - 0.5*x)) + 4.45977; + Func f1 = x => x / (1 - x) - 5 * Math.Log(0.4 * (1 - x) / (0.4 - 0.5 * x)) + 4.45977; double r = FindRoots.OfFunction(f1, 0, 0.79, 1e-10); Assert.AreEqual(0.757396293891, r, 1e-6); Assert.AreEqual(0, f1(r), 1e-6); @@ -146,7 +147,7 @@ public void Oneeq8() const double b = 40; const double c = 200; - return a*v*v + b*Math.Pow(v, 7/4) - c; + return a * v * v + b * Math.Pow(v, 7 / 4) - c; }; double x = FindRoots.OfFunction(f1, 0.01, 1); @@ -158,7 +159,7 @@ public void Oneeq8() public void StackOverflow39935588() { // Roots at -2, 2 - Func f1 = x => (x - 3.0)*(x - 4.0); + Func f1 = x => (x - 3.0) * (x - 4.0); Assert.AreEqual(3.0, FindRoots.OfFunction(f1, -2.0, 3.5), 1e-10); Assert.AreEqual(4.0, FindRoots.OfFunction(f1, 3.5, 5.5), 1e-10); Assert.AreEqual(0.0, f1(FindRoots.OfFunction(f1, -2.0, 5.5, 1e-14)), 1e-14); @@ -177,7 +178,7 @@ void AssertComplexEqual(Complex expected, Complex actual, double delta) public void QuadraticExpanded(double u, double v, double t) { // t*(x+u)*(x+v) = t*u*v + t*(u+v)*x + t*x^2 - double c = t*u*v, b = t*(u + v), a = t; + double c = t * u * v, b = t * (u + v), a = t; var x = FindRoots.Quadratic(c, b, a); Complex x1 = x.Item1, x2 = x.Item2; @@ -187,8 +188,8 @@ public void QuadraticExpanded(double u, double v, double t) || x1.AlmostEqualRelative(r2, 1e-14) && x2.AlmostEqualRelative(r1, 1e-14)); // Verify they really are roots - AssertComplexEqual(Complex.Zero, c + b*x1 + a*x1*x1, 1e-14); - AssertComplexEqual(Complex.Zero, c + b*x2 + a*x2*x2, 1e-14); + AssertComplexEqual(Complex.Zero, c + b * x1 + a * x1 * x1, 1e-14); + AssertComplexEqual(Complex.Zero, c + b * x2 + a * x2 * x2, 1e-14); } [TestCase(1d, 1d, 1d, -0.5, -0.866025403784439, -0.5, 0.866025403784439)] @@ -207,8 +208,101 @@ public void QuadraticExpected(double c, double b, double a, double x1R, double x || x1.AlmostEqualRelative(r2, 1e-14) && x2.AlmostEqualRelative(r1, 1e-14)); // Verify they really are roots - AssertComplexEqual(Complex.Zero, c + b*x1 + a*x1*x1, 1e-14); - AssertComplexEqual(Complex.Zero, c + b*x2 + a*x2*x2, 1e-14); + AssertComplexEqual(Complex.Zero, c + b * x1 + a * x1 * x1, 1e-14); + AssertComplexEqual(Complex.Zero, c + b * x2 + a * x2 * x2, 1e-14); + } + + public static (Complex, Complex) FasterQuadratic(double c, double b, double a) + { + Complex x1, x2; + if (a == 0.0) + { + if (b == 0.0) + { + x1 = new Complex(double.NaN, double.NaN); // Complex.NaN; + x2 = x1; + } + else + { + x1 = new Complex(-c / b, 0.0); + x2 = x1; + } + } + else + { + a = 1.0 / a; + b = -0.5 * b * a; + c = c * a; + double delta = b * b - c; + if (delta < 0.0) + { + double sqrtDelta = Math.Sqrt(-delta); + x1 = new Complex(b, sqrtDelta); + x2 = new Complex(b, -sqrtDelta); + } + else + { + double sqrtDelta = Math.Sqrt(delta); + x1 = new Complex(b + sqrtDelta, 0.0); + x2 = new Complex(b - sqrtDelta, 0.0); + } + } + return (x1, x2); + } + + [TestCase(70_000_000)] + public void TestQuadraticSpeed(int N) + { + Complex x1, x2; + double sum = 0.0; + double a, b, c; + Stopwatch stopwatch = new Stopwatch(); + for (int i = 11; i < N; i++) + { + a = -0.1 * i + 1.0; + b = 0.3464 * i + 5.0; + c = -0.3 * i + 7.0; + (x1, x2) = FindRoots.Quadratic(c, b, a); + sum += x1.Real + x2.Imaginary; + } + Console.WriteLine($"sum={sum}"); + stopwatch.Stop(); + Console.WriteLine($"FindRoots.Quadratic time: {stopwatch.ElapsedMilliseconds * 0.001:F3}s"); + + sum = 0.0; + stopwatch.Restart(); + for (int i = 11; i < N; i++) + { + a = -0.1 * i + 1.0; + b = 0.3464 * i + 5.0; + c = -0.3 * i + 7.0; + (x1, x2) = FasterQuadratic(c, b, a); + sum += x1.Real + x2.Imaginary; + } + Console.WriteLine($"sum={sum}"); + stopwatch.Stop(); + Console.WriteLine($"FasterQuadratic time: {stopwatch.ElapsedMilliseconds * 0.001:F3}s"); + } + + [TestCase(7_000_000)] + public void TestFasterQuadraticCorrectness(int N) + { + Complex x1, x2, x3, x4; + double a, b, c, b_a, c_a; + for (int i = 11; i < N; i++) + { + a = -0.1 * i + 1.0; + b = 0.3464 * i + 5.0; + c = -0.3 * i + 7.0; + b_a = -b / a; + c_a = c / a; + (x1, x2) = FasterQuadratic(c, b, a); + (x3, x4) = FindRoots.Quadratic(c, b, a); + AssertComplexEqual(x1 + x2, b_a, 1e-14); + AssertComplexEqual(x1 * x2, c_a, 1e-14); + AssertComplexEqual(x3 + x4, x1 + x2, 1e-14); + AssertComplexEqual(x3 * x4, x1 * x2, 1e-14); + } } } } diff --git a/src/Numerics/FindRoots.cs b/src/Numerics/FindRoots.cs index 0b91930d1..2ca4302c5 100644 --- a/src/Numerics/FindRoots.cs +++ b/src/Numerics/FindRoots.cs @@ -43,7 +43,7 @@ public static class FindRoots /// Maximum number of iterations. Example: 100. public static double OfFunction(Func f, double lowerBound, double upperBound, double accuracy = 1e-8, int maxIterations = 100) { - if (!ZeroCrossingBracketing.ExpandReduce(f, ref lowerBound, ref upperBound, 1.6, maxIterations, maxIterations*10)) + if (!ZeroCrossingBracketing.ExpandReduce(f, ref lowerBound, ref upperBound, 1.6, maxIterations, maxIterations * 10)) { throw new NonConvergenceException("The algorithm has failed, exceeded the number of iterations allowed or there is no root within the provided bounds."); } @@ -86,17 +86,40 @@ public static double OfFunctionDerivative(Func f, Func public static (Complex, Complex) Quadratic(double c, double b, double a) { - if (b == 0d) + Complex x1, x2; + if (a == 0.0) { - var t = new Complex(-c/a, 0d).SquareRoot(); - return (t, -t); + if (b == 0.0) + { + x1 = Complex.Zero / Complex.Zero; // Complex.NaN; + x2 = x1; + } + else + { + x1 = new Complex(-c / b, 0.0); + x2 = x1; + } } - - var q = b > 0d - ? -0.5*(b + new Complex(b*b - 4*a*c, 0d).SquareRoot()) - : -0.5*(b - new Complex(b*b - 4*a*c, 0d).SquareRoot()); - - return (q/a, c/q); + else + { + a = 1.0 / a; + b = -0.5 * b * a; + c = c * a; + double delta = b * b - c; + if (delta < 0.0) + { + double sqrtDelta = Math.Sqrt(-delta); + x1 = new Complex(b, sqrtDelta); + x2 = new Complex(b, -sqrtDelta); + } + else + { + double sqrtDelta = Math.Sqrt(delta); + x1 = new Complex(b + sqrtDelta, 0.0); + x2 = new Complex(b - sqrtDelta, 0.0); + } + } + return (x1, x2); } /// @@ -143,16 +166,16 @@ public static double[] ChebychevPolynomialFirstKind(int degree, double intervalB } // transform to map to [-1..1] interval - double location = 0.5*(intervalBegin + intervalEnd); - double scale = 0.5*(intervalEnd - intervalBegin); + double location = 0.5 * (intervalBegin + intervalEnd); + double scale = 0.5 * (intervalEnd - intervalBegin); // evaluate first kind chebychev nodes - double angleFactor = Constants.Pi/(2*degree); + double angleFactor = Constants.Pi / (2 * degree); var samples = new double[degree]; for (int i = 0; i < samples.Length; i++) { - samples[i] = location + scale*Math.Cos(((2*i) + 1)*angleFactor); + samples[i] = location + scale * Math.Cos(((2 * i) + 1) * angleFactor); } return samples; } @@ -172,16 +195,16 @@ public static double[] ChebychevPolynomialSecondKind(int degree, double interval } // transform to map to [-1..1] interval - double location = 0.5*(intervalBegin + intervalEnd); - double scale = 0.5*(intervalEnd - intervalBegin); + double location = 0.5 * (intervalBegin + intervalEnd); + double scale = 0.5 * (intervalEnd - intervalBegin); // evaluate second kind chebychev nodes - double angleFactor = Constants.Pi/(degree + 1); + double angleFactor = Constants.Pi / (degree + 1); var samples = new double[degree]; for (int i = 0; i < samples.Length; i++) { - samples[i] = location + scale*Math.Cos((i + 1)*angleFactor); + samples[i] = location + scale * Math.Cos((i + 1) * angleFactor); } return samples; } diff --git a/src/Numerics/SpecialFunctions/Factorial.cs b/src/Numerics/SpecialFunctions/Factorial.cs index 893001f54..6fbd7e5d9 100644 --- a/src/Numerics/SpecialFunctions/Factorial.cs +++ b/src/Numerics/SpecialFunctions/Factorial.cs @@ -39,6 +39,52 @@ public partial class SpecialFunctions { 1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200, 1307674368000, 20922789888000, 355687428096000, 6402373705728000, 1.21645100408832E+17, 2.43290200817664E+18, 5.109094217170944E+19, 1.1240007277776077E+21, 2.5852016738884978E+22, 6.2044840173323941E+23, 1.5511210043330986E+25, 4.0329146112660565E+26, 1.0888869450418352E+28, 3.0488834461171384E+29, 8.8417619937397008E+30, 2.6525285981219103E+32, 8.2228386541779224E+33, 2.6313083693369352E+35, 8.6833176188118859E+36, 2.9523279903960412E+38, 1.0333147966386144E+40, 3.7199332678990118E+41, 1.3763753091226343E+43, 5.2302261746660104E+44, 2.0397882081197442E+46, 8.1591528324789768E+47, 3.3452526613163803E+49, 1.4050061177528798E+51, 6.0415263063373834E+52, 2.6582715747884485E+54, 1.1962222086548019E+56, 5.5026221598120885E+57, 2.5862324151116818E+59, 1.2413915592536073E+61, 6.0828186403426752E+62, 3.0414093201713376E+64, 1.5511187532873822E+66, 8.0658175170943877E+67, 4.2748832840600255E+69, 2.3084369733924138E+71, 1.2696403353658276E+73, 7.1099858780486348E+74, 4.0526919504877221E+76, 2.3505613312828789E+78, 1.3868311854568986E+80, 8.3209871127413916E+81, 5.0758021387722484E+83, 3.1469973260387939E+85, 1.9826083154044401E+87, 1.2688693218588417E+89, 8.2476505920824715E+90, 5.4434493907744307E+92, 3.6471110918188683E+94, 2.4800355424368305E+96, 1.711224524281413E+98, 1.197857166996989E+100, 8.5047858856786218E+101, 6.1234458376886077E+103, 4.4701154615126834E+105, 3.3078854415193856E+107, 2.4809140811395391E+109, 1.8854947016660498E+111, 1.4518309202828584E+113, 1.1324281178206295E+115, 8.9461821307829729E+116, 7.1569457046263779E+118, 5.7971260207473655E+120, 4.7536433370128398E+122, 3.9455239697206569E+124, 3.314240134565352E+126, 2.8171041143805494E+128, 2.4227095383672724E+130, 2.1077572983795269E+132, 1.8548264225739836E+134, 1.6507955160908452E+136, 1.4857159644817607E+138, 1.3520015276784023E+140, 1.24384140546413E+142, 1.1567725070816409E+144, 1.0873661566567424E+146, 1.0329978488239052E+148, 9.916779348709491E+149, 9.6192759682482062E+151, 9.426890448883242E+153, 9.3326215443944096E+155, 9.3326215443944102E+157, 9.4259477598383536E+159, 9.6144667150351211E+161, 9.9029007164861754E+163, 1.0299016745145622E+166, 1.0813967582402903E+168, 1.1462805637347078E+170, 1.2265202031961373E+172, 1.3246418194518284E+174, 1.4438595832024928E+176, 1.5882455415227421E+178, 1.7629525510902437E+180, 1.9745068572210728E+182, 2.2311927486598123E+184, 2.5435597334721862E+186, 2.9250936934930141E+188, 3.3931086844518965E+190, 3.969937160808719E+192, 4.6845258497542883E+194, 5.5745857612076033E+196, 6.6895029134491239E+198, 8.09429852527344E+200, 9.8750442008335976E+202, 1.2146304367025325E+205, 1.5061417415111404E+207, 1.8826771768889254E+209, 2.3721732428800459E+211, 3.0126600184576582E+213, 3.8562048236258025E+215, 4.9745042224772855E+217, 6.4668554892204716E+219, 8.4715806908788174E+221, 1.1182486511960039E+224, 1.4872707060906852E+226, 1.9929427461615181E+228, 2.6904727073180495E+230, 3.6590428819525472E+232, 5.0128887482749898E+234, 6.9177864726194859E+236, 9.6157231969410859E+238, 1.346201247571752E+241, 1.8981437590761701E+243, 2.6953641378881614E+245, 3.8543707171800706E+247, 5.5502938327393013E+249, 8.0479260574719866E+251, 1.1749972043909099E+254, 1.7272458904546376E+256, 2.5563239178728637E+258, 3.8089226376305671E+260, 5.7133839564458505E+262, 8.6272097742332346E+264, 1.3113358856834518E+267, 2.0063439050956811E+269, 3.0897696138473489E+271, 4.7891429014633912E+273, 7.4710629262828905E+275, 1.1729568794264138E+278, 1.8532718694937338E+280, 2.9467022724950369E+282, 4.714723635992059E+284, 7.5907050539472148E+286, 1.2296942187394488E+289, 2.0044015765453015E+291, 3.2872185855342945E+293, 5.423910666131586E+295, 9.0036917057784329E+297, 1.5036165148649983E+300, 2.5260757449731969E+302, 4.2690680090047027E+304, 7.257415615307994E+306 }; + static readonly double[] _factorialLnCache = new double[] + { + 0.0, 0.0,0.6931471805599453094, 1.791759469228055001, + 3.17805383034794562, 4.787491742782045994, 6.579251212010100995, 8.5251613610654143, + 10.60460290274525023, 12.80182748008146961, 15.1044125730755153, 17.50230784587388584, + 19.98721449566188615, 22.55216385312342289, 25.1912211827386815, 27.89927138384089157, + 30.6718601060806728, 33.50507345013688888, 36.39544520803305358, 39.33988418719949404, + 42.33561646075348503, 45.38013889847690803, 48.47118135183522388, 51.60667556776437357, + 54.78472939811231919, 58.00360522298051994, 61.26170176100200198, 64.55753862700633106, + 67.88974313718153498, 71.25703896716800901, 74.65823634883016439, 78.09222355331531063, + 81.55795945611503718, 85.05446701758151741, 88.5808275421976788, 92.13617560368709248, + 95.71969454214320248, 99.33061245478742693, 102.9681986145138127, 106.6317602606434591, + 110.3206397147573954, 114.0342117814617032, 117.7718813997450715, 121.533081515438634, + 125.3172711493568951, 129.1239336391272149, 132.9525750356163099, 136.8027226373263685, + 140.6739236482342594, 144.565743946344886, 148.4777669517730321, 152.4095925844973578, + 156.3608363030787852, 160.331128216630907, 164.3201122631951814, 168.3274454484276523, + 172.3527971391628016, 176.3958484069973517, 180.456291417543771, 184.5338288614494905, + 188.6281734236715912, 192.7390472878449024, 196.866181672889994, 201.0093163992815267, + 205.1681994826411985, 209.3425867525368356, 213.5322414945632612, 217.7369341139542272, + 221.956441819130334, 226.1905483237275933, 230.4390435657769523, 234.7017234428182677, + 238.978389561834323, 243.2688490029827142, 247.5729140961868839, 251.8904022097231944, + 256.2211355500095255, 260.5649409718632093, 264.921649798552801, 269.2910976510198225, + 273.6731242856937042, 278.0675734403661429, 282.474292687630396, 286.893133295426994, + 291.3239500942703076, 295.766601350760624, 300.2209486470141318, 304.6868567656687155, + 309.1641935801469219, 313.6528299498790618, 318.1526396202093268, 322.6634991267261769, + 327.1852877037752172, 331.7178871969284731, 336.261181979198477, 340.8150588707990179, + 345.3794070622668541, 349.9541180407702369, 354.5390855194408088, 359.1342053695753988, + 363.7393755555634901, 368.3544960724047496, 372.9794688856890207, 377.6141978739186564, + 382.2585887730600291, 386.9125491232175525, 391.5759882173296196, 396.2488170517915258, + 400.9309482789157455, 405.6222961611448892, 410.3227765269373054, 415.0323067282496396, + 419.7508055995447341, 424.4781934182570747, 429.2143918666515701, 433.9593239950148202, + 438.7129141861211848, 443.475088120918941, 448.2457727453846057, 453.0248962384961351, + 457.8123879812781811, 462.6081785268749222, 467.4121995716081787, 472.2243839269805962, + 477.0446654925856331, 481.8729792298879342, 486.7092611368394122, 491.5534482232980035, + 496.4054784872176207, 501.2652908915792928, 506.1328253420348752, 511.0080226652360267, + 515.8908245878223976, 520.7811737160441514, 525.6790135159950627, 530.5842882944334922, + 535.4969431801695442, 540.4169241059976691, 545.3441777911548738, 550.2786517242855656, + 555.2202941468948698, 560.1690540372730381, 565.1248810948742989, 570.0877257251342061, + 575.0575390247102068, 580.0342727671307812, 585.0178793888391176, 590.0083119756178539, + 595.005524249381969, 600.0094705553274281, 605.0201058494236839, 610.0373856862386082, + 615.0612662070848846, 620.09170412847732, 625.1286567308909492, 630.1720818478101958, + 635.2219378550597329, 640.2781836604080409, 645.3407786934350077, 650.4096828956552392, + 655.4848567108890662, 660.5662610758735292, 665.6538574111059132, 670.7476076119126756, + 675.847474039736874, 680.9534195136374546, 686.0654073019939978, 691.183401114410753, + 696.3073650938140119, 701.4372638087370854, 706.5730622457873471 + }; /// /// Computes the factorial function x -> x! of an integer number > 0. The function can represent all number up @@ -102,14 +148,9 @@ public static double FactorialLn(int x) throw new ArgumentOutOfRangeException(nameof(x), "Value must be positive (and not zero)."); } - if (x <= 1) - { - return 0d; - } - - if (x < _factorialCache.Length) + if (x < _factorialLnCache.Length) { - return Math.Log(_factorialCache[x]); + return _factorialLnCache[x]; } return GammaLn(x + 1.0);