diff --git a/MathNet.Numerics.Portable.sln.DotSettings b/MathNet.Numerics.Portable.sln.DotSettings index 1cf7b657..7377cde6 100644 --- a/MathNet.Numerics.Portable.sln.DotSettings +++ b/MathNet.Numerics.Portable.sln.DotSettings @@ -12,4 +12,6 @@ False True True - False \ No newline at end of file + False + <data /> + <data><IncludeFilters /><ExcludeFilters /></data> \ No newline at end of file diff --git a/MathNet.Numerics.sln.DotSettings b/MathNet.Numerics.sln.DotSettings index 75f0e0a2..d7e15e37 100644 --- a/MathNet.Numerics.sln.DotSettings +++ b/MathNet.Numerics.sln.DotSettings @@ -14,4 +14,6 @@ False False True - False \ No newline at end of file + False + <data /> + <data><IncludeFilters /><ExcludeFilters /></data> \ No newline at end of file diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index c46a500a..d53df1bc 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -104,6 +104,9 @@ + + + diff --git a/src/Numerics/SpecialFunctions/Erf.cs b/src/Numerics/SpecialFunctions/Erf.cs index 4a62d00f..ca996993 100644 --- a/src/Numerics/SpecialFunctions/Erf.cs +++ b/src/Numerics/SpecialFunctions/Erf.cs @@ -204,7 +204,7 @@ namespace MathNet.Numerics double[] n = new[] { 0.00337916709551257388990745, -0.00073695653048167948530905, -0.374732337392919607868241, 0.0817442448733587196071743, -0.0421089319936548595203468, 0.0070165709512095756344528, -0.00495091255982435110337458, 0.000871646599037922480317225 }; double[] d = new[] { 1, -0.218088218087924645390535, 0.412542972725442099083918, -0.0841891147873106755410271, 0.0655338856400241519690695, -0.0120019604454941768171266, 0.00408165558926174048329689, -0.000615900721557769691924509 }; - result = (z * 1.125) + (z * EvaluatePolynomial(n, z) / EvaluatePolynomial(d, z)); + result = (z * 1.125) + (z * Evaluate.Polynomial(n, z) / Evaluate.Polynomial(d, z)); } } else if ((z < 110) || ((z < 110) && invert)) @@ -219,7 +219,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 5.582813374e-21 double[] n = new[] { -0.0361790390718262471360258, 0.292251883444882683221149, 0.281447041797604512774415, 0.125610208862766947294894, 0.0274135028268930549240776, 0.00250839672168065762786937 }; double[] d = new[] { 1, 1.8545005897903486499845, 1.43575803037831418074962, 0.582827658753036572454135, 0.124810476932949746447682, 0.0113724176546353285778481 }; - r = EvaluatePolynomial(n, z - 0.5) / EvaluatePolynomial(d, z - 0.5); + r = Evaluate.Polynomial(n, z - 0.5) / Evaluate.Polynomial(d, z - 0.5); b = 0.3440242112F; } else if (z < 1.25) @@ -227,7 +227,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 4.01854729e-21 double[] n = new[] { -0.0397876892611136856954425, 0.153165212467878293257683, 0.191260295600936245503129, 0.10276327061989304213645, 0.029637090615738836726027, 0.0046093486780275489468812, 0.000307607820348680180548455 }; double[] d = new[] { 1, 1.95520072987627704987886, 1.64762317199384860109595, 0.768238607022126250082483, 0.209793185936509782784315, 0.0319569316899913392596356, 0.00213363160895785378615014 }; - r = EvaluatePolynomial(n, z - 0.75) / EvaluatePolynomial(d, z - 0.75); + r = Evaluate.Polynomial(n, z - 0.75) / Evaluate.Polynomial(d, z - 0.75); b = 0.419990927F; } else if (z < 2.25) @@ -235,7 +235,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 2.866005373e-21 double[] n = new[] { -0.0300838560557949717328341, 0.0538578829844454508530552, 0.0726211541651914182692959, 0.0367628469888049348429018, 0.00964629015572527529605267, 0.00133453480075291076745275, 0.778087599782504251917881e-4 }; double[] d = new[] { 1, 1.75967098147167528287343, 1.32883571437961120556307, 0.552528596508757581287907, 0.133793056941332861912279, 0.0179509645176280768640766, 0.00104712440019937356634038, -0.106640381820357337177643e-7 }; - r = EvaluatePolynomial(n, z - 1.25) / EvaluatePolynomial(d, z - 1.25); + r = Evaluate.Polynomial(n, z - 1.25) / Evaluate.Polynomial(d, z - 1.25); b = 0.4898625016F; } else if (z < 3.5) @@ -243,7 +243,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 1.045355789e-21 double[] n = new[] { -0.0117907570137227847827732, 0.014262132090538809896674, 0.0202234435902960820020765, 0.00930668299990432009042239, 0.00213357802422065994322516, 0.00025022987386460102395382, 0.120534912219588189822126e-4 }; double[] d = new[] { 1, 1.50376225203620482047419, 0.965397786204462896346934, 0.339265230476796681555511, 0.0689740649541569716897427, 0.00771060262491768307365526, 0.000371421101531069302990367 }; - r = EvaluatePolynomial(n, z - 2.25) / EvaluatePolynomial(d, z - 2.25); + r = Evaluate.Polynomial(n, z - 2.25) / Evaluate.Polynomial(d, z - 2.25); b = 0.5317370892F; } else if (z < 5.25) @@ -251,7 +251,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 8.300028706e-22 double[] n = new[] { -0.00546954795538729307482955, 0.00404190278731707110245394, 0.0054963369553161170521356, 0.00212616472603945399437862, 0.000394984014495083900689956, 0.365565477064442377259271e-4, 0.135485897109932323253786e-5 }; double[] d = new[] { 1, 1.21019697773630784832251, 0.620914668221143886601045, 0.173038430661142762569515, 0.0276550813773432047594539, 0.00240625974424309709745382, 0.891811817251336577241006e-4, -0.465528836283382684461025e-11 }; - r = EvaluatePolynomial(n, z - 3.5) / EvaluatePolynomial(d, z - 3.5); + r = Evaluate.Polynomial(n, z - 3.5) / Evaluate.Polynomial(d, z - 3.5); b = 0.5489973426F; } else if (z < 8) @@ -259,7 +259,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 1.700157534e-21 double[] n = new[] { -0.00270722535905778347999196, 0.0013187563425029400461378, 0.00119925933261002333923989, 0.00027849619811344664248235, 0.267822988218331849989363e-4, 0.923043672315028197865066e-6 }; double[] d = new[] { 1, 0.814632808543141591118279, 0.268901665856299542168425, 0.0449877216103041118694989, 0.00381759663320248459168994, 0.000131571897888596914350697, 0.404815359675764138445257e-11 }; - r = EvaluatePolynomial(n, z - 5.25) / EvaluatePolynomial(d, z - 5.25); + r = Evaluate.Polynomial(n, z - 5.25) / Evaluate.Polynomial(d, z - 5.25); b = 0.5571740866F; } else if (z < 11.5) @@ -267,7 +267,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 3.002278011e-22 double[] n = new[] { -0.00109946720691742196814323, 0.000406425442750422675169153, 0.000274499489416900707787024, 0.465293770646659383436343e-4, 0.320955425395767463401993e-5, 0.778286018145020892261936e-7 }; double[] d = new[] { 1, 0.588173710611846046373373, 0.139363331289409746077541, 0.0166329340417083678763028, 0.00100023921310234908642639, 0.24254837521587225125068e-4 }; - r = EvaluatePolynomial(n, z - 8) / EvaluatePolynomial(d, z - 8); + r = Evaluate.Polynomial(n, z - 8) / Evaluate.Polynomial(d, z - 8); b = 0.5609807968F; } else if (z < 17) @@ -275,7 +275,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 6.741114695e-21 double[] n = new[] { -0.00056907993601094962855594, 0.000169498540373762264416984, 0.518472354581100890120501e-4, 0.382819312231928859704678e-5, 0.824989931281894431781794e-7 }; double[] d = new[] { 1, 0.339637250051139347430323, 0.043472647870310663055044, 0.00248549335224637114641629, 0.535633305337152900549536e-4, -0.117490944405459578783846e-12 }; - r = EvaluatePolynomial(n, z - 11.5) / EvaluatePolynomial(d, z - 11.5); + r = Evaluate.Polynomial(n, z - 11.5) / Evaluate.Polynomial(d, z - 11.5); b = 0.5626493692F; } else if (z < 24) @@ -283,7 +283,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 7.802346984e-22 double[] n = new[] { -0.000241313599483991337479091, 0.574224975202501512365975e-4, 0.115998962927383778460557e-4, 0.581762134402593739370875e-6, 0.853971555085673614607418e-8 }; double[] d = new[] { 1, 0.233044138299687841018015, 0.0204186940546440312625597, 0.000797185647564398289151125, 0.117019281670172327758019e-4 }; - r = EvaluatePolynomial(n, z - 17) / EvaluatePolynomial(d, z - 17); + r = Evaluate.Polynomial(n, z - 17) / Evaluate.Polynomial(d, z - 17); b = 0.5634598136F; } else if (z < 38) @@ -291,7 +291,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 2.414228989e-22 double[] n = new[] { -0.000146674699277760365803642, 0.162666552112280519955647e-4, 0.269116248509165239294897e-5, 0.979584479468091935086972e-7, 0.101994647625723465722285e-8 }; double[] d = new[] { 1, 0.165907812944847226546036, 0.0103361716191505884359634, 0.000286593026373868366935721, 0.298401570840900340874568e-5 }; - r = EvaluatePolynomial(n, z - 24) / EvaluatePolynomial(d, z - 24); + r = Evaluate.Polynomial(n, z - 24) / Evaluate.Polynomial(d, z - 24); b = 0.5638477802F; } else if (z < 60) @@ -299,7 +299,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 5.896543869e-24 double[] n = new[] { -0.583905797629771786720406e-4, 0.412510325105496173512992e-5, 0.431790922420250949096906e-6, 0.993365155590013193345569e-8, 0.653480510020104699270084e-10 }; double[] d = new[] { 1, 0.105077086072039915406159, 0.00414278428675475620830226, 0.726338754644523769144108e-4, 0.477818471047398785369849e-6 }; - r = EvaluatePolynomial(n, z - 38) / EvaluatePolynomial(d, z - 38); + r = Evaluate.Polynomial(n, z - 38) / Evaluate.Polynomial(d, z - 38); b = 0.5640528202F; } else if (z < 85) @@ -307,7 +307,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 3.080612264e-21 double[] n = new[] { -0.196457797609229579459841e-4, 0.157243887666800692441195e-5, 0.543902511192700878690335e-7, 0.317472492369117710852685e-9 }; double[] d = new[] { 1, 0.052803989240957632204885, 0.000926876069151753290378112, 0.541011723226630257077328e-5, 0.535093845803642394908747e-15 }; - r = EvaluatePolynomial(n, z - 60) / EvaluatePolynomial(d, z - 60); + r = Evaluate.Polynomial(n, z - 60) / Evaluate.Polynomial(d, z - 60); b = 0.5641309023F; } else @@ -315,7 +315,7 @@ namespace MathNet.Numerics // Worst case absolute error found: 8.094633491e-22 double[] n = new[] { -0.789224703978722689089794e-5, 0.622088451660986955124162e-6, 0.145728445676882396797184e-7, 0.603715505542715364529243e-10 }; double[] d = new[] { 1, 0.0375328846356293715248719, 0.000467919535974625308126054, 0.193847039275845656900547e-5 }; - r = EvaluatePolynomial(n, z - 85) / EvaluatePolynomial(d, z - 85); + r = Evaluate.Polynomial(n, z - 85) / Evaluate.Polynomial(d, z - 85); b = 0.5641584396F; } @@ -409,7 +409,7 @@ namespace MathNet.Numerics double[] P = new[] { -0.000508781949658280665617, -0.00836874819741736770379, 0.0334806625409744615033, -0.0126926147662974029034, -0.0365637971411762664006, 0.0219878681111168899165, 0.00822687874676915743155, -0.00538772965071242932965 }; double[] Q = new[] { 1, -0.970005043303290640362, -1.56574558234175846809, 1.56221558398423026363, 0.662328840472002992063, -0.71228902341542847553, -0.0527396382340099713954, 0.0795283687341571680018, -0.00233393759374190016776, 0.000886216390456424707504 }; double g = p * (p + 10); - double r = EvaluatePolynomial(P, p) / EvaluatePolynomial(Q, p); + double r = Evaluate.Polynomial(P, p) / Evaluate.Polynomial(Q, p); result = (g * Y) + (g * r); } else if (q >= 0.25) @@ -431,7 +431,7 @@ namespace MathNet.Numerics double[] Q = new[] { 1, 6.24264124854247537712, 3.9713437953343869095, -28.6608180499800029974, -20.1432634680485188801, 48.5609213108739935468, 10.8268667355460159008, -22.6436933413139721736, 1.72114765761200282724 }; double g = Math.Sqrt(-2 * Math.Log(q)); double xs = q - 0.25; - double r = EvaluatePolynomial(P, xs) / EvaluatePolynomial(Q, xs); + double r = Evaluate.Polynomial(P, xs) / Evaluate.Polynomial(Q, xs); result = g / (Y + r); } else @@ -463,7 +463,7 @@ namespace MathNet.Numerics double[] P = new[] { -0.131102781679951906451, -0.163794047193317060787, 0.117030156341995252019, 0.387079738972604337464, 0.337785538912035898924, 0.142869534408157156766, 0.0290157910005329060432, 0.00214558995388805277169, -0.679465575181126350155e-6, 0.285225331782217055858e-7, -0.681149956853776992068e-9 }; double[] Q = new[] { 1, 3.46625407242567245975, 5.38168345707006855425, 4.77846592945843778382, 2.59301921623620271374, 0.848854343457902036425, 0.152264338295331783612, 0.01105924229346489121 }; double xs = x - 1.125; - double R = EvaluatePolynomial(P, xs) / EvaluatePolynomial(Q, xs); + double R = Evaluate.Polynomial(P, xs) / Evaluate.Polynomial(Q, xs); result = (Y * x) + (R * x); } else if (x < 6) @@ -473,7 +473,7 @@ namespace MathNet.Numerics double[] P = new[] { -0.0350353787183177984712, -0.00222426529213447927281, 0.0185573306514231072324, 0.00950804701325919603619, 0.00187123492819559223345, 0.000157544617424960554631, 0.460469890584317994083e-5, -0.230404776911882601748e-9, 0.266339227425782031962e-11 }; double[] Q = new[] { 1, 1.3653349817554063097, 0.762059164553623404043, 0.220091105764131249824, 0.0341589143670947727934, 0.00263861676657015992959, 0.764675292302794483503e-4 }; double xs = x - 3; - double R = EvaluatePolynomial(P, xs) / EvaluatePolynomial(Q, xs); + double R = Evaluate.Polynomial(P, xs) / Evaluate.Polynomial(Q, xs); result = (Y * x) + (R * x); } else if (x < 18) @@ -483,7 +483,7 @@ namespace MathNet.Numerics double[] P = new[] { -0.0167431005076633737133, -0.00112951438745580278863, 0.00105628862152492910091, 0.000209386317487588078668, 0.149624783758342370182e-4, 0.449696789927706453732e-6, 0.462596163522878599135e-8, -0.281128735628831791805e-13, 0.99055709973310326855e-16 }; double[] Q = new[] { 1, 0.591429344886417493481, 0.138151865749083321638, 0.0160746087093676504695, 0.000964011807005165528527, 0.275335474764726041141e-4, 0.282243172016108031869e-6 }; double xs = x - 6; - double R = EvaluatePolynomial(P, xs) / EvaluatePolynomial(Q, xs); + double R = Evaluate.Polynomial(P, xs) / Evaluate.Polynomial(Q, xs); result = (Y * x) + (R * x); } else if (x < 44) @@ -493,7 +493,7 @@ namespace MathNet.Numerics double[] P = new[] { -0.0024978212791898131227, -0.779190719229053954292e-5, 0.254723037413027451751e-4, 0.162397777342510920873e-5, 0.396341011304801168516e-7, 0.411632831190944208473e-9, 0.145596286718675035587e-11, -0.116765012397184275695e-17 }; double[] Q = new[] { 1, 0.207123112214422517181, 0.0169410838120975906478, 0.000690538265622684595676, 0.145007359818232637924e-4, 0.144437756628144157666e-6, 0.509761276599778486139e-9 }; double xs = x - 18; - double R = EvaluatePolynomial(P, xs) / EvaluatePolynomial(Q, xs); + double R = Evaluate.Polynomial(P, xs) / Evaluate.Polynomial(Q, xs); result = (Y * x) + (R * x); } else @@ -503,31 +503,12 @@ namespace MathNet.Numerics double[] P = new[] { -0.000539042911019078575891, -0.28398759004727721098e-6, 0.899465114892291446442e-6, 0.229345859265920864296e-7, 0.225561444863500149219e-9, 0.947846627503022684216e-12, 0.135880130108924861008e-14, -0.348890393399948882918e-21 }; double[] Q = new[] { 1, 0.0845746234001899436914, 0.00282092984726264681981, 0.468292921940894236786e-4, 0.399968812193862100054e-6, 0.161809290887904476097e-8, 0.231558608310259605225e-11 }; double xs = x - 44; - double R = EvaluatePolynomial(P, xs) / EvaluatePolynomial(Q, xs); + double R = Evaluate.Polynomial(P, xs) / Evaluate.Polynomial(Q, xs); result = (Y * x) + (R * x); } } return s * result; } - - /// - /// A helper function to evaluate polynomials fast. - /// - /// The coefficients of the polynomial. - /// The location where to evaluate the polynomial at. - /// the evaluation of the polynomial. - private static double EvaluatePolynomial(double[] poly, double z) - { - int count = poly.Length; - double sum = poly[count - 1]; - for (int i = count - 2; i >= 0; --i) - { - sum *= z; - sum += poly[i]; - } - - return sum; - } } } diff --git a/src/Numerics/SpecialFunctions/Evaluate.cs b/src/Numerics/SpecialFunctions/Evaluate.cs new file mode 100644 index 00000000..096ca82e --- /dev/null +++ b/src/Numerics/SpecialFunctions/Evaluate.cs @@ -0,0 +1,247 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2012 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +// +// CERN - European Laboratory for Particle Physics +// http://www.docjar.com/html/api/cern/jet/math/Bessel.java.html +// Copyright 1999 CERN - European Laboratory for Particle Physics. +// Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose +// is hereby granted without fee, provided that the above copyright notice appear in all copies and +// that both that copyright notice and this permission notice appear in supporting documentation. +// CERN makes no representations about the suitability of this software for any purpose. +// It is provided "as is" without expressed or implied warranty. +// TOMS757 - Uncommon Special Functions (Fortran77) by Allan McLeod +// http://people.sc.fsu.edu/~jburkardt/f77_src/toms757/toms757.html +// Wei Wu +// Cephes Math Library, Stephen L. Moshier +// ALGLIB 2.0.1, Sergey Bochkanov +// + +// ReSharper disable CheckNamespace +namespace MathNet.Numerics +// ReSharper restore CheckNamespace +{ + using System; + + /// + /// Evaluation functions, useful for function approximation. + /// + public static class Evaluate + { + /// + /// Evaluate polynomials. + /// + /// The coefficients of the polynomial. + /// The location where to evaluate the polynomial at. + /// the evaluation of the polynomial. + public static double Polynomial(double[] coefficients, double z) + { + int count = coefficients.Length; + double sum = coefficients[count - 1]; + for (int i = count - 2; i >= 0; --i) + { + sum *= z; + sum += coefficients[i]; + } + + return sum; + } + + /// + /// Numerically stable series summation + /// + /// provides the summands sequentially + /// Sum + internal static double Series(Func nextSummand) + { + double compensation = 0.0; + double current; + const double factor = 1 << 16; + + double sum = nextSummand(); + + do + { + // Kahan Summation + // NOTE (ruegg): do NOT optimize. Now, how to tell that the compiler? + current = nextSummand(); + double y = current - compensation; + double t = sum + y; + compensation = t - sum; + compensation -= y; + sum = t; + } + while (Math.Abs(sum) < Math.Abs(factor * current)); + + return sum; + } + + /// Evaluates the series of Chebyshev polynomials Ti at argument x/2. + /// The series is given by + /// 
+        ///       N-1
+        ///        - '
+        /// y  =   >   coef[i] T (x/2)
+        ///        -            i
+        ///       i=0
+        ///
+ /// Coefficients are stored in reverse order, i.e. the zero + /// order term is last in the array. Note N is the number of + /// coefficients, not the order. + ///

