diff --git a/docs/src/index.md b/docs/src/index.md index 5d6f38d4..e0b7c45f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -19,6 +19,8 @@ libraries. | [`invdigamma(x)`](@ref SpecialFunctions.invdigamma) | [invdigamma function](http://bariskurt.com/calculating-the-inverse-of-digamma-function/) (i.e. inverse of `digamma` function at `x` using fixed-point iteration algorithm) | | [`trigamma(x)`](@ref SpecialFunctions.trigamma) | [trigamma function](https://en.wikipedia.org/wiki/Trigamma_function) (i.e the logarithmic second derivative of `gamma` at `x`) | | [`polygamma(m,x)`](@ref SpecialFunctions.polygamma) | [polygamma function](https://en.wikipedia.org/wiki/Polygamma_function) (i.e the (m+1)-th derivative of the `lgamma` function at `x`) | +| [`ellipk(x)`](@ref SpecialFunctions.ellipk) | [complete elliptic integral of 1st kind](https://en.wikipedia.org/wiki/Elliptic_integral) (i.e evaluates complete elliptic integral of 1st kind at `x`) | +| [`ellipe(x)`](@ref SpecialFunctions.ellipe) | [complete elliptic integral of 2nd kind](https://en.wikipedia.org/wiki/Elliptic_integral) (i.e evaluates complete elliptic integral of 2nd kind at `x`) | | [`eta(x)`](@ref SpecialFunctions.eta) | [Dirichlet eta function](https://en.wikipedia.org/wiki/Dirichlet_eta_function) at `x` | | [`zeta(x)`](@ref SpecialFunctions.zeta) | [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function) at `x` | | [`airyai(z)`](@ref SpecialFunctions.airyai) | [Airy Ai function](https://en.wikipedia.org/wiki/Airy_function) at `z` | diff --git a/docs/src/special.md b/docs/src/special.md index b6eb5d2d..88b67414 100644 --- a/docs/src/special.md +++ b/docs/src/special.md @@ -44,6 +44,8 @@ SpecialFunctions.besseli SpecialFunctions.besselix SpecialFunctions.besselk SpecialFunctions.besselkx +SpecialFunctions.ellipk +SpecialFunctions.ellipe SpecialFunctions.eta SpecialFunctions.zeta SpecialFunctions.gamma diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 88c9c52f..9dd9200d 100644 --- a/src/SpecialFunctions.jl +++ b/src/SpecialFunctions.jl @@ -35,6 +35,8 @@ export bessely1, besselyx, dawson, + ellipk, + ellipe, erf, erfc, erfcinv, @@ -57,12 +59,13 @@ export include("bessel.jl") include("erf.jl") +include("ellip.jl") include("sincosint.jl") include("gamma.jl") include("deprecated.jl") for f in (:digamma, :erf, :erfc, :erfcinv, :erfcx, :erfi, :erfinv, - :eta, :gamma, :invdigamma, :lfactorial, :lgamma, :trigamma) + :eta, :gamma, :invdigamma, :lfactorial, :lgamma, :trigamma, :ellipk, :ellipe) @eval $(f)(::Missing) = missing end for f in (:beta, :lbeta) diff --git a/src/ellip.jl b/src/ellip.jl new file mode 100644 index 00000000..f30952c9 --- /dev/null +++ b/src/ellip.jl @@ -0,0 +1,202 @@ +using Base.Math: @horner + +#Using piecewise approximation polynomial as given in +#'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions' +# Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr +# DOI 10.1007/s10569-009-9228-z +#Link : https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf +# +#For m<0 , followed Fukushima, Toshio. (2014). Precise, compact, and fast computation of complete elliptic integrals by piecewise minimax rational function approximation. Journal of Computational and Applied Mathematics. 282. 10.13140/2.1.1946.6245. +# Link: https://www.researchgate.net/profile/Toshio_Fukushima/publication/267330394_Precise_compact_and_fast_computation_of_complete_elliptic_integrals_by_piecewise_minimax_rational_function_approximation/links/544b81a40cf2d6347f43074f/Precise-compact-and-fast-computation-of-complete-elliptic-integrals-by-piecewise-minimax-rational-function-approximation.pdf?origin=publication_detail +#Also suggested in this paper that we should consider domain only from (-inf,1]. + + + +""" + ellipk(x) + + DLMF : https://dlmf.nist.gov/19.2#E4 , https://dlmf.nist.gov/19.2#E8 + Wiki : https://en.wikipedia.org/wiki/Elliptic_integral +Computes Complete Elliptic Integral of 1st kind at `x`-> K(x)--- given by: ``\\int_{0}^{\\pi/2} \\frac{1}{\\sqrt{1 - x(\\sin \\theta )^{2}}} d\\theta`` +""" +function ellipk(a::Float64) + flag = false + if a<0.0 + x=a/(a-1) #dealing with negative args + flag=true + elseif a>=0.0 + x=a + flag=false + end + + if x==0.0 + return pi/2 + + elseif x==1.0 + return Inf + + elseif x>1.0 + throw(DomainError(x,"`x` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]")) + + elseif 0.0 <= x < 0.1 #Table 2 from paper + t= x-0.05 + t= @horner(t, 1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 , 0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 , 0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 , 0.085842591595413900 , 0.081541118718303215) + + elseif 0.1 <= x <0.2 #Table 3 + t=x-0.15 + t= @horner(t , 1.635256732264579992 , 0.471190626148732291 , 0.309728410831499587 , 0.252208311773135699 , 0.226725623219684650 , 0.215774446729585976 , 0.213108771877348910 , 0.216029124605188282 , 0.223255831633057896 , 0.234180501294209925 , 0.248557682972264071 , 0.266363809892617521) + + elseif 0.2 <= x < 0.3 #Table 4 + t=x-0.25 + t= @horner(t , 1.685750354812596043 , 0.541731848613280329 , 0.401524438390690257 , 0.369642473420889090 , 0.376060715354583645 , 0.405235887085125919 , 0.453294381753999079 , 0.520518947651184205 , 0.609426039204995055 , 0.724263522282908870 , 0.871013847709812357 , 1.057652872753547036) + + elseif 0.3 <= x < 0.4 #Table 5 + t= x-0.35 + t= @horner(t , 1.744350597225613243 , 0.634864275371935304 , 0.539842564164445538 , 0.571892705193787391 , 0.670295136265406100 , 0.832586590010977199 , 1.073857448247933265 , 1.422091460675497751 , 1.920387183402304829 , 2.632552548331654201 , 3.652109747319039160 , 5.115867135558865806 , 7.224080007363877411) + + elseif 0.4 <= x < 0.5 #Table 6 + t=x-0.45 + t= @horner(t, 1.813883936816982644 , 0.763163245700557246 , 0.761928605321595831 , 0.951074653668427927 , 1.315180671703161215 , 1.928560693477410941 , 2.937509342531378755 , 4.594894405442878062 , 7.330071221881720772 , 11.87151259742530180 , 19.45851374822937738 , 32.20638657246426863 , 53.73749198700554656 , 90.27388602940998849) + + elseif 0.5 <= x < 0.6 #Table 7 + t= x-0.55 + t= @horner(t , 1.898924910271553526 , 0.950521794618244435 , 1.151077589959015808 , 1.750239106986300540 , 2.952676812636875180 , 5.285800396121450889 , 9.832485716659979747 , 18.78714868327559562 , 36.61468615273698145 , 72.45292395127771801 , 145.1079577347069102 , 293.4786396308497026 , 598.3851815055010179 , 1228.420013075863451 , 2536.529755382764488) + elseif 0.6 <= x < 0.7 #Table 8 + t=x-0.65 + t =@horner(t , 2.007598398424376302 , 1.248457231212347337 , 1.926234657076479729 , 3.751289640087587680 , 8.119944554932045802 , 18.66572130873555361 , 44.60392484291437063 , 109.5092054309498377 , 274.2779548232413480,697.5598008606326163 , 1795.716014500247129 , 4668.381716790389910 , 12235.76246813664335 , 32290.17809718320818 , 85713.07608195964685 , 228672.1890493117096 , 612757.2711915852774) + + elseif 0.7 <= x <0.8 #Table 9 + t=x-0.75 + t= @horner(t, 2.156515647499643235 , 1.791805641849463243 , 3.826751287465713147 , 10.38672468363797208 , 31.40331405468070290 , 100.9237039498695416 , 337.3268282632272897 , 1158.707930567827917 , 4060.990742193632092 , 14454.00184034344795 , 52076.66107599404803 , 189493.6591462156887 , 695184.5762413896145 , 2567994.048255284686 , 9541921.966748386322 , 35634927.44218076174 , 133669298.4612040871 , 503352186.6866284541 , 1901975729.538660119 , 7208915015.330103756) + + elseif 0.8 <= x <0.85 #Table 10 + t = x-0.825 + t =@horner(t , 2.318122621712510589 , 2.616920150291232841 , 7.897935075731355823 , 30.50239715446672327 , 131.4869365523528456 , 602.9847637356491617 , 2877.024617809972641 , 14110.51991915180325 , 70621.44088156540229 , 358977.2665825309926 , 1847238.263723971684 , 9600515.416049214109 , 50307677.08502366879 , 265444188.6527127967 , 1408862325.028702687 , 7515687935.373774627) + + elseif 0.85 <= x <0.9 #Table 11 + t= x-0.875 + t =@horner(t, 2.473596173751343912 , 3.727624244118099310 , 15.60739303554930496 , 84.12850842805887747 , 506.9818197040613935 , 3252.277058145123644 , 21713.24241957434256 , 149037.0451890932766 , 1043999.331089990839 , 7427974.817042038995 , 53503839.67558661151 , 389249886.9948708474 , 2855288351.100810619 , 21090077038.76684053 , 156699833947.7902014 , 1170222242422.439893 , 8777948323668.937971 , 66101242752484.95041 , 499488053713388.7989 , 37859743397240299.20) + + elseif x>=0.9 + td= 1-x + td1=td-0.05 + qd = @horner(td, 0.0, (1.0/16.0), (1.0/32.0), (21.0/1024.0), (31.0/2048.0), (6257.0/524288.0), (10293.0/1048576.0), (279025.0/33554432.0), (483127.0/67108864.0), (435506703.0/68719476736.0), (776957575.0/137438953472.0) , (22417045555.0/4398046511104.0) , (40784671953.0/8796093022208.0) , (9569130097211.0/2251799813685248.0) , (17652604545791.0/4503599627370496.0)) + + kdm = @horner(td1 , 1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 , 0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 , 0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 , 0.085842591595413900 , 0.081541118718303215) + km = -Base.log(qd) * (kdm/pi) + t = km + end + if flag == false + return t + elseif flag == true + ans = t / sqrt(1.0-a) + return ans + end +end + + + +#Using piecewise approximation polynomial as given in +#'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions' +# Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr +# DOI 10.1007/s10569-009-9228-z +#Link : https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf +# +#For m<0 , followed Fukushima, Toshio. (2014). Precise, compact, and fast computation of complete elliptic integrals by piecewise minimax rational function approximation. Journal of Computational and Applied Mathematics. 282. 10.13140/2.1.1946.6245. +# Link: https://www.researchgate.net/profile/Toshio_Fukushima/publication/267330394_Precise_compact_and_fast_computation_of_complete_elliptic_integrals_by_piecewise_minimax_rational_function_approximation/links/544b81a40cf2d6347f43074f/Precise-compact-and-fast-computation-of-complete-elliptic-integrals-by-piecewise-minimax-rational-function-approximation.pdf?origin=publication_detail +#Also suggested in this paper that we should consider domain only from (-inf,1]. + + +""" + ellipe(x) + + DLMF : https://dlmf.nist.gov/19.2#E5 , https://dlmf.nist.gov/19.2#E8 + Wiki : https://en.wikipedia.org/wiki/Elliptic_integral +Computes the Complete Elliptic Integral of 2nd kind at `x` ->E(x)---gven by: ``\\int_{0}^{\\pi/2} \\sqrt{1 - x(\\sin \\theta )^{2}} d\\theta`` +""" +function ellipe(a::Float64) + flag=false + if a<0.0 + x=a/(a-1) #for dealing with negative args + flag=true + elseif a>=0.0 + x=a + flag=false + + end + + if x==0.0 + return pi/2 + elseif x==1.0 + return 1.0 + + elseif x>1.0 + throw(DomainError(x,"`x` must lie between -Inf and 1 ---- Domain : (-inf,1.0]")) + + elseif 0.0 <= x < 0.1 #Table 2 from paper + t= x-0.05 + t= @horner(t , +1.550973351780472328 ,-0.400301020103198524 ,-0.078498619442941939 ,-0.034318853117591992 ,-0.019718043317365499 ,-0.013059507731993309 ,-0.009442372874146547 ,-0.007246728512402157 ,-0.005807424012956090 ,-0.004809187786009338) + + elseif 0.1 <= x < 0.2 #Table 3 + t=x-0.15 + t= @horner(t , +1.510121832092819728 , -0.417116333905867549 , -0.090123820404774569 , -0.043729944019084312, -0.027965493064761785 , -0.020644781177568105 , -0.016650786739707238 , -0.014261960828842520 , -0.012759847429264803 , -0.011799303775587354 , -0.011197445703074968) + + elseif 0.2 <= x < 0.3 #Table 4 + t=x-0.25 + t= @horner(t , 1.467462209339427155 , -0.436576290946337775 , -0.105155557666942554 , -0.057371843593241730 , -0.041391627727340220 , -0.034527728505280841 , -0.031495443512532783 , -0.030527000890325277 , -0.030916984019238900 , -0.032371395314758122 , -0.034789960386404158) + + elseif 0.3 <= x < 0.4 #Table 5 + t= x-0.35 + t= @horner(t, +1.422691133490879171, -0.459513519621048674, -0.125250539822061878, -0.078138545094409477, -0.064714278472050002, -0.062084339131730311, -0.065197032815572477,-0.072793895362578779, -0.084959075171781003, -0.102539850131045997, -0.127053585157696036, -0.160791120691274606) + + elseif 0.4 <= x < 0.5 #Table 6 + t=x-0.45 + t= @horner(t ,+1.375401971871116291 ,-0.487202183273184837 ,-0.153311701348540228 ,-0.111849444917027833 ,-0.108840952523135768 ,-0.122954223120269076 ,-0.152217163962035047 ,-0.200495323642697339 ,-0.276174333067751758 ,-0.393513114304375851 ,-0.575754406027879147 ,-0.860523235727239756 ,-1.308833205758540162) + + elseif 0.5 <= x < 0.6 #Table 7 + t= x-0.55 + t= @horner(t ,+1.325024497958230082 ,-0.521727647557566767 ,-0.194906430482126213 ,-0.171623726822011264 ,-0.202754652926419141 ,-0.278798953118534762 ,-0.420698457281005762 ,-0.675948400853106021 ,-1.136343121839229244 ,-1.976721143954398261 ,-3.531696773095722506 ,-6.446753640156048150 ,-11.97703130208884026) + elseif 0.6 <= x < 0.7 #Table 8 + t=x-0.65 + t= @horner(t, +1.270707479650149744, -0.566839168287866583, -0.262160793432492598, -0.292244173533077419, -0.440397840850423189, -0.774947641381397458, -1.498870837987561088, -3.089708310445186667, -6.667595903381001064, -14.89436036517319078, -34.18120574251449024, -80.15895841905397306, -191.3489480762984920, -463.5938853480342030, -1137.380822169360061) + + elseif 0.7 <= x < 0.8 #Table 9 + t=x-0.75 + t= @horner(t, +1.211056027568459525, -0.630306413287455807, -0.387166409520669145, -0.592278235311934603, -1.237555584513049844, -3.032056661745247199, -8.181688221573590762, -23.55507217389693250, -71.04099935893064956, -221.8796853192349888, -712.1364793277635425, -2336.125331440396407, -7801.945954775964673, -26448.19586059191933, -90799.48341621365251, -315126.0406449163424, -1104011.344311591159) + + elseif 0.8 <= x < 0.85 #Table 10 + t = x-0.825 + t= @horner(t, +1.161307152196282836, -0.701100284555289548, -0.580551474465437362, -1.243693061077786614, -3.679383613496634879, -12.81590924337895775, -49.25672530759985272, -202.1818735434090269, -869.8602699308701437, -3877.005847313289571, -17761.70710170939814, -83182.69029154232061, -396650.4505013548170, -1920033.413682634405) + + elseif 0.85 <= x < 0.9 #Table 11 + t= x-0.875 + t = @horner(t, +1.124617325119752213, -0.770845056360909542, -0.844794053644911362, -2.490097309450394453, -10.23971741154384360, -49.74900546551479866, -267.0986675195705196, -1532.665883825229947, -9222.313478526091951, -57502.51612140314030, -368596.1167416106063, -2415611.088701091428, -16120097.81581656797, -109209938.5203089915, -749380758.1942496220, -5198725846.725541393, -36409256888.12139973) + + elseif x>=0.9 + td= 1-x + td1=td-0.05 + + kdm = @horner(td1 , 1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 , 0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711,0.108938811574293531 , 0.098853409871592910 ,0.091439629201749751 ,0.085842591595413900 ,0.081541118718303215) + edm = @horner(td1 ,+1.550973351780472328 ,-0.400301020103198524 ,-0.078498619442941939 ,-0.034318853117591992,-0.019718043317365499 ,-0.013059507731993309 ,-0.009442372874146547 ,-0.007246728512402157 ,-0.005807424012956090 ,-0.004809187786009338) + hdm = kdm -edm + km = ellipk(Float64(x)) + #em = km + (pi/2 - km*edm)/kdm + em = (pi/2 + km*hdm)/kdm #to avoid precision loss near 1 + t= em + end + if flag == false + return t + elseif flag == true + return t * sqrt(1.0-a) + + end +end + +for f in (:ellipk,:ellipe) + @eval begin + ($f)(x::Float16) = Float16(($f)(Float64(x))) + ($f)(x::Float32) = Float32(($f)(Float64(x))) + ($f)(x::Real) = ($f)(float(x)) + ($f)(x::AbstractFloat) = throw(MethodError($f,(x,""))) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 84312173..92bc8f76 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -61,7 +61,39 @@ relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x))) @test_throws MethodError erf(big(1.0)*im) @test_throws MethodError erfi(big(1.0)) end - +@testset "elliptic integrals" begin +#Computed using Wolframalpha EllipticK and EllipticE functions. + @test ellipk(0) ≈ 1.570796326794896619231322 rtol=2*eps() + @test ellipk(0.92) ≈ 2.683551406315229344 rtol=2*eps() + @test ellipk(0.5) ≈ 1.854074677301371918 rtol=2*eps() + @test ellipk(0.01) ≈ 1.57474556151735595 rtol=2*eps() + @test ellipk(0.45) ≈ 1.81388393681698264 rtol=2*eps() + @test ellipk(-0.5) ≈ 1.41573720842595619 rtol=2*eps() + @test ellipk(0.75) ≈ 2.15651564749964323 rtol=2*eps() + @test ellipk(0.17) ≈ 1.6448064907988806 rtol=2*eps() + @test ellipk(0.25) ≈ 1.685750354812596 rtol=2*eps() + @test ellipk(0.69) ≈ 2.0608816467301313 rtol=2*eps() + @test ellipk(0.84) ≈ 2.3592635547450067 rtol=2*eps() + @test ellipe(0.15) ≈ 1.5101218320928197 rtol=2*eps() + @test ellipe(0.21) ≈ 1.4847605813318776 rtol=2*eps() + @test ellipe(0.42) ≈ 1.3898829914929717 rtol=2*eps() + @test ellipe(0.66) ≈ 1.2650125751607508 rtol=2*eps() + @test ellipe(0.76) ≈ 1.2047136418292115 rtol=2*eps() + @test ellipe(0.865) ≈ 1.1322436887003925 rtol=2*eps() + @test ellipe(0) ≈ 1.570796326794896619231322 rtol=2*eps() + @test ellipe(0.8) ≈ 1.17848992432783852 rtol=2*eps() + @test ellipe(0.5) ≈ 1.3506438810476755 rtol=2*eps() + @test ellipe(0.01) ≈ 1.5668619420216682 rtol=2*eps() + @test ellipe(0.99) ≈ 1.0159935450252239 rtol=2*eps() + @test ellipe(-0.1) ≈ 1.6093590249375295 rtol=2*eps() + @test ellipe(0.3) ≈ 1.4453630644126652 rtol=2*eps() + @test ellipe(1.0) ≈ 1.00 + @test ellipk(1.0)==Inf + @test_throws MethodError ellipk(BigFloat(0.5)) + @test_throws MethodError ellipe(BigFloat(-1)) + @test_throws DomainError ellipe(Float16(2.0)) + @test_throws DomainError ellipe(Float32(2.5)) +end @testset "sine and cosine integrals" begin # Computed via wolframalpha.com: SinIntegral[SetPrecision[Table[x,{x, 1,20,1}],20]] and CosIntegral[SetPrecision[Table[x,{x, 1,20,1}],20]] sinintvals = [0.9460830703671830149, 1.605412976802694849, 1.848652527999468256, 1.75820313894905306, 1.54993124494467414, 1.4246875512805065, 1.4545966142480936, 1.5741868217069421, 1.665040075829602, 1.658347594218874, 1.578306806945727416, 1.504971241526373371, 1.499361722862824564, 1.556211050077665054, 1.618194443708368739, 1.631302268270032886, 1.590136415870701122, 1.536608096861185462, 1.518630031769363932, 1.548241701043439840]