diff --git a/docs/src/functions_list.md b/docs/src/functions_list.md index 5890373f..203603ee 100644 --- a/docs/src/functions_list.md +++ b/docs/src/functions_list.md @@ -143,6 +143,7 @@ jinc ```@docs ellipk ellipe +ellipke ``` ## Zeta and Related Functions diff --git a/docs/src/functions_overview.md b/docs/src/functions_overview.md index 6cea9b02..2e203c99 100644 --- a/docs/src/functions_overview.md +++ b/docs/src/functions_overview.md @@ -108,6 +108,7 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi |:-------- |:----------- | | [`ellipk(m)`](@ref SpecialFunctions.ellipk) | [complete elliptic integral of 1st kind](https://en.wikipedia.org/wiki/Elliptic_integral#Notational_variants) ``K(m)`` | | [`ellipe(m)`](@ref SpecialFunctions.ellipe) | [complete elliptic integral of 2nd kind](https://en.wikipedia.org/wiki/Elliptic_integral#Complete_elliptic_integral_of_the_second_kind) ``E(m)`` | +| [`ellipke(m)`](@ref SpecialFunctions.ellipke) | complete elliptic integrals ``K(m)`` and ``E(m)`` simultaneously | ## Zeta and Related Functions diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 8b6c289e..fd2fa30b 100644 --- a/src/SpecialFunctions.jl +++ b/src/SpecialFunctions.jl @@ -47,6 +47,7 @@ export dawson, ellipk, ellipe, + ellipke, erf, erfc, erfcinv, diff --git a/src/ellip.jl b/src/ellip.jl index 153323ea..7b033782 100644 --- a/src/ellip.jl +++ b/src/ellip.jl @@ -1,8 +1,6 @@ using Base.Math: @horner - - @doc raw""" ellipk(m) @@ -35,89 +33,96 @@ As suggested in this paper, the domain is restricted to ``(-\infty,1]``. """ ellipk(m::Real) = _ellipk(float(m)) +# Wrapper function that handles negative m and special cases function _ellipk(m::Float64) isnan(m) && return NaN - (m === -Inf) && return 0.0 - flag_is_m_neg = m < 0.0 - if flag_is_m_neg - x = m / (m-1) #dealing with negative args + if m >= 0 + return _ellipk_nonneg(m) else - x = m + m === -Inf && return 0.0 + # Negative m transformation: x = m/(m-1), K(m) = K(x)/sqrt(1-m) + x = m / (m - 1) + # When m is extremely negative, x rounds to 1.0; K(m) → 0 as m → -∞ + x == 1.0 && return 0.0 + return _ellipk_nonneg(x) / sqrt(1 - m) end +end - if x == 0.0 - return Float64(halfπ) +# Implementation of ellipk for m >= 0 +@inline function _ellipk_nonneg(m::Float64) + # Edge cases + m == 0.0 && return Float64(halfπ) + m == 1.0 && return Inf + m > 1.0 && throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]")) - elseif x == 1.0 - return flag_is_m_neg ? 0.0 : Inf + # Compute index and delegate to table lookup + idx = unsafe_trunc(Int, m * 10) + return @inbounds _ellipk_horner_table(idx, m) +end - elseif x > 1.0 - throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]")) +# Pure polynomial lookup for ellipk using linear if-else structure +# LLVM optimizes this to efficient switch/jump/binary-search table +@inline function _ellipk_horner_table(idx::Int, m::Float64) + @boundscheck (0 <= idx <= 9) || throw(ArgumentError("idx=$idx out of valid range 0:9")) - elseif 0.0 < x < 0.1 #Table 2 from paper - t = x-0.05 - t = @horner(t, + if idx == 0 # 0.0 < m < 0.1, Table 2 + t = m - 0.05 + return @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 , + elseif idx == 1 # 0.1 <= m < 0.2, Table 3 + t = m - 0.15 + return @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 , + elseif idx == 2 # 0.2 <= m < 0.3, Table 4 + t = m - 0.25 + return @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 , + elseif idx == 3 # 0.3 <= m < 0.4, Table 5 + t = m - 0.35 + return @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, + elseif idx == 4 # 0.4 <= m < 0.5, Table 6 + t = m - 0.45 + return @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 , + elseif idx == 5 # 0.5 <= m < 0.6, Table 7 + t = m - 0.55 + return @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 , + elseif idx == 6 # 0.6 <= m < 0.7, Table 8 + t = m - 0.65 + return @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, + elseif idx == 7 # 0.7 <= m < 0.8, Table 9 + t = m - 0.75 + return @horner(t, 2.156515647499643235 , 1.791805641849463243 , 3.826751287465713147 , 10.38672468363797208 , 31.40331405468070290 , 100.9237039498695416 , 337.3268282632272897 , 1158.707930567827917 , 4060.990742193632092 , @@ -125,57 +130,45 @@ function _ellipk(m::Float64) 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) - - else - @assert 0.9 <= x < 1 - td = 1-x - td1 = td-0.05 + elseif idx == 8 # 0.8 <= m < 0.9 + if m < 0.85 # Table 10, 0.8 <= m < 0.85 + t = m - 0.825 + return @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) + else # Table 11, 0.85 <= m < 0.9 + t = m - 0.875 + return @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) + end + else # 0.9 <= m < 1.0, special log case + td = 1 - m + 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 , + (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 * invπ) - t = km - end - - if flag_is_m_neg - ans = t / sqrt(1.0-m) - return ans - else - return t + return -Base.log(qd) * (kdm * invπ) end end - @doc raw""" ellipe(m) @@ -208,141 +201,236 @@ As suggested in this paper, the domain is restricted to ``(-\infty,1]``. """ ellipe(m::Real) = _ellipe(float(m)) +# Wrapper function that handles negative m and special cases function _ellipe(m::Float64) isnan(m) && return NaN - (m === -Inf) && return Inf - - flag_is_m_neg = m < 0.0 - if flag_is_m_neg - x = m / (m-1) #dealing with negative args + + if m >= 0 + return _ellipe_nonneg(m) else - x = m + m === -Inf && return Inf + # Negative m transformation: x = m/(m-1), E(m) = E(x)*sqrt(1-m) + x = m / (m - 1) + # When m is extremely negative, x rounds to 1.0; E(m) → ∞ as m → -∞ + x == 1.0 && return Inf + return _ellipe_nonneg(x) * sqrt(1 - m) end +end - if x == 0.0 - return Float64(halfπ) - elseif x == 1.0 - return flag_is_m_neg ? Inf : 1.0 +# Core implementation of ellipe for m >= 0 +@inline function _ellipe_nonneg(m::Float64) + # Edge cases + m == 0.0 && return Float64(halfπ) + m == 1.0 && return 1.0 + m > 1.0 && throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]")) - elseif x > 1.0 - throw(DomainError(m,"`m` must lie between -Inf and 1 ---- Domain : (-inf,1.0]")) + # Compute index and delegate to table lookup + idx = unsafe_trunc(Int, m * 10) + return @inbounds _ellipe_horner_table(idx, m) +end + +# Pure polynomial lookup for ellipe using linear if-else structure +# LLVM optimizes this to efficient switch/jump table +@inline function _ellipe_horner_table(idx::Int, m::Float64) + @boundscheck (0 <= idx <= 9) || throw(ArgumentError("idx=$idx out of valid range 0:9")) - elseif 0.0 < x < 0.1 #Table 2 from paper - t = x-0.05 - t = @horner(t , + if idx == 0 # 0.0 < m < 0.1, Table 2 + t = m - 0.05 + return @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 , + elseif idx == 1 # 0.1 <= m < 0.2, Table 3 + t = m - 0.15 + return @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 , + elseif idx == 2 # 0.2 <= m < 0.3, Table 4 + t = m - 0.25 + return @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, + elseif idx == 3 # 0.3 <= m < 0.4, Table 5 + t = m - 0.35 + return @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 , + elseif idx == 4 # 0.4 <= m < 0.5, Table 6 + t = m - 0.45 + return @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 , + elseif idx == 5 # 0.5 <= m < 0.6, Table 7 + t = m - 0.55 + return @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, + elseif idx == 6 # 0.6 <= m < 0.7, Table 8 + t = m - 0.65 + return @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, + elseif idx == 7 # 0.7 <= m < 0.8, Table 9 + t = m - 0.75 + return @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) - - else - @assert 0.9 <= x < 1.0 - td = 1-x - td1 = td-0.05 - - kdm = @horner(td1 , + elseif idx == 8 # 0.8 <= m < 0.9 + if m < 0.85 # Table 10, 0.8 <= m < 0.85 + t = m - 0.825 + return @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) + else # Table 11, 0.85 <= m < 0.9 + t = m - 0.875 + return @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) + end + else # 0.9 <= m < 1.0, special case using ellipk + td = 1 - m + 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 , + 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 = (halfπ + km*hdm) / kdm #to avoid precision loss near 1 - t = em + km = @inbounds _ellipk_horner_table(9, m) + return (halfπ + km * hdm) / kdm end - if flag_is_m_neg - return t * sqrt(1.0-m) +end + +@doc raw""" + ellipke(m) + +Computes both Complete Elliptic Integrals ``K(m)`` and ``E(m)`` +simultaneously, returning them as a tuple `(K, E)`. + +This is more efficient than calling `ellipk(m)` and `ellipe(m)` +separately when both values are needed. +See also: [`ellipk(m)`](@ref), [`ellipe(m)`](@ref). + +# Arguments +- `m`: parameter ``m``, restricted to the domain ``(-\infty, 1]`` + +# Returns +- `(K, E)::Tuple`: Complete elliptic integrals of first and second kind +""" +ellipke(m::Real) = _ellipke(float(m)) + +# Wrapper function that handles negative m and special cases +function _ellipke(m::Float64) + isnan(m) && return (NaN, NaN) + + if m >= 0 + return _ellipke_nonneg(m) else - return t + m === -Inf && return (0.0, Inf) + # Negative m transformation: x = m/(m-1) + # K(m) = K(x) / sqrt(1-m), E(m) = E(x) * sqrt(1-m) + x = m / (m - 1) + # When m is extremely negative, x rounds to 1.0 + x == 1.0 && return (0.0, Inf) + K, E = _ellipke_nonneg(x) + sqrt_1_m = sqrt(1 - m) + return (K / sqrt_1_m, E * sqrt_1_m) + end +end + +# Core implementation of combined ellipke for m >= 0 +# Returns (K, E) tuple, computing both simultaneously for efficiency +@inline function _ellipke_nonneg(m::Float64) + # Edge cases + m == 0.0 && return (Float64(halfπ), Float64(halfπ)) + m == 1.0 && return (Inf, 1.0) + m > 1.0 && throw(DomainError(m, "`m` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]")) + + # Compute index and delegate to table lookup + idx = unsafe_trunc(Int, m * 10) + + # Special case for idx == 9 (m >= 0.9): compute K and E together + # to share the expensive log calculation + if idx == 9 + return _ellipke_near_one(m) end + + return @inbounds (_ellipk_horner_table(idx, m), _ellipe_horner_table(idx, m)) +end + +# Specialized path for m >= 0.9 where K and E share the log calculation +@inline function _ellipke_near_one(m::Float64) + td = 1 - m + 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) + edm = @horner(td1, + +1.550973351780472328 , -0.400301020103198524 , -0.078498619442941939 , + -0.034318853117591992 , -0.019718043317365499 , -0.013059507731993309 , + -0.009442372874146547 , -0.007246728512402157 , -0.005807424012956090 , + -0.004809187786009338) + # K calculation (single log call) + K = -Base.log(qd) * (kdm * invπ) + # E calculation reusing K + hdm = kdm - edm + E = (halfπ + K * hdm) / kdm + return (K, E) end # Support for Float32 and Float16 for internalf in (:_ellipk, :_ellipe), T in (:Float16, :Float32) @eval $internalf(x::$T) = $T($internalf(Float64(x))) end + +for T in (:Float16, :Float32) + @eval function _ellipke(x::$T) + K, E = _ellipke(Float64(x)) + return ($T(K), $T(E)) + end +end + +# Missing propagation for ellipke +# Returns (missing, missing) to preserve tuple structure, following Base.sincos convention +# Usage: any(ismissing, ellipke(x)) to check for missing +ellipke(::Missing) = (missing, missing) diff --git a/test/ellip.jl b/test/ellip.jl index b22c4fd7..5b0b7dce 100644 --- a/test/ellip.jl +++ b/test/ellip.jl @@ -48,3 +48,147 @@ @test_throws DomainError ellipe(1.2) end end + +@testset "elliptic edge cases" begin + # Boundary value tests at each polynomial table transition (idx = 0,1,2,...,9) + # These test the boundaries where the binary tree branching changes paths + @testset "ellipk boundary values" begin + # Table boundaries: 0.0, 0.1, 0.2, ..., 0.9, 1.0 + # Reference values from current implementation (regression test baseline) + @test ellipk(0.0) ≈ 1.5707963267948966 rtol=2*eps() # π/2 + @test ellipk(0.1) ≈ 1.6124413487202192 rtol=2*eps() + @test ellipk(0.2) ≈ 1.6596235986105277 rtol=2*eps() + @test ellipk(0.3) ≈ 1.7138894481787912 rtol=2*eps() + @test ellipk(0.4) ≈ 1.7775193714912532 rtol=2*eps() + @test ellipk(0.5) ≈ 1.8540746773013721 rtol=2*eps() + @test ellipk(0.6) ≈ 1.9495677498060260 rtol=2*eps() + @test ellipk(0.7) ≈ 2.0753631352924691 rtol=2*eps() + @test ellipk(0.8) ≈ 2.2572053268208530 rtol=2*eps() + @test ellipk(0.85) ≈ 2.3890164863255676 rtol=2*eps() # Table 10-11 boundary + @test ellipk(0.9) ≈ 2.5780921133481733 rtol=2*eps() + @test ellipk(0.95) ≈ 2.9083372484445515 rtol=2*eps() + @test ellipk(0.99) ≈ 3.6956373629898747 rtol=2*eps() + + # Near-boundary tests using prevfloat/nextfloat + @test ellipk(prevfloat(0.1)) ≈ ellipk(0.1) rtol=1e-10 + @test ellipk(nextfloat(0.1)) ≈ ellipk(0.1) rtol=1e-10 + @test ellipk(prevfloat(0.5)) ≈ ellipk(0.5) rtol=1e-10 + @test ellipk(nextfloat(0.5)) ≈ ellipk(0.5) rtol=1e-10 + @test ellipk(prevfloat(0.9)) ≈ ellipk(0.9) rtol=1e-10 + @test ellipk(nextfloat(0.9)) ≈ ellipk(0.9) rtol=1e-10 + end + + @testset "ellipe boundary values" begin + # Table boundaries: 0.0, 0.1, 0.2, ..., 0.9, 1.0 + # Reference values from current implementation (regression test baseline) + @test ellipe(0.0) ≈ 1.5707963267948966 rtol=2*eps() # π/2 + @test ellipe(0.1) ≈ 1.5307576368977631 rtol=2*eps() + @test ellipe(0.2) ≈ 1.4890350580958527 rtol=2*eps() + @test ellipe(0.3) ≈ 1.4453630644126654 rtol=2*eps() + @test ellipe(0.4) ≈ 1.3993921388974322 rtol=2*eps() + @test ellipe(0.5) ≈ 1.3506438810476753 rtol=2*eps() + @test ellipe(0.6) ≈ 1.2984280350469133 rtol=2*eps() + @test ellipe(0.7) ≈ 1.2416705679458224 rtol=2*eps() + @test ellipe(0.8) ≈ 1.1784899243278386 rtol=2*eps() + @test ellipe(0.85) ≈ 1.1433957918831656 rtol=2*eps() # Table 10-11 boundary + @test ellipe(0.9) ≈ 1.1047747327040722 rtol=2*eps() + @test ellipe(0.95) ≈ 1.0604737277662784 rtol=2*eps() + @test ellipe(0.99) ≈ 1.0159935450252240 rtol=2*eps() + + # Near-boundary tests + @test ellipe(prevfloat(0.1)) ≈ ellipe(0.1) rtol=1e-10 + @test ellipe(nextfloat(0.1)) ≈ ellipe(0.1) rtol=1e-10 + @test ellipe(prevfloat(0.5)) ≈ ellipe(0.5) rtol=1e-10 + @test ellipe(nextfloat(0.5)) ≈ ellipe(0.5) rtol=1e-10 + @test ellipe(prevfloat(0.9)) ≈ ellipe(0.9) rtol=1e-10 + @test ellipe(nextfloat(0.9)) ≈ ellipe(0.9) rtol=1e-10 + end + + @testset "negative m values" begin + # Negative m uses transformation: x = m/(m-1), K(m) = K(x)/sqrt(1-m) + # Reference values from current implementation (regression test baseline) + @test ellipk(-0.1) ≈ 1.5335928197134567 rtol=2*eps() + @test ellipk(-0.5) ≈ 1.4157372084259563 rtol=2*eps() + @test ellipk(-1.0) ≈ 1.3110287771460600 rtol=2*eps() + @test ellipk(-2.0) ≈ 1.1714200841467699 rtol=2*eps() + @test ellipk(-10.0) ≈ 0.7908718902387385 rtol=2*eps() + + @test ellipe(-0.1) ≈ 1.6093590249375296 rtol=2*eps() + @test ellipe(-0.5) ≈ 1.7517712756948178 rtol=2*eps() + @test ellipe(-1.0) ≈ 1.9100988945138559 rtol=2*eps() + @test ellipe(-2.0) ≈ 2.1844381427462012 rtol=2*eps() + @test ellipe(-10.0) ≈ 3.6391380384177681 rtol=2*eps() + end + + @testset "special edge cases" begin + # Values very close to 1.0 + @test ellipk(prevfloat(1.0)) > 10.0 # Should be very large + @test ellipe(prevfloat(1.0)) ≈ 1.0 rtol=1e-6 # Should be close to 1 + + # Very small positive values + @test ellipk(1e-15) ≈ π/2 rtol=1e-10 + @test ellipe(1e-15) ≈ π/2 rtol=1e-10 + + # Consistency: E(m) ≤ π/2 for m ∈ [0,1] + for m in 0.0:0.1:1.0 + @test ellipe(m) <= π/2 + eps() + end + + # Consistency: K(m) ≥ π/2 for m ∈ [0,1) + for m in 0.0:0.1:0.9 + @test ellipk(m) >= π/2 - eps() + end + end +end + + +@testset "ellipke" begin + @testset "bit-wise equality - all branches" begin + # ellipke is implemented independently (not calling ellipk/ellipe internally) + # so bit-wise comparison is a valid oracle test + # Test values covering ALL polynomial branches (idx 0-9) + for m in [0.0, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.825, 0.875, 0.95, 0.99] + K, E = ellipke(m) + @test K === ellipk(m) # bit-wise equality + @test E === ellipe(m) # bit-wise equality + end + end + + @testset "boundary cases m=0 and m=1" begin + @test ellipke(0.0) === (ellipk(0.0), ellipe(0.0)) # (π/2, π/2) + @test ellipke(1.0) === (ellipk(1.0), ellipe(1.0)) # (Inf, 1.0) + end + + @testset "negative m values" begin + # Bit-wise equality with ellipk/ellipe (independent implementations) + for m in [-0.1, -0.5, -1.0, -2.0, -10.0, -1e5] + @test ellipke(m) === (ellipk(m), ellipe(m)) + end + + # Extreme negative values to test x==1.0 branch path + @test ellipke(-1e30) === (ellipk(-1e30), ellipe(-1e30)) + + # Negative infinity special case + @test ellipke(-Inf) === (ellipk(-Inf), ellipe(-Inf)) # (0.0, Inf) + end + + @testset "NaN handling" begin + @test ellipke(NaN) === (NaN, NaN) + end + + @testset "DomainError for m > 1" begin + @test_throws DomainError ellipke(1.1) + @test_throws DomainError ellipke(2.0) + end + + @testset "Missing" begin + @test ellipke(missing) === (missing, missing) + end + + @testset "Float16 and Float32 support" begin + for T in (Float16, Float32) + m = T(0.5) + @test ellipke(m) === (ellipk(m), ellipe(m)) + end + end +end