+ /// If coefficients are for the interval a to b, x must + /// have been transformed to x -> 2(2x - b - a)/(b-a) before + /// entering the routine. This maps x from (a, b) to (-1, 1), + /// over which the Chebyshev polynomials are defined. + ///

+ /// If the coefficients are for the inverted interval, in + /// which (a, b) is mapped to (1/b, 1/a), the transformation + /// required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity, + /// this becomes x -> 4a/x - 1. + ///

+ /// SPEED: + ///

+ /// Taking advantage of the recurrence properties of the + /// Chebyshev polynomials, the routine requires one more + /// addition per loop than evaluating a nested polynomial of + /// the same degree. + ///

+ /// The coefficients of the polynomial. + /// + /// Argument to the polynomial. + /// + /// The number of coefficients. + /// + /// + /// Reference: https://bpm2.svn.codeplex.com/svn/Common.Numeric/Arithmetic.cs + ///

+ /// Marked as Deprecated in + /// http://people.apache.org/~isabel/mahout_site/mahout-matrix/apidocs/org/apache/mahout/jet/math/Arithmetic.html + /// + internal static double ChebyshevA(double[] coefficients, double x) + { + // TODO: Unify, normalize, then make public + + double b2; + + int p = 0; + + double b0 = coefficients[p++]; + double b1 = 0.0; + int i = coefficients.Length - 1; + + do + { + b2 = b1; + b1 = b0; + b0 = x * b1 - b2 + coefficients[p++]; + } + while (--i > 0); + + return (0.5 * (b0 - b2)); + } + + ///

