From 1e96e2dcbb53b9f94b7121a9b09a29f6928fd411 Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Mon, 17 Dec 2018 00:38:37 +0530 Subject: [PATCH 01/10] Added Complete Elliptic Functions of 1st and 2nd kind --- src/ellip.jl | 210 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 210 insertions(+) create mode 100644 src/ellip.jl diff --git a/src/ellip.jl b/src/ellip.jl new file mode 100644 index 00000000..6aefc153 --- /dev/null +++ b/src/ellip.jl @@ -0,0 +1,210 @@ +using Base.Math: @horner, libm +using Base.MPFR: ROUNDING_MODE +""" +Computes Complete Elliptic function of 1st kind.------- K(m) ; + +Principal domain : 0<=m<=1 + +Integral: from 0 to pi/2 1 + ----------------------- d(theta) + (1-m*(sin theta)^2)^1/2 + +Using piecewise approximation polynomial as given in + +'Fast Computation of Complete Elliptic Integrals and +Jacobian Elliptic Functions' + by Toshio Fukushima +Link : https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf + +For m<0 , followed: 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]. +""" +function elliptic1(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 + +""" +Computes Complete Elliptic function of 2nd kind ------ E(m) + +Principal domain: 0<=m<=1 + +Integral: from 0 to pi/2: (1 - m* (sin theta)^2)^1/2 d(theta) = E(m) + + +Using piecewise approximation polynomial as given in + +'Fast Computation of Complete Elliptic Integrals and +Jacobian Elliptic Functions' + by Toshio Fukushima +Link : https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf + +For m<0 , followed: 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]. +""" + +function elliptic2(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 = elliptic1(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 (:elliptic1,:elliptic2) + @eval begin + ($f)(x::Float32) = Float32($f(Float64(x))) + ($f)(x::Real) = ($f)(float(x)) + ($f)(a::Float16) = Float16($f(Float32(a))) + ($f)(a::Integer) = ($f(Float64(a))) + #($f)(a::Complex{Float64}) = throw(DomainError(a,"Currently , only (-inf,1] domain is supported")) + end +end + + + + + + From 3d1af2241c765685538ed59ff5d3d79e34e24327 Mon Sep 17 00:00:00 2001 From: Sumegh Roychowdhury <37850881+Sumegh-git@users.noreply.github.com> Date: Mon, 17 Dec 2018 00:44:41 +0530 Subject: [PATCH 02/10] Updated docs with elliptic1 and elliptic2 --- docs/src/index.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/src/index.md b/docs/src/index.md index 5d6f38d4..456f7df6 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`) | +| [`elliptic1(x)`](@ref SpecialFunctions.elliptic1) | [complete elliptic integral of 1st kind](https://en.wikipedia.org/wiki/Elliptic_integral) (i.e evaluates complete elliptic integral of 1st kind at `x`) | +| [`elliptic2(x)`](@ref SpecialFunctions.elliptic2) | [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` | From c770b931f5ba2e90d18733f0ef133d2dfaa75209 Mon Sep 17 00:00:00 2001 From: Sumegh Roychowdhury <37850881+Sumegh-git@users.noreply.github.com> Date: Mon, 17 Dec 2018 10:48:33 +0530 Subject: [PATCH 03/10] Added elliptic1 and elliptic2. --- docs/src/special.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/src/special.md b/docs/src/special.md index b6eb5d2d..9b3002bd 100644 --- a/docs/src/special.md +++ b/docs/src/special.md @@ -44,6 +44,8 @@ SpecialFunctions.besseli SpecialFunctions.besselix SpecialFunctions.besselk SpecialFunctions.besselkx +SpecialFunctions.elliptic1 +SpecialFunctions.elliptic2 SpecialFunctions.eta SpecialFunctions.zeta SpecialFunctions.gamma From 0a26c34fa1c6698004fcd602c6de3b9ad7b0ba9c Mon Sep 17 00:00:00 2001 From: Sumegh Roychowdhury <37850881+Sumegh-git@users.noreply.github.com> Date: Mon, 17 Dec 2018 12:41:53 +0530 Subject: [PATCH 04/10] Updated ellip.jl Made all required changes. --- src/ellip.jl | 72 ++++++++++++++++++++-------------------------------- 1 file changed, 27 insertions(+), 45 deletions(-) diff --git a/src/ellip.jl b/src/ellip.jl index 6aefc153..ccba0276 100644 --- a/src/ellip.jl +++ b/src/ellip.jl @@ -1,23 +1,19 @@ -using Base.Math: @horner, libm -using Base.MPFR: ROUNDING_MODE -""" -Computes Complete Elliptic function of 1st kind.------- K(m) ; - -Principal domain : 0<=m<=1 - -Integral: from 0 to pi/2 1 - ----------------------- d(theta) - (1-m*(sin theta)^2)^1/2 +using Base.Math: @horner -Using piecewise approximation polynomial as given in +#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]. -'Fast Computation of Complete Elliptic Integrals and -Jacobian Elliptic Functions' - by Toshio Fukushima -Link : https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf +""" + elliptic1(x) -For m<0 , followed: 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]. +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 elliptic1(a::Float64) flag = false @@ -92,29 +88,24 @@ function elliptic1(a::Float64) ans = t / sqrt(1.0-a) return ans end - - end -""" -Computes Complete Elliptic function of 2nd kind ------ E(m) - -Principal domain: 0<=m<=1 - -Integral: from 0 to pi/2: (1 - m* (sin theta)^2)^1/2 d(theta) = E(m) -Using piecewise approximation polynomial as given in +#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]. -'Fast Computation of Complete Elliptic Integrals and -Jacobian Elliptic Functions' - by Toshio Fukushima -Link : https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf - -For m<0 , followed: 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]. """ - + elliptic2(x) +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 elliptic2(a::Float64) flag=false if a<0.0 @@ -181,7 +172,7 @@ function elliptic2(a::Float64) 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 = elliptic1(Float64(x)) - #em = km + (pi/2 - km*edm)/kdm + #em = km + (pi/2 - km*edm)/kdm em = (pi/2 + km*hdm)/kdm #to avoid precision loss near 1 t= em end @@ -195,16 +186,7 @@ end for f in (:elliptic1,:elliptic2) @eval begin - ($f)(x::Float32) = Float32($f(Float64(x))) + ($f)(x::AbstractFloat) = throw(MethodError($f,(x,""))) ($f)(x::Real) = ($f)(float(x)) - ($f)(a::Float16) = Float16($f(Float32(a))) - ($f)(a::Integer) = ($f(Float64(a))) - #($f)(a::Complex{Float64}) = throw(DomainError(a,"Currently , only (-inf,1] domain is supported")) end end - - - - - - From 053da6bdcc3983a397fe7186e44a957eb1b3a919 Mon Sep 17 00:00:00 2001 From: Sumegh Roychowdhury <37850881+Sumegh-git@users.noreply.github.com> Date: Mon, 17 Dec 2018 12:49:09 +0530 Subject: [PATCH 05/10] Added wiki and dlmf link for ellip.jl --- src/ellip.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/ellip.jl b/src/ellip.jl index ccba0276..c98d78f9 100644 --- a/src/ellip.jl +++ b/src/ellip.jl @@ -9,6 +9,8 @@ using Base.Math: @horner #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]. +#Wiki : https://en.wikipedia.org/wiki/Elliptic_integral +#DLMF: https://dlmf.nist.gov/19 """ elliptic1(x) @@ -101,6 +103,8 @@ end #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]. +#Wiki : https://en.wikipedia.org/wiki/Elliptic_integral +#DLMF: https://dlmf.nist.gov/19 """ elliptic2(x) From 9cefeb48f9bb81c96557d716ca150084cec9c055 Mon Sep 17 00:00:00 2001 From: Sumegh Roychowdhury <37850881+Sumegh-git@users.noreply.github.com> Date: Mon, 17 Dec 2018 19:15:51 +0530 Subject: [PATCH 06/10] Updated src/SpecialFunctions.jl --- src/SpecialFunctions.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 88c9c52f..a685a9fe 100644 --- a/src/SpecialFunctions.jl +++ b/src/SpecialFunctions.jl @@ -35,6 +35,8 @@ export bessely1, besselyx, dawson, + elliptic1, + elliptic2, 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, :elliptic1, :elliptic2) @eval $(f)(::Missing) = missing end for f in (:beta, :lbeta) From d8659f9287f1d325cb3ec5a28ec7f6456aafea02 Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Tue, 18 Dec 2018 02:35:32 +0530 Subject: [PATCH 07/10] Renamed , Added tests from Wolframalpha --- docs/src/index.md | 4 ++-- docs/src/special.md | 4 ++-- src/SpecialFunctions.jl | 6 +++--- src/ellip.jl | 28 +++++++++++++++++----------- test/runtests.jl | 20 +++++++++++++++++++- 5 files changed, 43 insertions(+), 19 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 456f7df6..e0b7c45f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -19,8 +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`) | -| [`elliptic1(x)`](@ref SpecialFunctions.elliptic1) | [complete elliptic integral of 1st kind](https://en.wikipedia.org/wiki/Elliptic_integral) (i.e evaluates complete elliptic integral of 1st kind at `x`) | -| [`elliptic2(x)`](@ref SpecialFunctions.elliptic2) | [complete elliptic integral of 2nd kind](https://en.wikipedia.org/wiki/Elliptic_integral) (i.e evaluates complete elliptic integral of 2nd kind 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 9b3002bd..88b67414 100644 --- a/docs/src/special.md +++ b/docs/src/special.md @@ -44,8 +44,8 @@ SpecialFunctions.besseli SpecialFunctions.besselix SpecialFunctions.besselk SpecialFunctions.besselkx -SpecialFunctions.elliptic1 -SpecialFunctions.elliptic2 +SpecialFunctions.ellipk +SpecialFunctions.ellipe SpecialFunctions.eta SpecialFunctions.zeta SpecialFunctions.gamma diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index a685a9fe..9dd9200d 100644 --- a/src/SpecialFunctions.jl +++ b/src/SpecialFunctions.jl @@ -35,8 +35,8 @@ export bessely1, besselyx, dawson, - elliptic1, - elliptic2, + ellipk, + ellipe, erf, erfc, erfcinv, @@ -65,7 +65,7 @@ include("gamma.jl") include("deprecated.jl") for f in (:digamma, :erf, :erfc, :erfcinv, :erfcx, :erfi, :erfinv, - :eta, :gamma, :invdigamma, :lfactorial, :lgamma, :trigamma, :elliptic1, :elliptic2) + :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 index c98d78f9..f30952c9 100644 --- a/src/ellip.jl +++ b/src/ellip.jl @@ -9,15 +9,17 @@ using Base.Math: @horner #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]. -#Wiki : https://en.wikipedia.org/wiki/Elliptic_integral -#DLMF: https://dlmf.nist.gov/19 + + """ - elliptic1(x) + 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 elliptic1(a::Float64) +function ellipk(a::Float64) flag = false if a<0.0 x=a/(a-1) #dealing with negative args @@ -103,14 +105,16 @@ end #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]. -#Wiki : https://en.wikipedia.org/wiki/Elliptic_integral -#DLMF: https://dlmf.nist.gov/19 + """ - elliptic2(x) + 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 elliptic2(a::Float64) +function ellipe(a::Float64) flag=false if a<0.0 x=a/(a-1) #for dealing with negative args @@ -175,7 +179,7 @@ function elliptic2(a::Float64) 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 = elliptic1(Float64(x)) + 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 @@ -188,9 +192,11 @@ function elliptic2(a::Float64) end end -for f in (:elliptic1,:elliptic2) +for f in (:ellipk,:ellipe) @eval begin - ($f)(x::AbstractFloat) = throw(MethodError($f,(x,""))) + ($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..decc0451 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -61,7 +61,25 @@ 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 + @test ellipk(0.92) ≈ 2.683551406315229344 + @test ellipk(0.5) ≈ 1.854074677301371918 + @test ellipk(0.01) ≈ 1.57474556151735595 + @test ellipk(0.45) ≈ 1.81388393681698264 + @test ellipk(-0.5) ≈ 1.41573720842595619 rtol=2*eps() + @test ellipk(0.75) ≈ 2.15651564749964323 rtol=2*eps() + @test ellipe(0) ≈ 1.570796326794896619231322 + @test ellipe(0.8) ≈ 1.17848992432783852 + @test ellipe(0.5) ≈ 1.3506438810476755 rtol=2*eps() + @test ellipe(0.01) ≈ 1.5668619420216682 + @test ellipe(0.99) ≈ 1.0159935450252239 + @test ellipe(-0.1) ≈ 1.6093590249375295 rtol=2*eps() + @test ellipe(0.3) ≈ 1.4453630644126652 rtol=2*eps() + @test_throws MethodError ellipk(BigFloat(0.5)) + @test_throws MethodError ellipe(BigFloat(-1)) +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] From 88e0e774b061f776b7033e62c372487d4fb9ae85 Mon Sep 17 00:00:00 2001 From: Sumegh Roychowdhury <37850881+Sumegh-git@users.noreply.github.com> Date: Tue, 18 Dec 2018 03:10:03 +0530 Subject: [PATCH 08/10] Updated runtests.jl --- test/runtests.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index decc0451..6c8d82fe 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -63,18 +63,18 @@ relerrc(z, x) = max(relerr(real(z),real(x)), relerr(imag(z),imag(x))) end @testset "elliptic integrals" begin #Computed using Wolframalpha EllipticK and EllipticE functions. - @test ellipk(0) ≈ 1.570796326794896619231322 - @test ellipk(0.92) ≈ 2.683551406315229344 - @test ellipk(0.5) ≈ 1.854074677301371918 - @test ellipk(0.01) ≈ 1.57474556151735595 - @test ellipk(0.45) ≈ 1.81388393681698264 + @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 ellipe(0) ≈ 1.570796326794896619231322 - @test ellipe(0.8) ≈ 1.17848992432783852 + @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 - @test ellipe(0.99) ≈ 1.0159935450252239 + @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_throws MethodError ellipk(BigFloat(0.5)) From 930a084745f68029c2d826e600c9574b7137c32e Mon Sep 17 00:00:00 2001 From: Sumegh Roychowdhury <37850881+Sumegh-git@users.noreply.github.com> Date: Wed, 19 Dec 2018 04:07:54 +0530 Subject: [PATCH 09/10] Updated runtests.jl (1) --- test/runtests.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 6c8d82fe..70770cce 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -70,6 +70,17 @@ end @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.15) ≈ 1.635256732264579 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 ellipk(Float32(0.899)) ≈ 2.573411 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() @@ -77,8 +88,11 @@ end @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)) 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]] From a6423822aa4932dacf6945fa37583e5f7c6235da Mon Sep 17 00:00:00 2001 From: Sumegh Roychowdhury <37850881+Sumegh-git@users.noreply.github.com> Date: Wed, 19 Dec 2018 04:18:00 +0530 Subject: [PATCH 10/10] Final tests. --- test/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 70770cce..92bc8f76 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -70,11 +70,10 @@ end @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.15) ≈ 1.635256732264579 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 ellipk(Float32(0.899)) ≈ 2.573411 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() @@ -93,6 +92,7 @@ end @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]]