Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add derivatives for ellipk and ellipe #370

Merged
merged 3 commits into from
Jan 26, 2022
Merged

Add derivatives for ellipk and ellipe #370

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8be3895
Add derivatives for `ellipk` and `ellipe`
devmotion Dec 16, 2021
f92faf9
Bump version
devmotion Dec 16, 2021
615d577
Fix derivatives at 0
devmotion Dec 16, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "SpecialFunctions"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.0.0"
version = "2.1.0"

[deps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
Expand Down
4 changes: 4 additions & 0 deletions src/chainrules.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -194,6 +194,10 @@ ChainRulesCore.@scalar_rule(expinti(x), exp(x) / x)
ChainRulesCore.@scalar_rule(sinint(x), sinc(invπ * x))
ChainRulesCore.@scalar_rule(cosint(x), cos(x) / x)

# elliptic integrals
ChainRulesCore.@scalar_rule(ellipk(m), (ellipe(m) / (1 - m) - Ω) / (2 * m))
ChainRulesCore.@scalar_rule(ellipe(m), (Ω - ellipk(m)) / (2 * m))

# non-holomorphic functions
function ChainRulesCore.frule((_, Δν, Δx), ::typeof(besselix), ν::Number, x::Number)
# primal
Expand Down
18 changes: 9 additions & 9 deletions src/ellip.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,9 @@ For ``m<0``, followed by
> <https://www.researchgate.net/publication/267330394>
As suggested in this paper, the domain is restricted to ``(-\infty,1]``.
"""
function ellipk(m::Float64)
ellipk(m::Real) = _ellipk(float(m))

function _ellipk(m::Float64)
flag_is_m_neg = false
if m < 0.0
x = m / (m-1) #dealing with negative args
Expand Down Expand Up @@ -214,7 +216,9 @@ For ``m<0``, followed by
> <https://www.researchgate.net/publication/267330394>
As suggested in this paper, the domain is restricted to ``(-\infty,1]``.
"""
function ellipe(m::Float64)
ellipe(m::Real) = _ellipe(float(m))

function _ellipe(m::Float64)
flag_is_m_neg = false
if m < 0.0
x = m / (m-1) #dealing with negative args
Expand Down Expand Up @@ -346,11 +350,7 @@ function ellipe(m::Float64)
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
# Support for Float32 and Float16
for internalf in (:_ellipk, :_ellipe), T in (:Float16, :Float32)
@eval $internalf(x::$T) = $T($internalf(Float64(x)))
end
6 changes: 6 additions & 0 deletions test/chainrules.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,12 @@
test_scalar(airyaiprimex, x)
end
end

# finite differencing fails for x = 0
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is wrong, I misinterpreted the ChainRules tests. I always mix up what "expected" and "actual" refer to.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added the limit of the derivative at 0 as a special case.

if x isa Real && x < 1 && !iszero(x)
test_scalar(ellipk, x)
test_scalar(ellipe, x)
end
end
end

Expand Down