+ /// Summation of Chebyshev polynomials, using the Clenshaw method with Reinsch modification. + /// + /// The no. of terms in the sequence. + /// The coefficients of the Chebyshev series, length n+1. + /// The value at which the series is to be evaluated. + /// + /// ORIGINAL AUTHOR: + /// Dr. Allan J. MacLeod; Dept. of Mathematics and Statistics, University of Paisley; High St., PAISLEY, SCOTLAND + /// REFERENCES: + /// "An error analysis of the modified Clenshaw method for evaluating Chebyshev and Fourier series" + /// J. Oliver, J.I.M.A., vol. 20, 1977, pp379-391 + /// + internal static double ChebyshevSum(int n, double[] coefficients, double x) + { + // TODO: Unify, normalize, then make public + + // If |x| < 0.6 use the standard Clenshaw method + if (Math.Abs(x) < 0.6) + { + double u0 = 0.0; + double u1 = 0.0; + double u2 = 0.0; + double xx = x + x; + + for (int i = n; i >= 0; i--) + { + u2 = u1; + u1 = u0; + u0 = xx * u1 + coefficients[i] - u2; + } + + return (u0 - u2) / 2.0; + } + + // If ABS ( T ) > = 0.6 use the Reinsch modification + // T > = 0.6 code + if (x > 0.0) + { + double u1 = 0.0; + double d1 = 0.0; + double d2 = 0.0; + double xx = (x - 0.5) - 0.5; + xx = xx + xx; + + for (int i = n; i >= 0; i--) + { + d2 = d1; + double u2 = u1; + d1 = xx * u2 + coefficients[i] + d2; + u1 = d1 + u2; + } + + return (d1 + d2) / 2.0; + } + else + { + // T < = -0.6 code + double u1 = 0.0; + double d1 = 0.0; + double d2 = 0.0; + double xx = (x + 0.5) + 0.5; + xx = xx + xx; + + for (int i = n; i >= 0; i--) + { + d2 = d1; + double u2 = u1; + d1 = xx * u2 + coefficients[i] - d2; + u1 = d1 - u2; + } + + return (d1 - d2) / 2.0; + } + } + } +} diff --git a/src/Numerics/SpecialFunctions/ModifiedBessel.cs b/src/Numerics/SpecialFunctions/ModifiedBessel.cs new file mode 100644 index 00000000..ec4e3716 --- /dev/null +++ b/src/Numerics/SpecialFunctions/ModifiedBessel.cs @@ -0,0 +1,287 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2012 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +// +// CERN - European Laboratory for Particle Physics +// http://www.docjar.com/html/api/cern/jet/math/Bessel.java.html +// Copyright 1999 CERN - European Laboratory for Particle Physics. +// Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose +// is hereby granted without fee, provided that the above copyright notice appear in all copies and +// that both that copyright notice and this permission notice appear in supporting documentation. +// CERN makes no representations about the suitability of this software for any purpose. +// It is provided "as is" without expressed or implied warranty. +// TOMS757 - Uncommon Special Functions (Fortran77) by Allan McLeod +// http://people.sc.fsu.edu/~jburkardt/f77_src/toms757/toms757.html +// Wei Wu +// Cephes Math Library, Stephen L. Moshier +// ALGLIB 2.0.1, Sergey Bochkanov +// + +// ReSharper disable CheckNamespace +namespace MathNet.Numerics +// ReSharper restore CheckNamespace +{ + using System; + + /// + /// This partial implementation of the SpecialFunctions class contains all methods related to the modified bessel function. + /// + public static partial class SpecialFunctions + { + /// + /// ************************************** + /// COEFFICIENTS FOR METHODS bessi0 * + /// ************************************** + /// + /// Chebyshev coefficients for exp(-x) I0(x) + /// in the interval [0, 8]. + /// + /// lim(x->0){ exp(-x) I0(x) } = 1. + /// + private static readonly double[] BesselI0A = new[] { -4.41534164647933937950e-18, 3.33079451882223809783e-17, -2.43127984654795469359e-16, 1.71539128555513303061e-15, -1.16853328779934516808e-14, 7.67618549860493561688e-14, -4.85644678311192946090e-13, 2.95505266312963983461e-12, -1.72682629144155570723e-11, 9.67580903537323691224e-11, -5.18979560163526290666e-10, 2.65982372468238665035e-9, -1.30002500998624804212e-8, 6.04699502254191894932e-8, -2.67079385394061173391e-7, 1.11738753912010371815e-6, -4.41673835845875056359e-6, 1.64484480707288970893e-5, -5.75419501008210370398e-5, 1.88502885095841655729e-4, -5.76375574538582365885e-4, 1.63947561694133579842e-3, -4.32430999505057594430e-3, 1.05464603945949983183e-2, -2.37374148058994688156e-2, 4.93052842396707084878e-2, -9.49010970480476444210e-2, 1.71620901522208775349e-1, -3.04682672343198398683e-1, 6.76795274409476084995e-1 }; + + /// Chebyshev coefficients for exp(-x) sqrt(x) I0(x) + /// in the inverted interval [8, infinity]. + /// + /// lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi). + /// + private static readonly double[] BesselI0B = new[] { -7.23318048787475395456e-18, -4.83050448594418207126e-18, 4.46562142029675999901e-17, 3.46122286769746109310e-17, -2.82762398051658348494e-16, -3.42548561967721913462e-16, 1.77256013305652638360e-15, 3.81168066935262242075e-15, -9.55484669882830764870e-15, -4.15056934728722208663e-14, 1.54008621752140982691e-14, 3.85277838274214270114e-13, 7.18012445138366623367e-13, -1.79417853150680611778e-12, -1.32158118404477131188e-11, -3.14991652796324136454e-11, 1.18891471078464383424e-11, 4.94060238822496958910e-10, 3.39623202570838634515e-9, 2.26666899049817806459e-8, 2.04891858946906374183e-7, 2.89137052083475648297e-6, 6.88975834691682398426e-5, 3.36911647825569408990e-3, 8.04490411014108831608e-1 }; + + /// + /// ************************************** + /// COEFFICIENTS FOR METHODS bessi1 * + /// ************************************** + /// + /// Chebyshev coefficients for exp(-x) I1(x) / x + /// in the interval [0, 8]. + /// + /// lim(x->0){ exp(-x) I1(x) / x } = 1/2. + /// + private static readonly double[] BesselI1A = new[] { 2.77791411276104639959e-18, -2.11142121435816608115e-17, 1.55363195773620046921e-16, -1.10559694773538630805e-15, 7.60068429473540693410e-15, -5.04218550472791168711e-14, 3.22379336594557470981e-13, -1.98397439776494371520e-12, 1.17361862988909016308e-11, -6.66348972350202774223e-11, 3.62559028155211703701e-10, -1.88724975172282928790e-9, 9.38153738649577178388e-9, -4.44505912879632808065e-8, 2.00329475355213526229e-7, -8.56872026469545474066e-7, 3.47025130813767847674e-6, -1.32731636560394358279e-5, 4.78156510755005422638e-5, -1.61760815825896745588e-4, 5.12285956168575772895e-4, -1.51357245063125314899e-3, 4.15642294431288815669e-3, -1.05640848946261981558e-2, 2.47264490306265168283e-2, -5.29459812080949914269e-2, 1.02643658689847095384e-1, -1.76416518357834055153e-1, 2.52587186443633654823e-1 }; + + /// Chebyshev coefficients for exp(-x) sqrt(x) I1(x) + /// in the inverted interval [8, infinity]. + /// + /// lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi). + /// + private static readonly double[] BesselI1B = new[] { 7.51729631084210481353e-18, 4.41434832307170791151e-18, -4.65030536848935832153e-17, -3.20952592199342395980e-17, 2.96262899764595013876e-16, 3.30820231092092828324e-16, -1.88035477551078244854e-15, -3.81440307243700780478e-15, 1.04202769841288027642e-14, 4.27244001671195135429e-14, -2.10154184277266431302e-14, -4.08355111109219731823e-13, -7.19855177624590851209e-13, 2.03562854414708950722e-12, 1.41258074366137813316e-11, 3.25260358301548823856e-11, -1.89749581235054123450e-11, -5.58974346219658380687e-10, -3.83538038596423702205e-9, -2.63146884688951950684e-8, -2.51223623787020892529e-7, -3.88256480887769039346e-6, -1.10588938762623716291e-4, -9.76109749136146840777e-3, 7.78576235018280120474e-1 }; + + /// + /// ************************************** + /// COEFFICIENTS FOR METHODS bessk0, bessk0e * + /// ************************************** + /// + /// Chebyshev coefficients for K0(x) + log(x/2) I0(x) + /// in the interval [0, 2]. The odd order coefficients are all + /// zero; only the even order coefficients are listed. + /// + /// lim(x->0){ K0(x) + log(x/2) I0(x) } = -EUL. + /// + private static readonly double[] BesselK0A = new[] { 1.37446543561352307156e-16, 4.25981614279661018399e-14, 1.03496952576338420167e-11, 1.90451637722020886025e-9, 2.53479107902614945675e-7, 2.28621210311945178607e-5, 1.26461541144692592338e-3, 3.59799365153615016266e-2, 3.44289899924628486886e-1, -5.35327393233902768720e-1 }; + + /// Chebyshev coefficients for exp(x) sqrt(x) K0(x) + /// in the inverted interval [2, infinity]. + /// + /// lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2). + /// + private static readonly double[] BesselK0B = new[] { 5.30043377268626276149e-18, -1.64758043015242134646e-17, 5.21039150503902756861e-17, -1.67823109680541210385e-16, 5.51205597852431940784e-16, -1.84859337734377901440e-15, 6.34007647740507060557e-15, -2.22751332699166985548e-14, 8.03289077536357521100e-14, -2.98009692317273043925e-13, 1.14034058820847496303e-12, -4.51459788337394416547e-12, 1.85594911495471785253e-11, -7.95748924447710747776e-11, 3.57739728140030116597e-10, -1.69753450938905987466e-9, 8.57403401741422608519e-9, -4.66048989768794782956e-8, 2.76681363944501510342e-7, -1.83175552271911948767e-6, 1.39498137188764993662e-5, -1.28495495816278026384e-4, 1.56988388573005337491e-3, -3.14481013119645005427e-2, 2.44030308206595545468e0 }; + + /// + /// ************************************** + /// COEFFICIENTS FOR METHODS bessk1, bessk1e * + /// ************************************** + /// + /// Chebyshev coefficients for x(K1(x) - log(x/2) I1(x)) + /// in the interval [0, 2]. + /// + /// lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1. + /// + private static readonly double[] BesselK1A = new[] { -7.02386347938628759343e-18, -2.42744985051936593393e-15, -6.66690169419932900609e-13, -1.41148839263352776110e-10, -2.21338763073472585583e-8, -2.43340614156596823496e-6, -1.73028895751305206302e-4, -6.97572385963986435018e-3, -1.22611180822657148235e-1, -3.53155960776544875667e-1, 1.52530022733894777053e0 }; + + /// Chebyshev coefficients for exp(x) sqrt(x) K1(x) + /// in the interval [2, infinity]. + /// + /// lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2). + /// + private static readonly double[] BesselK1B = new[] { -5.75674448366501715755e-18, 1.79405087314755922667e-17, -5.68946255844285935196e-17, 1.83809354436663880070e-16, -6.05704724837331885336e-16, 2.03870316562433424052e-15, -7.01983709041831346144e-15, 2.47715442448130437068e-14, -8.97670518232499435011e-14, 3.34841966607842919884e-13, -1.28917396095102890680e-12, 5.13963967348173025100e-12, -2.12996783842756842877e-11, 9.21831518760500529508e-11, -4.19035475934189648750e-10, 2.01504975519703286596e-9, -1.03457624656780970260e-8, 5.74108412545004946722e-8, -3.50196060308781257119e-7, 2.40648494783721712015e-6, -1.93619797416608296024e-5, 1.95215518471351631108e-4, -2.85781685962277938680e-3, 1.03923736576817238437e-1, 2.72062619048444266945e0 }; + + /// Returns the modified Bessel function of first kind, order 0 of the argument. + ///

+ /// The function is defined as i0(x) = j0( ix ). + ///

+ /// The range is partitioned into the two intervals [0, 8] and + /// (8, infinity). Chebyshev polynomial expansions are employed + /// in each interval. + ///

+ /// The value to compute the bessel function of. + /// + public static double BesselI0(double x) + { + if (x < 0) + { + x = -x; + } + if (x <= 8.0) + { + double y = (x / 2.0) - 2.0; + return (Math.Exp(x) * Evaluate.ChebyshevA(BesselI0A, y)); + } + + double x1 = 32.0 / x - 2.0; + return (Math.Exp(x) * Evaluate.ChebyshevA(BesselI0B, x1) / Math.Sqrt(x)); + } + + /// Returns the modified Bessel function of first kind, + /// order 1 of the argument. + ///

+ /// The function is defined as i1(x) = -i j1( ix ). + ///

+ /// The range is partitioned into the two intervals [0, 8] and + /// (8, infinity). Chebyshev polynomial expansions are employed + /// in each interval. + ///

+ /// The value to compute the bessel function of. + /// + public static double BesselI1(double x) + { + double z = Math.Abs(x); + if (z <= 8.0) + { + double y = (z / 2.0) - 2.0; + z = Evaluate.ChebyshevA(BesselI1A, y) * z * Math.Exp(z); + } + else + { + double x1 = 32.0 / z - 2.0; + z = Math.Exp(z) * Evaluate.ChebyshevA(BesselI1B, x1) / Math.Sqrt(z); + } + if (x < 0.0) + { + z = -z; + } + return z; + } + + /// Returns the modified Bessel function of the second kind + /// of order 0 of the argument. + ///

+ /// The range is partitioned into the two intervals [0, 8] and + /// (8, infinity). Chebyshev polynomial expansions are employed + /// in each interval. + ///

+ /// The value to compute the bessel function of. + /// + public static double BesselK0(double x) + { + if (x <= 0.0) + { + throw new ArithmeticException(); + } + if (x <= 2.0) + { + double y = x * x - 2.0; + return Evaluate.ChebyshevA(BesselK0A, y) - Math.Log(0.5 * x) * BesselI0(x); + } + + double z = 8.0 / x - 2.0; + return Math.Exp(-x) * Evaluate.ChebyshevA(BesselK0B, z) / Math.Sqrt(x); + } + + /// Returns the exponentially scaled modified Bessel function + /// of the second kind of order 0 of the argument. + /// + /// The value to compute the bessel function of. + /// + public static double BesselK0e(double x) + { + if (x <= 0.0) + { + throw new ArithmeticException(); + } + if (x <= 2.0) + { + double y = x * x - 2.0; + return Evaluate.ChebyshevA(BesselK0A, y) - Math.Log(0.5 * x) * BesselI0(x) * Math.Exp(x); + } + + double x1 = 8.0 / x - 2.0; + return Evaluate.ChebyshevA(BesselK0B, x1) / Math.Sqrt(x); + } + + /// Returns the modified Bessel function of the second kind + /// of order 1 of the argument. + ///

+ /// The range is partitioned into the two intervals [0, 2] and + /// (2, infinity). Chebyshev polynomial expansions are employed + /// in each interval. + ///

+ /// The value to compute the bessel function of. + /// + public static double BesselK1(double x) + { + double z = 0.5 * x; + if (z <= 0.0) + { + throw new ArithmeticException(); + } + if (x <= 2.0) + { + double y = x * x - 2.0; + return Math.Log(z) * BesselI1(x) + Evaluate.ChebyshevA(BesselK1A, y) / x; + } + + double x1 = 8.0 / x - 2.0; + return Math.Exp(-x) * Evaluate.ChebyshevA(BesselK1B, x1) / Math.Sqrt(x); + } + + /// Returns the exponentially scaled modified Bessel function + /// of the second kind of order 1 of the argument. + ///

+ /// k1e(x) = exp(x) * k1(x). + ///

+ /// The value to compute the bessel function of. + /// + public static double BesselK1e(double x) + { + if (x <= 0.0) + { + throw new ArithmeticException(); + } + if (x <= 2.0) + { + double y = x * x - 2.0; + return Math.Log(0.5 * x) * BesselI1(x) + Evaluate.ChebyshevA(BesselK1A, y) / x * Math.Exp(x); + } + + double x1 = 8.0 / x - 2.0; + return Evaluate.ChebyshevA(BesselK1B, x1) / Math.Sqrt(x); + } + } +} diff --git a/src/Numerics/SpecialFunctions/ModifiedStruve.cs b/src/Numerics/SpecialFunctions/ModifiedStruve.cs new file mode 100644 index 00000000..504c5357 --- /dev/null +++ b/src/Numerics/SpecialFunctions/ModifiedStruve.cs @@ -0,0 +1,556 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2012 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +// +// CERN - European Laboratory for Particle Physics +// http://www.docjar.com/html/api/cern/jet/math/Bessel.java.html +// Copyright 1999 CERN - European Laboratory for Particle Physics. +// Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose +// is hereby granted without fee, provided that the above copyright notice appear in all copies and +// that both that copyright notice and this permission notice appear in supporting documentation. +// CERN makes no representations about the suitability of this software for any purpose. +// It is provided "as is" without expressed or implied warranty. +// TOMS757 - Uncommon Special Functions (Fortran77) by Allan McLeod +// http://people.sc.fsu.edu/~jburkardt/f77_src/toms757/toms757.html +// Wei Wu +// Cephes Math Library, Stephen L. Moshier +// ALGLIB 2.0.1, Sergey Bochkanov +// + +// ReSharper disable CheckNamespace +namespace MathNet.Numerics +// ReSharper restore CheckNamespace +{ + using System; + + /// + /// This partial implementation of the SpecialFunctions class contains all methods related to the modified bessel function. + /// + public static partial class SpecialFunctions + { + /// + /// Returns the modified Struve function of order 0. + /// + /// The value to compute the function of. + /// + public static double StruveL0(double x) + { + //*********************************************************************72 + // + //c STRVL0 calculates the modified Struve function of order 0. + // + // DESCRIPTION: + // + // This function calculates the modified Struve function of + // order 0, denoted L0(x), defined as the solution of the + // second-order equation + // + // x*D(Df) + Df - x*f = 2x/pi + // + // This subroutine is set up to work on IEEE machines. + // For other machines, you should retrieve the code + // from the general MISCFUN archive. + // + // + // ERROR RETURNS: + // + // If the value of |XVALUE| is too large, the result + // would cause an floating-pt overflow. An error message + // is printed and the function returns the value of + // sign(XVALUE)*XMAX where XMAX is the largest possible + // floating-pt argument. + // + // + // MACHINE-DEPENDENT PARAMETERS: + // + // NTERM1 - INTEGER - The no. of terms for the array ARL0. + // The recommended value is such that + // ABS(ARL0(NTERM1)) < EPS/100 + // + // NTERM2 - INTEGER - The no. of terms for the array ARL0AS. + // The recommended value is such that + // ABS(ARL0AS(NTERM2)) < EPS/100 + // + // NTERM3 - INTEGER - The no. of terms for the array AI0ML0. + // The recommended value is such that + // ABS(AI0ML0(NTERM3)) < EPS/100 + // + // XLOW - DOUBLE PRECISION - The value of x below which L0(x) = 2*x/pi + // to machine precision. The recommended value is + // 3*SQRT(EPS) + // + // XHIGH1 - DOUBLE PRECISION - The value beyond which the Chebyshev series + // in the asymptotic expansion of I0 - L0 gives + // 1.0 to machine precision. The recommended value + // is SQRT( 30/EPSNEG ) + // + // XHIGH2 - DOUBLE PRECISION - The value beyond which the Chebyshev series + // in the asymptotic expansion of I0 gives 1.0 + // to machine precision. The recommended value + // is 28 / EPSNEG + // + // XMAX - DOUBLE PRECISION - The value of XMAX, where XMAX is the + // largest possible floating-pt argument. + // This is used to prevent overflow. + // + // For values of EPS, EPSNEG and XMAX the user should refer + // to the file MACHCON.TXT + // + // The machine-arithmetic constants are given in DATA + // statements. + // + // + // INTRINSIC FUNCTIONS USED: + // + // EXP , LOG , SQRT + // + // + // OTHER MISCFUN SUBROUTINES USED: + // + // CHEVAL , ERRPRN + // + // + // AUTHOR: + // DR. ALLAN J. MACLEOD + // DEPT. OF MATHEMATICS AND STATISTICS + // UNIVERSITY OF PAISLEY + // HIGH ST. + // PAISLEY + // SCOTLAND + // PA1 2BE + // + // (e-mail: macl_ms0@paisley.ac.uk ) + // + // + // LATEST REVISION: + // 12 JANUARY, 1996 + // + // + + if (x < 0.0) + { + return -StruveL0(-x); + } + + const double LNR2PI = 0.91893853320467274178; + const double TWOBPI = 0.63661977236758134308; + + double[] ARL0 = new double[28]; + ARL0[0] = 0.42127458349979924863; + ARL0[1] = -0.33859536391220612188; + ARL0[2] = 0.21898994812710716064; + ARL0[3] = -0.12349482820713185712; + ARL0[4] = 0.6214209793866958440e-1; + ARL0[5] = -0.2817806028109547545e-1; + ARL0[6] = 0.1157419676638091209e-1; + ARL0[7] = -0.431658574306921179e-2; + ARL0[8] = 0.146142349907298329e-2; + ARL0[9] = -0.44794211805461478e-3; + ARL0[10] = 0.12364746105943761e-3; + ARL0[11] = -0.3049028334797044e-4; + ARL0[12] = 0.663941401521146e-5; + ARL0[13] = -0.125538357703889e-5; + ARL0[14] = 0.20073446451228e-6; + ARL0[15] = -0.2588260170637e-7; + ARL0[16] = 0.241143742758e-8; + ARL0[17] = -0.10159674352e-9; + ARL0[18] = -0.1202430736e-10; + ARL0[19] = 0.262906137e-11; + ARL0[20] = -0.15313190e-12; + ARL0[21] = -0.1574760e-13; + ARL0[22] = 0.315635e-14; + ARL0[23] = -0.4096e-16; + ARL0[24] = -0.3620e-16; + ARL0[25] = 0.239e-17; + ARL0[26] = 0.36e-18; + ARL0[27] = -0.4e-19; + + double[] ARL0AS = new double[16]; + ARL0AS[0] = 2.00861308235605888600; + ARL0AS[1] = 0.403737966500438470e-2; + ARL0AS[2] = -0.25199480286580267e-3; + ARL0AS[3] = 0.1605736682811176e-4; + ARL0AS[4] = -0.103692182473444e-5; + ARL0AS[5] = 0.6765578876305e-7; + ARL0AS[6] = -0.444999906756e-8; + ARL0AS[7] = 0.29468889228e-9; + ARL0AS[8] = -0.1962180522e-10; + ARL0AS[9] = 0.131330306e-11; + ARL0AS[10] = -0.8819190e-13; + ARL0AS[11] = 0.595376e-14; + ARL0AS[12] = -0.40389e-15; + ARL0AS[13] = 0.2651e-16; + ARL0AS[14] = -0.208e-17; + ARL0AS[15] = 0.11e-18; + + double[] AI0ML0 = new double[24]; + AI0ML0[0] = 2.00326510241160643125; + AI0ML0[1] = 0.195206851576492081e-2; + AI0ML0[2] = 0.38239523569908328e-3; + AI0ML0[3] = 0.7534280817054436e-4; + AI0ML0[4] = 0.1495957655897078e-4; + AI0ML0[5] = 0.299940531210557e-5; + AI0ML0[6] = 0.60769604822459e-6; + AI0ML0[7] = 0.12399495544506e-6; + AI0ML0[8] = 0.2523262552649e-7; + AI0ML0[9] = 0.504634857332e-8; + AI0ML0[10] = 0.97913236230e-9; + AI0ML0[11] = 0.18389115241e-9; + AI0ML0[12] = 0.3376309278e-10; + AI0ML0[13] = 0.611179703e-11; + AI0ML0[14] = 0.108472972e-11; + AI0ML0[15] = 0.18861271e-12; + AI0ML0[16] = 0.3280345e-13; + AI0ML0[17] = 0.565647e-14; + AI0ML0[18] = 0.93300e-15; + AI0ML0[19] = 0.15881e-15; + AI0ML0[20] = 0.2791e-16; + AI0ML0[21] = 0.389e-17; + AI0ML0[22] = 0.70e-18; + AI0ML0[23] = 0.16e-18; + + // MACHINE-DEPENDENT VALUES (Suitable for IEEE-arithmetic machines) + const int NTERM1 = 25; const int NTERM2 = 14; const int NTERM3 = 21; + const double XLOW = 4.4703484e-8; const double XMAX = 1.797693e308; + const double XHIGH1 = 5.1982303e8; const double XHIGH2 = 2.5220158e17; + + // Code for |xvalue| <= 16 + if (x <= 16.0) + { + if (x < XLOW) + { + return TWOBPI * x; + } + + double T = (4.0 * x - 24.0) / (x + 24.0); + return TWOBPI * x * Evaluate.ChebyshevSum(NTERM1, ARL0, T) * Math.Exp(x); + } + + // Code for |xvalue| > 16 + double ch1; + if (x > XHIGH2) + { + ch1 = 1.0; + } + else + { + double T = (x - 28.0) / (4.0 - x); + ch1 = Evaluate.ChebyshevSum(NTERM2, ARL0AS, T); + } + + double ch2; + if (x > XHIGH1) + { + ch2 = 1.0; + } + else + { + double xsq = x * x; + double T = (800.0 - xsq) / (288.0 + xsq); + ch2 = Evaluate.ChebyshevSum(NTERM3, AI0ML0, T); + } + + double test = Math.Log(ch1) - LNR2PI - Math.Log(x) / 2.0 + x; + if (test > Math.Log(XMAX)) + { + throw new ArithmeticException("ERROR IN MISCFUN FUNCTION STRVL0: ARGUMENT CAUSES OVERFLOW"); + } + + return Math.Exp(test) - TWOBPI * ch2 / x; + } + + /// + /// Returns the modified Struve function of order 1. + /// + /// The value to compute the function of. + /// + public static double StruveL1(double x) + { + //*********************************************************************72 + // + //c STRVL1 calculates the modified Struve function of order 1. + // + // DESCRIPTION: + // + // This function calculates the modified Struve function of + // order 1, denoted L1(x), defined as the solution of + // + // x*x*D(Df) + x*Df - (x*x+1)f = 2*x*x/pi + // + // This subroutine is set up to work on IEEE machines. + // For other machines, you should retrieve the code + // from the general MISCFUN archive. + // + // + // ERROR RETURNS: + // + // If the value of |XVALUE| is too large, the result + // would cause an floating-pt overflow. An error message + // is printed and the function returns the value of + // sign(XVALUE)*XMAX where XMAX is the largest possible + // floating-pt argument. + // + // + // MACHINE-DEPENDENT PARAMETERS: + // + // NTERM1 - INTEGER - The no. of terms for the array ARL1. + // The recommended value is such that + // ABS(ARL1(NTERM1)) < EPS/100 + // + // NTERM2 - INTEGER - The no. of terms for the array ARL1AS. + // The recommended value is such that + // ABS(ARL1AS(NTERM2)) < EPS/100 + // + // NTERM3 - INTEGER - The no. of terms for the array AI1ML1. + // The recommended value is such that + // ABS(AI1ML1(NTERM3)) < EPS/100 + // + // XLOW1 - DOUBLE PRECISION - The value of x below which + // L1(x) = 2*x*x/(3*pi) + // to machine precision. The recommended + // value is SQRT(15*EPS) + // + // XLOW2 - DOUBLE PRECISION - The value of x below which L1(x) set to 0.0. + // This is used to prevent underflow. The + // recommended value is + // SQRT(5*XMIN) + // + // XHIGH1 - DOUBLE PRECISION - The value of |x| above which the Chebyshev + // series in the asymptotic expansion of I1 + // equals 1.0 to machine precision. The + // recommended value is SQRT( 30 / EPSNEG ). + // + // XHIGH2 - DOUBLE PRECISION - The value of |x| above which the Chebyshev + // series in the asymptotic expansion of I1 - L1 + // equals 1.0 to machine precision. The recommended + // value is 30 / EPSNEG. + // + // XMAX - DOUBLE PRECISION - The value of XMAX, where XMAX is the + // largest possible floating-pt argument. + // This is used to prevent overflow. + // + // For values of EPS, EPSNEG, XMIN, and XMAX the user should refer + // to the file MACHCON.TXT + // + // The machine-arithmetic constants are given in DATA + // statements. + // + // + // INTRINSIC FUNCTIONS USED: + // + // EXP , LOG , SQRT + // + // + // OTHER MISCFUN SUBROUTINES USED: + // + // CHEVAL , ERRPRN + // + // + // AUTHOR: + // DR. ALLAN J. MACLEOD + // DEPT. OF MATHEMATICS AND STATISTICS + // UNIVERSITY OF PAISLEY + // HIGH ST. + // PAISLEY + // SCOTLAND + // PA1 2BE + // + // (e-mail: macl_ms0@paisley.ac.uk ) + // + // + // LATEST UPDATE: + // 12 JANUARY, 1996 + // + // + + if (x < 0.0) + { + return StruveL1(-x); + } + + const double LNR2PI = 0.91893853320467274178; + const double PI3BY2 = 4.71238898038468985769; + const double TWOBPI = 0.63661977236758134308; + + double[] ARL1 = new double[27]; + ARL1[0] = 0.38996027351229538208; + ARL1[1] = -0.33658096101975749366; + ARL1[2] = 0.23012467912501645616; + ARL1[3] = -0.13121594007960832327; + ARL1[4] = 0.6425922289912846518e-1; + ARL1[5] = -0.2750032950616635833e-1; + ARL1[6] = 0.1040234148637208871e-1; + ARL1[7] = -0.350532294936388080e-2; + ARL1[8] = 0.105748498421439717e-2; + ARL1[9] = -0.28609426403666558e-3; + ARL1[10] = 0.6925708785942208e-4; + ARL1[11] = -0.1489693951122717e-4; + ARL1[12] = 0.281035582597128e-5; + ARL1[13] = -0.45503879297776e-6; + ARL1[14] = 0.6090171561770e-7; + ARL1[15] = -0.623543724808e-8; + ARL1[16] = 0.38430012067e-9; + ARL1[17] = 0.790543916e-11; + ARL1[18] = -0.489824083e-11; + ARL1[19] = 0.46356884e-12; + ARL1[20] = 0.684205e-14; + ARL1[21] = -0.569748e-14; + ARL1[22] = 0.35324e-15; + ARL1[23] = 0.4244e-16; + ARL1[24] = -0.644e-17; + ARL1[25] = -0.21e-18; + ARL1[26] = 0.9e-19; + + double[] ARL1AS = new double[17]; + ARL1AS[0] = 1.97540378441652356868; + ARL1AS[1] = -0.1195130555088294181e-1; + ARL1AS[2] = 0.33639485269196046e-3; + ARL1AS[3] = -0.1009115655481549e-4; + ARL1AS[4] = 0.30638951321998e-6; + ARL1AS[5] = -0.953704370396e-8; + ARL1AS[6] = 0.29524735558e-9; + ARL1AS[7] = -0.951078318e-11; + ARL1AS[8] = 0.28203667e-12; + ARL1AS[9] = -0.1134175e-13; + ARL1AS[10] = 0.147e-17; + ARL1AS[11] = -0.6232e-16; + ARL1AS[12] = -0.751e-17; + ARL1AS[13] = -0.17e-18; + ARL1AS[14] = 0.51e-18; + ARL1AS[15] = 0.23e-18; + ARL1AS[16] = 0.5e-19; + + double[] AI1ML1 = new double[26]; + AI1ML1[0] = 1.99679361896789136501; + AI1ML1[1] = -0.190663261409686132e-2; + AI1ML1[2] = -0.36094622410174481e-3; + AI1ML1[3] = -0.6841847304599820e-4; + AI1ML1[4] = -0.1299008228509426e-4; + AI1ML1[5] = -0.247152188705765e-5; + AI1ML1[6] = -0.47147839691972e-6; + AI1ML1[7] = -0.9020819982592e-7; + AI1ML1[8] = -0.1730458637504e-7; + AI1ML1[9] = -0.332323670159e-8; + AI1ML1[10] = -0.63736421735e-9; + AI1ML1[11] = -0.12180239756e-9; + AI1ML1[12] = -0.2317346832e-10; + AI1ML1[13] = -0.439068833e-11; + AI1ML1[14] = -0.82847110e-12; + AI1ML1[15] = -0.15562249e-12; + AI1ML1[16] = -0.2913112e-13; + AI1ML1[17] = -0.543965e-14; + AI1ML1[18] = -0.101177e-14; + AI1ML1[19] = -0.18767e-15; + AI1ML1[20] = -0.3484e-16; + AI1ML1[21] = -0.643e-17; + AI1ML1[22] = -0.118e-17; + AI1ML1[23] = -0.22e-18; + AI1ML1[24] = -0.4e-19; + AI1ML1[25] = -0.1e-19; + + // MACHINE-DEPENDENT VALUES (Suitable for IEEE-arithmetic machines) + const int NTERM1 = 24; const int NTERM2 = 13; const int NTERM3 = 22; + const double XLOW1 = 5.7711949e-8; const double XLOW2 = 3.3354714e-154; const double XMAX = 1.797693e308; + const double XHIGH1 = 5.19823025e8; const double XHIGH2 = 2.7021597e17; + + // CODE FOR |x| <= 16 + if (x <= 16.0) + { + if (x <= XLOW2) + { + return 0.0; + } + + double xsq = x * x; + if (x < XLOW1) + { + return xsq / PI3BY2; + } + + double t = (4.0 * x - 24.0) / (x + 24.0); + return xsq * Evaluate.ChebyshevSum(NTERM1, ARL1, t) * Math.Exp(x) / PI3BY2; + } + + // CODE FOR |x| > 16 + double ch1; + if (x > XHIGH2) + { + ch1 = 1.0; + } + else + { + double t = (x - 30.0) / (2.0 - x); + ch1 = Evaluate.ChebyshevSum(NTERM2, ARL1AS, t); + } + + double ch2; + if (x > XHIGH1) + { + ch2 = 1.0; + } + else + { + double xsq = x * x; + double t = (800.0 - xsq) / (288.0 + xsq); + ch2 = Evaluate.ChebyshevSum(NTERM3, AI1ML1, t); + } + + double test = Math.Log(ch1) - LNR2PI - Math.Log(x) / 2.0 + x; + if (test > Math.Log(XMAX)) + { + throw new ArithmeticException("ERROR IN MISCFUN FUNCTION STRVL1: ARGUMENT CAUSES OVERFLOW"); + } + + return Math.Exp(test) - TWOBPI * ch2; + } + + /// + /// Returns the difference between the Bessel I0 and Struve L0 functions. + /// + /// The value to compute the function of. + /// + public static double BesselI0MStruveL0(double x) + { + // TODO: way off for large x (e.g. 100) - needs direct approximation + return BesselI0(x) - StruveL0(x); + } + + /// + /// Returns the difference between the Bessel I1 and Struve L1 functions. + /// + /// The value to compute the function of. + /// + public static double BesselI1MStruveL1(double x) + { + // TODO: way off for large x (e.g. 100) - needs direct approximation + return BesselI1(x) - StruveL1(x); + } + } +} diff --git a/src/Numerics/SpecialFunctions/Stability.cs b/src/Numerics/SpecialFunctions/Stability.cs index 7cf72d85..8955091a 100644 --- a/src/Numerics/SpecialFunctions/Stability.cs +++ b/src/Numerics/SpecialFunctions/Stability.cs @@ -58,7 +58,7 @@ namespace MathNet.Numerics // Series Expansion to x^k / k! int k = 0; double term = 1.0; - return Series( + return Evaluate.Series( () => { k++; @@ -164,34 +164,5 @@ namespace MathNet.Numerics return 0f; } - - /// - /// Numerically stable series summation - /// - /// provides the summands sequentially - /// Sum - private static double Series(Func nextSummand) - { - double compensation = 0.0; - double current; - const double factor = 1 << 16; - - double sum = nextSummand(); - - do - { - // Kahan Summation - // NOTE (ruegg): do NOT optimize. Now, how to tell that the compiler? - current = nextSummand(); - double y = current - compensation; - double t = sum + y; - compensation = t - sum; - compensation -= y; - sum = t; - } - while (Math.Abs(sum) < Math.Abs(factor * current)); - - return sum; - } } } diff --git a/src/Portable/Portable.csproj b/src/Portable/Portable.csproj index f44e3c78..f6c13134 100644 --- a/src/Portable/Portable.csproj +++ b/src/Portable/Portable.csproj @@ -990,6 +990,9 @@ SpecialFunctions\Erf.cs + + SpecialFunctions\Evaluate.cs + SpecialFunctions\Factorial.cs @@ -1002,6 +1005,12 @@ SpecialFunctions\Logistic.cs + + SpecialFunctions\ModifiedBessel.cs + + + SpecialFunctions\ModifiedStruve.cs + SpecialFunctions\Stability.cs diff --git a/src/UnitTests/SpecialFunctionsTests/ErfTests.cs b/src/UnitTests/SpecialFunctionsTests/ErfTests.cs index bca2f68d..e0a385e7 100644 --- a/src/UnitTests/SpecialFunctionsTests/ErfTests.cs +++ b/src/UnitTests/SpecialFunctionsTests/ErfTests.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com +// // Copyright (c) 2009-2010 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND diff --git a/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs b/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs index 903ddef0..7bd777b8 100644 --- a/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs +++ b/src/UnitTests/SpecialFunctionsTests/FactorialTest.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com +// // Copyright (c) 2009-2010 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND diff --git a/src/UnitTests/SpecialFunctionsTests/GammaTests.cs b/src/UnitTests/SpecialFunctionsTests/GammaTests.cs index 690e8e83..94ea8681 100644 --- a/src/UnitTests/SpecialFunctionsTests/GammaTests.cs +++ b/src/UnitTests/SpecialFunctionsTests/GammaTests.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com +// // Copyright (c) 2009-2010 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -32,6 +36,7 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests /// /// Gamma functions tests. /// + [TestFixture] public class GammaTests { /// diff --git a/src/UnitTests/SpecialFunctionsTests/ModifiedBesselTests.cs b/src/UnitTests/SpecialFunctionsTests/ModifiedBesselTests.cs new file mode 100644 index 00000000..4bbedede --- /dev/null +++ b/src/UnitTests/SpecialFunctionsTests/ModifiedBesselTests.cs @@ -0,0 +1,120 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2012 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests +{ + using System; + using NUnit.Framework; + + /// + /// Modified Bessel functions tests. + /// + [TestFixture] + public class BodifiedBesselTests + { + [Test] + public void BesselI0Approx([Range(-3.75, 3.75, 0.25)] double x) + { + // Approx by Abramowitz/Stegun 9.8.1 + Assert.AreEqual(Evaluate.Polynomial(new[] { 1.0, 0.0, 3.5156229, 0.0, 3.0899424, 0.0, 1.2067492, 0.0, 0.2659732, 0.0, 0.0360768, 0.0, 0.0045813 }, x / 3.75), SpecialFunctions.BesselI0(x), 1e-7); + } + + [TestCase(0.0, 1.0)] + [TestCase(0.005, 1.000006250009766)] + [TestCase(0.5, 1.063483370741324)] + [TestCase(1.5, 1.646723189772891)] + [TestCase(10.0, 2815.716628466254)] + [TestCase(100.0, 1.073751707131074e+42)] + [TestCase(-0.005, 1.000006250009766)] + [TestCase(-10.0, 2815.716628466254)] + public void BesselI0Exact(double x, double expected) + { + AssertHelpers.AlmostEqual(expected, SpecialFunctions.BesselI0(x), 14); + } + + [Test] + public void BesselI1Approx([Range(-3.75, 3.75, 0.25)] double x) + { + // Approx by Abramowitz/Stegun 9.8.3 + Assert.AreEqual(Evaluate.Polynomial(new[] { 0.5, 0.0, 0.87890594, 0.0, 0.51498869, 0.0, 0.15084934, 0.0, 0.02658733, 0.0, 0.00301532, 0.0, 0.00032411 }, x / 3.75) * x, SpecialFunctions.BesselI1(x), 1e-8); + } + + [TestCase(0.0, 0.0)] + [TestCase(0.005, 0.002500007812508138)] + [TestCase(0.5, 0.2578943053908963)] + [TestCase(1.5, 0.9816664285779076)] + [TestCase(10.0, 2670.988303701255)] + [TestCase(100.0, 1.068369390338162e+42)] + [TestCase(-0.005, -0.002500007812508138)] + [TestCase(-10.0, -2670.988303701255)] + public void BesselI1Exact(double x, double expected) + { + AssertHelpers.AlmostEqual(expected, SpecialFunctions.BesselI1(x), 14); + } + + [Test] + public void BesselK0Approx([Range(0.20, 2.0, 0.20)] double x) + { + // Approx by Abramowitz/Stegun 9.8.5 + Assert.AreEqual(Evaluate.Polynomial(new[] { -Math.Log(x/2.0)*SpecialFunctions.BesselI0(x)-0.57721566, 0.0, 0.42278420, 0.0, 0.23069756, 0.0, 0.03488590, 0.0, 0.00262698, 0.0, 0.00010750, 0.0, 0.00000740 }, x / 2.0), SpecialFunctions.BesselK0(x), 1e-8); + } + + [TestCase(1e-10, 23.14178244559887)] + [TestCase(1e-5, 11.62885698094436)] + [TestCase(0.005, 5.414288971329485)] + [TestCase(0.5, 0.9244190712276659)] + [TestCase(1.5, 0.2138055626475257)] + [TestCase(10.0, 0.00001778006231616765)] + [TestCase(100.0, 4.656628229175902e-45)] + public void BesselK0Exact(double x, double expected) + { + AssertHelpers.AlmostEqual(expected, SpecialFunctions.BesselK0(x), 14); + } + + [Test] + public void BesselK1Approx([Range(0.20, 2.0, 0.20)] double x) + { + // Approx by Abramowitz/Stegun 9.8.7 + Assert.AreEqual(Evaluate.Polynomial(new[] { x * Math.Log(x / 2.0) * SpecialFunctions.BesselI1(x) + 1.0, 0.0, 0.15443144, 0.0, -0.67278579, 0.0, -0.18156897, 0.0, -0.01919402, 0.0, -0.00110404, 0.0, -0.00004686 }, x / 2.0), SpecialFunctions.BesselK1(x) * x, 1e-8); + } + + [TestCase(1e-10, 1.0e+10)] + [TestCase(1e-5, 99999.99993935572)] + [TestCase(0.005, 199.9852143257300)] + [TestCase(0.5, 1.656441120003301)] + [TestCase(1.5, 0.2773878004568438)] + [TestCase(10.0, 0.00001864877345382558)] + [TestCase(100.0, 4.679853735636909e-45)] + public void BesselK1Exact(double x, double expected) + { + AssertHelpers.AlmostEqual(expected, SpecialFunctions.BesselK1(x), 14); + } + } +} diff --git a/src/UnitTests/SpecialFunctionsTests/ModifiedStruveTests.cs b/src/UnitTests/SpecialFunctionsTests/ModifiedStruveTests.cs new file mode 100644 index 00000000..047a14de --- /dev/null +++ b/src/UnitTests/SpecialFunctionsTests/ModifiedStruveTests.cs @@ -0,0 +1,99 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2012 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests +{ + using NUnit.Framework; + + /// + /// Modified Struve functions tests. + /// + [TestFixture] + public class ModifiedStruveTests + { + [TestCase(0.0, 0.0)] + [TestCase(0.005, 0.003183107703788032)] + [TestCase(0.5, 0.3272406993941808)] + [TestCase(1.5, 1.216162510717182)] + [TestCase(10.0, 2815.652249374595)] + [TestCase(100.0, 1.073751707131074e+42)] + [TestCase(-0.005, -0.003183107703788032)] + [TestCase(-10.0, -2815.652249374595)] + public void StruveL0Exact(double x, double expected) + { + AssertHelpers.AlmostEqual(expected, SpecialFunctions.StruveL0(x), 14); + } + + [TestCase(0.0, 0.0)] + [TestCase(0.005, 5.305173611677443e-6)] + [TestCase(0.5, 0.05394218262352266)] + [TestCase(1.5, 0.5538569084469910)] + [TestCase(10.0, 2670.358285208483)] + [TestCase(100.0, 1.068369390338162e+42)] + [TestCase(-0.005, 5.305173611677443e-6)] + [TestCase(-10.0, 2670.358285208483)] + public void StruveL1Exact(double x, double expected) + { + AssertHelpers.AlmostEqual(expected, SpecialFunctions.StruveL1(x), 14); + } + + [TestCase(0.0, 1.0)] + [TestCase(0.1, 0.938769)] + [TestCase(0.4, 0.781198)] + [TestCase(2.0, 0.342152)] + [TestCase(2.4, 0.289765)] + [TestCase(4.5, 0.150279)] + [TestCase(4.9, 0.136938)] + [TestCase(10.0, 0.064379)] + [TestCase(20.0, 0.031912)] + //[TestCase(100.0, 0.006367)] Needs direct approximation + public void BesselI0MStruveL0Exact(double x, double expected) + { + // Abramowitz/Stegun Table 12.1, 12.2 + AssertHelpers.AlmostEqual(expected, SpecialFunctions.BesselI0MStruveL0(x), 5); + } + + [TestCase(0.0, 0.0)] + [TestCase(0.1, 0.047939)] + [TestCase(0.4, 0.169710)] + [TestCase(2.0, 0.487877)] + [TestCase(2.4, 0.521712)] + [TestCase(4.5, 0.600147)] + [TestCase(4.9, 0.606142)] + [TestCase(10.0, 0.630018)] + [TestCase(20.0, 0.635016)] + //[TestCase(100.0, 0.636556)] Needs direct approximation + public void BesselI1MStruveL1Exact(double x, double expected) + { + // Abramowitz/Stegun Table 12.1, 12.2 + AssertHelpers.AlmostEqual(expected, SpecialFunctions.BesselI1MStruveL1(x), 5); + } + } +} diff --git a/src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs b/src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs index bc9ddb92..88de0202 100644 --- a/src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs +++ b/src/UnitTests/SpecialFunctionsTests/SpecialFunctionsTests.cs @@ -3,7 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com +// // Copyright (c) 2009-2010 Math.NET +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -12,8 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index cbae7d4d..9e80ec1c 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -766,6 +766,8 @@ + +