Jacobi elliptic functions and complete elliptic integral of the first kind #79

ettersi · 2018-03-05T18:20:04Z

Implements the Jacobi elliptic functions and the complete elliptic integral of the first kind.

Both sets of functions work for real and complex arguments as well as all FloatXX and BigFloat types. The test sets involve comparison with Python's mpmath library (http://mpmath.org/) as well as checks that the results are within a few eps() of the exact result.

simonbyrne · 2018-03-05T19:39:54Z

Thanks! Is this based on the mpmath code? If so, that's fine as it's a compatible licence, but we need to acknowledge the original copyright.

simonbyrne · 2018-03-05T19:41:00Z

src/SpecialFunctions.jl

@@ -71,7 +71,12 @@ end
 export sinint,
       cosint

+export ellipj,
+       jss,jsc,jsd,jsn,jcs,jcc,jcd,jcn,jds,jdc,jdd,jdn,jns,jnc,jnd,jnn,


What are all these exports?

Ah, I see, they're generated by macros. The names are a bit cryptic: are these standard? Maybe prefix them (elliptic_jss) or stick them in a module (Elliptic.jss)?

There are 16 Jacobi elliptic functions, namely all pairs of the letters scdn. Since cd is already taken, I decided to put a j (for Jacobi) in front of the function to avoid name collision. Elliptic.jl solved this issue by putting these functions into a module Jacobi, but I'd say jsn(u,m) is more convenient than Jacobi.sn(u,m). Might be worth discussing, though.

ah ok, how about jacobisn? This matches nicely with besselk/besselj, etc.

jacobisn possibly makes it hard to distinguish the functions, and it wouldn't really save much typing compared to Jacobi.sn. It might be worth to get input from someone who uses these functions regularly, though.

Looking at other packages:

Mathematica: JacobiSN

Matlab & SciPy provides them all at once via ellipj function.

GSL provides them all via gsl_sf_elljac_e

We might also want to rename the K and iK functions. It's just too much of a hassle not to be able to use K for anything else. For reference:

Mathematica: EllipticK

Matlab: ellipticK

Scipy: ellipk

GSL: gsl_sf_ellint_Kcomp
I'd vote for ellipK, as the K is usually uppercase in mathematical notation. Definitely not ellipticK unless we rename ellipj.

Also, I am a bit unsure about providing iK and how it should be named. If you look at the code, it definitely makes sense to provide iK because K(m) is really just iK(1-m) and so K(1-m) = iK(1-(1-m)) which is just stupid. But AFAIK no other software package provides this, which means I might be doing something wrong. And the commonly used mathematical notation for K(1-m) seems to be K', but obviously we can't use that so we need to think of something else.

I'd still prefer to rename the functions to longer names.

src/elliptic.jl

ettersi · 2018-03-05T19:46:21Z

I skimmed through the mpmath codes (and a few others), but at most I took ideas, no code.

src/elliptic.jl

simonbyrne · 2018-03-06T00:00:49Z

src/elliptic.jl

+#----------------
+# Pick algorithm
+
+Base.@pure puresqrt(x) = sqrt(x)


Isn't sqrt already pure?

Actually, it isn't for BigFloat, but that is intentional, since we can change precision via setprecision. Using @pure here is incorrect, and will cause problems when computing at different precisions.

It also isn't for the FloatXX types:

julia> foo(T) = Val{sqrt(eps(T))}() foo (generic function with 1 method) julia> @code_warntype foo(Float64) Variables: #self# <optimized out> T <optimized out> Body: begin return ((Core.apply_type)(Main.Val, (Base.Math.sqrt_llvm)(2.220446049250313e-16)::Float64)::Type{Val{_}} where _)()::Val{_} where _ end::Val{_} where _ julia> foo(T) = Val{puresqrt(eps(T))}() foo (generic function with 1 method) julia> @code_warntype foo(Float64) Variables: #self# <optimized out> T <optimized out> Body: begin return $(QuoteNode(Val{1.4901161193847656e-8}())) end::Val{1.4901161193847656e-8}

But I agree, I should have done

puresqrt(x) = sqrt(x) Base.@pure puresqrt(x::Union{Float16,Float32,Float64}) = sqrt(x)

With this, I get

julia> foo(T) = Val{puresqrt(eps(T))}() foo (generic function with 1 method) julia> @code_warntype foo(Float64) Variables: #self# <optimized out> T <optimized out> Body: begin return $(QuoteNode(Val{1.4901161193847656e-8}())) end::Val{1.4901161193847656e-8} julia> @code_warntype foo(BigFloat) Variables: #self# <optimized out> T::Type{BigFloat} Body: begin SSAValue(0) = $(Expr(:invoke, MethodInstance for eps(::Type{BigFloat}), :(Main.eps), :(T))) return ((Core.apply_type)(Main.Val, $(Expr(:invoke, MethodInstance for sqrt(::BigFloat), :(Main.sqrt), SSAValue(0))))::Type{Val{_}} where _)()::Val{_} where _ end::Val{_} where _

Just out of curiosity, is there any fundamental obstacle to do this in Base, and possible for all numerical functions?

Actually, maybe it would be better to just do

Base.@pure puresqrt(x) = Float64(sqrt(x))

Then puresqrt is pure, and since the nsteps computations are done in Float64 anyway, I don't lose anything in terms of generality.

simonbyrne · 2018-03-06T00:04:25Z

Overall, this looks pretty good to me. It would be nice to have a few more comments, particularly as to how the generated functions work and why they're needed.

I'm sure @stevengj will have some thoughts as well.

simonbyrne · 2018-03-06T00:18:50Z

Actually, I am still slightly hesitant about all the generated functions. It would be nice if we could reduce their use a bit.

I think the ellipj_* functions could be written as simple iterators, does that impact performance?
shrinkm is a bit more difficult, but currently it is unstable for arbitrary precision types like BigFloat (well, it should be, see my comment on puresqrt). Perhaps you could have it return tuples for fixed-precision types (Float64, Float32, etc.) and arrays for BigFloat?

ettersi · 2018-03-06T10:43:18Z

The main reason for the generated functions is to avoid heap allocation. The structure of the algorithms is to first compute k_i for i = 1:nsteps and then use them in the order reverse(1:nsteps). Without generated functions, I would have to either

allocate an array on the heap, or
recompute them from the last k_i.
I don't like the heap allocation because it's performance impact can be a bit unpredictable, and recomputing them is probably slower and could affect numerical accuracy. Because of the "up and down" structure of the algorithm, I am not sure what would be gained if I implemented the ellipj_shrinkm and ellipj_growm functions as iterators.

Regarding the instability of BigFloats: as far as I understand, the only difference is that eps(BigFloat) is no longer pure, hence the number of Landen transformation steps has to be computed at runtime and the calls in ellipj_dispatch have to be resolved at runtime rather than compile-time. That seems like a small price to pay, though, especially since the alternative would be to do heap allocation of the k_i which would probably have a similar performance penalty. Making the shrinkm function switch between returning tuples or arrays would either lead to code duplication or at least complicate the rather simple code we have at the moment. Both alternatives seem undesirable for maintenance. If you are worried about the number of methods that will be compiled because of this generated functions business, I don't think that will be much of a concern. Most users probably don't bother changing the precision of BigFloats, and even if they do, the Landen transformation converges quadratically so even nsteps(1e-300, Complex) is just 8.

There is a slight hick-up in that currently we compute the number of Landen steps twice for BigFloats falling into the last two branches of ellipj_dispatch. I'll fix that.

simonbyrne · 2018-03-06T18:07:52Z

What I meant was

function ellipj_growm(sn,cn,dn, k)
    # [1, Sec 16.12] / https://dlmf.nist.gov/22.7.i
    for kk in reverse(k)
        sn,cn,dn = (1+kk)*sn/(1+kk*sn^2), cn*dn/(1+kk*sn^2), (1-kk*sn^2)/(1+kk*sn^2)
                     # ^ Use [1, 16.9.1]. Idea taken from [2]
    end
    return sn,cn,dn
end

Basically, generated functions should be considered a weapon of last resort: they are powerful, but do add considerable complexity, both for the compiler (i.e. they make precompilation more difficult) and for humans reading the code.

As for BigFloats: these are already heap allocated, so array vs tuple for them are a bit of a wash (though, again, benchmarking would be useful).

ettersi · 2018-03-06T21:39:48Z

Right, I agree I overused @generated for these two functions. Sorry about that, fixed it.

I also agree that there probably isn't much point in avoiding heap allocations for BigFloats, but the types we want to optimise for are Float64 and Float32 and there it probably matters. And I'd still say there is not much point in having two versions of the same code because of this.

src/elliptic.jl

ettersi · 2018-03-20T20:31:34Z

The only remaining issues with this PR seem to be about the names of the functions. Currently, we have

ellipj following Matlab and SciPy
jpq for p,q in (s,c,d,n). This is my own invention.
ellipK compared to ellipticK (Matlab) and ellipk (SciPy).

Some more feedback regarding this would be helpful.

There further was some discussion regarding my use of generated functions and the @pure annotation. The current usage seems both correct and reasonable to me, but please let me know if there are any doubts remaining.

MasonProtter · 2018-03-29T02:07:41Z

Is there a timeline for this to merge?

MasonProtter · 2018-03-29T02:36:22Z

I couldn't find an appropriate place in the code to comment on this so I'll do it here if thats okay. Currently the algorithm doesn't know to use the analytic continuation for real valued arguments > 1:

julia> ellipK(2.0)
ERROR: DomainError:
Stacktrace:
 [1] ellipiK_agm(::Float64) at /Users/mason/.julia/v0.6/SpecialFunctions/src/elliptic.jl:207
 [2] ellipK(::Float64) at /Users/mason/.julia/v0.6/SpecialFunctions/src/elliptic.jl:214

whereas giving a complex argument will cause it to use the correct analytic continuation

julia> ellipK(2.0+0im)
1.3110287771460596 - 1.3110287771460596im

Is this intended behaviour?

ettersi · 2018-03-29T07:49:06Z

This is intentional and consistent with Base.sqrt:

julia> sqrt(-1)
ERROR: DomainError:
sqrt will only return a complex result if called with a complex argument. Try sqrt(complex(x)).
Stacktrace:
 [1] sqrt(::Int64) at ./math.jl:434

We could and possible should print a similar error message for ellipK, but I couldn't quite figure out how to do this on 0.6. Advice from a more experienced Julia user would be appreciated. I guess in 0.7 I should just call Base.Math.throw_complex_domainerror()?

ettersi · 2018-06-28T03:15:07Z

Repeating Mason Protter's question: Is there a timeline for this to merge?

Datseris · 2018-09-08T20:51:09Z

If you ever do 2nd kind elliptic integrals, please tag me, so I can benchmark and compare with https://github.com/nolta/Elliptic.jl , which is my current dependency.

(assuming of course they ever get merged)

ettersi · 2019-02-21T16:48:58Z

Once again, is there any timeline on getting this merged? If not, why?

jebej · 2020-07-14T20:18:43Z

bump

simonbyrne · 2020-07-14T21:59:39Z

src/elliptic.jl

+Base.@pure function nsteps(m,ε)
+    i = 0
+    while abs(m) > ε
+        m = descstep(m)^2
+        i += 1
+    end
+    return i
+end
+Base.@pure nsteps(ε,::Type{<:Real}) = nsteps(0.5,ε) # Guarantees convergence in [-1,0.5]
+Base.@pure nsteps(ε,::Type{<:Complex}) = nsteps(0.5+sqrt(3)/2im,ε) # This is heuristic.


I still don't think the @pures here are necessary (unless you have benchmarks suggesting otherwise

These functions aren’t pure anyways...

Right, this needs to be looked into.

simonbyrne · 2020-07-14T21:59:46Z

src/elliptic.jl

+    # For all FloatXX types, this can be done at compile time, while for
+    # BigFloat this has to be done at runtime.
+    T = promote_type(typeof(u),typeof(m))
+    ε = puresqrt(Float64(eps(real(typeof(m)))))


Suggested change

ε = puresqrt(Float64(eps(real(typeof(m)))))

ε = sqrt(Float64(eps(real(typeof(m)))))

should be sufficient in most cases.

ettersi · 2020-07-16T09:39:05Z

Here's a list of remaining issues with this PR:

Decide on function names (see Jacobi elliptic functions and complete elliptic integral of the first kind #79 (review)). Fixable in 15 minutes once the decision is made.
Bring this PR up to date with master, and update to Julia 1. Should be doable within an hour or so.
Sort out the @pure (see Jacobi elliptic functions and complete elliptic integral of the first kind #79 (review)). Might need a couple of hours to do proper benchmarking and coming up with a fix.
Decide whether we want to keep this implementation of ellipk or the one from Added complete elliptic functions for 1st and 2nd kind #135.

Regarding the last point: the advantage of this implementation is that it supports complex arguments and arbitrary number types. On the other hand, #135 provides a uniform implementation for both ellipk and ellipe. Also, we should check whether there's a significant performance difference between the two codes.

Sorting this out might take a couple of days to read through the literature, benchmark, maybe implement ellipe in a similar style as the one in this PR, etc.

If anyone is willing to tackle these issues, please feel free to do so. I'd be happy to give a hand, but I won't invest any further time unless there's a clear path to how this PR might finally get merged.

stevengj · 2020-12-15T15:58:21Z

Sorry I haven't chimed in here, but @simonbyrne's comments all seem pretty reasonable to me and I had been waiting fo them to be addressed before I looked at this PR closely.

ptiede · 2022-01-12T19:37:25Z

src/elliptic.jl

+Base.@pure function nsteps(m,ε)
+    i = 0
+    while abs(m) > ε
+        m = descstep(m)^2
+        i += 1
+    end
+    return i
+end
+Base.@pure nsteps(ε,::Type{<:Real}) = nsteps(0.5,ε) # Guarantees convergence in [-1,0.5]
+Base.@pure nsteps(ε,::Type{<:Complex}) = nsteps(0.5+sqrt(3)/2im,ε) # This is heuristic.
+function ellipj_nsteps(u,m)
+    # Compute the number of Landen steps required to reach machine precision.
+    # For all FloatXX types, this can be done at compile time, while for
+    # BigFloat this has to be done at runtime.
+    T = promote_type(typeof(u),typeof(m))
+    ε = puresqrt(Float64(eps(real(typeof(m)))))
+    N = nsteps(ε,typeof(m))
+    return ellipj_dispatch(u,m,Val{N}())::NTuple{3,T}
+end


Hi. Sorry if this isn't in the correct place. I was hoping to this would be finished. From my reading of this pull request the major stumbling block in this pull request are these lines that use a lot of @pure. When I tried to test this I found that the @pure where needed to prevent dynamic dispatch. I tried to modify the code a bit to prevent this and wound up at the following

Suggested change

Base.@pure function nsteps(m,ε)

i = 0

while abs(m) > ε

m = descstep(m)^2

i += 1

end

return i

end

Base.@pure nsteps(ε,::Type{<:Real}) = nsteps(0.5,ε) # Guarantees convergence in [-1,0.5]

Base.@pure nsteps(ε,::Type{<:Complex}) = nsteps(0.5+sqrt(3)/2im,ε) # This is heuristic.

function ellipj_nsteps(u,m)

# Compute the number of Landen steps required to reach machine precision.

# For all FloatXX types, this can be done at compile time, while for

# BigFloat this has to be done at runtime.

T = promote_type(typeof(u),typeof(m))

ε = puresqrt(Float64(eps(real(typeof(m)))))

N = nsteps(ε,typeof(m))

return ellipj_dispatch(u,m,Val{N}())::NTuple{3,T}

end

Base.@pure function nsteps(M)

m = _convfac(M)

ε = _vareps(M)

i = 0

while abs(m) > ε

m = descstep(m)^2

i += 1

end

return i

end

@inline _convfac(::Type{<:Real}) = 0.5

@inline _convfac(::Type{<:Complex}) = 0.5 + sqrt(3)/2im

@inline _vareps(::Type{<:Complex{T}}) where {T} = Float64(sqrt(eps(T)))

@inline _vareps(T::Type{<:Real}) = Float64(sqrt(eps(T)))

function ellipj_nsteps(u,m)

# Compute the number of Landen steps required to reach machine precision.

# For all FloatXX types, this can be done at compile time, while for

# BigFloat this has to be done at runtime.

T = promote_type(typeof(u),typeof(m))

N = nsteps(T)

return ellipj_dispatch(u,m,Val{N}())::NTuple{3,T}

end

This seems to work in my testing. For versions < 1.7 I still needed the @pure to prevent dynamic dispatch. I am not sure if this usage of @pure is alright, but it seems to be more inline with what I read in the docs. The good news is that for 1.7 everything infers correctly so the @pure is not needed.

What's the chances that this gets implemented soonish? I can try to take over this pull request @ettersi if any additional work is needed. I just need this implementation of the jacobi functions for some research I am doing.

Simon Etter added 8 commits March 1, 2018 14:14

Add ellipj implementation

396fc72

Pre-cleanup

863bd9b

Post cleanup

68fff8a

Finish Jacobi elliptic functions

36f2262

Add jpq functions

8f22ff7

Include in SpecialFunctions and add tests

3825f65

Update

2660d16

Update documentation

de7d83b

This was referenced Mar 5, 2018

Support for extended precision (BigFloats)? nolta/Elliptic.jl#11

Open

elliptic integrals with complex argument nolta/Elliptic.jl#3

Closed

Fix 0.7 deprecations and erros.

6f0d23d

simonbyrne reviewed Mar 5, 2018

View reviewed changes

src/elliptic.jl Show resolved Hide resolved

Fix test error on 32bit systems (fingers crossed)

2d406dc

simonbyrne reviewed Mar 5, 2018

View reviewed changes

src/elliptic.jl Outdated Show resolved Hide resolved

Simon Etter added 3 commits March 5, 2018 22:28

Factor trigonometric and hyperbolic functions into variables

9f16249

Align equations

dbc4bd4

Add references to DLMF

3b34edd

simonbyrne reviewed Mar 6, 2018

View reviewed changes

Make sure nsteps gets called only once even for BigFloats.

47d4a55

Remove some @generated.

0eee6b2

simonbyrne reviewed Mar 6, 2018

View reviewed changes

src/elliptic.jl Outdated Show resolved Hide resolved

simonbyrne reviewed Mar 6, 2018

View reviewed changes

src/elliptic.jl Show resolved Hide resolved

Simon Etter added 3 commits March 6, 2018 22:46

Optimise evaluation of tanh(u) and sech(u)

f4f00b9

Rename (i)K -> ellip(i)K

613fd93

Fix puresqrt()

37e5999

simonbyrne assigned stevengj Mar 23, 2018

simonbyrne reviewed Jul 14, 2020

View reviewed changes

ptiede reviewed Jan 12, 2022

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Jacobi elliptic functions and complete elliptic integral of the first kind #79

Jacobi elliptic functions and complete elliptic integral of the first kind #79

ettersi commented Mar 5, 2018

simonbyrne commented Mar 5, 2018

simonbyrne Mar 5, 2018

simonbyrne Mar 5, 2018

ettersi Mar 5, 2018

simonbyrne Mar 5, 2018

ettersi Mar 5, 2018

simonbyrne Mar 5, 2018

ettersi Mar 6, 2018

simonbyrne Jul 14, 2020

ettersi commented Mar 5, 2018

simonbyrne Mar 6, 2018

simonbyrne Mar 6, 2018 •

edited

Loading

ettersi Mar 6, 2018

ettersi Mar 6, 2018

simonbyrne commented Mar 6, 2018

simonbyrne commented Mar 6, 2018

ettersi commented Mar 6, 2018

simonbyrne commented Mar 6, 2018

ettersi commented Mar 6, 2018

ettersi commented Mar 20, 2018

MasonProtter commented Mar 29, 2018

MasonProtter commented Mar 29, 2018

ettersi commented Mar 29, 2018 •

edited

Loading

ettersi commented Jun 28, 2018

Datseris commented Sep 8, 2018 •

edited

Loading

ettersi commented Feb 21, 2019

jebej commented Jul 14, 2020

simonbyrne Jul 14, 2020

MasonProtter Jul 14, 2020 •

edited

Loading

ettersi Jul 16, 2020

simonbyrne Jul 14, 2020

ettersi commented Jul 16, 2020 •

edited

Loading

stevengj commented Dec 15, 2020

ptiede Jan 12, 2022

	ε = puresqrt(Float64(eps(real(typeof(m)))))
	ε = sqrt(Float64(eps(real(typeof(m)))))

Jacobi elliptic functions and complete elliptic integral of the first kind #79

Are you sure you want to change the base?

Jacobi elliptic functions and complete elliptic integral of the first kind #79

Conversation

ettersi commented Mar 5, 2018

simonbyrne commented Mar 5, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ettersi commented Mar 5, 2018

Choose a reason for hiding this comment

simonbyrne Mar 6, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

simonbyrne commented Mar 6, 2018

simonbyrne commented Mar 6, 2018

ettersi commented Mar 6, 2018

simonbyrne commented Mar 6, 2018

ettersi commented Mar 6, 2018

ettersi commented Mar 20, 2018

MasonProtter commented Mar 29, 2018

MasonProtter commented Mar 29, 2018

ettersi commented Mar 29, 2018 • edited Loading

ettersi commented Jun 28, 2018

Datseris commented Sep 8, 2018 • edited Loading

ettersi commented Feb 21, 2019

jebej commented Jul 14, 2020

Choose a reason for hiding this comment

MasonProtter Jul 14, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ettersi commented Jul 16, 2020 • edited Loading

stevengj commented Dec 15, 2020

Choose a reason for hiding this comment

simonbyrne Mar 6, 2018 •

edited

Loading

ettersi commented Mar 29, 2018 •

edited

Loading

Datseris commented Sep 8, 2018 •

edited

Loading

MasonProtter Jul 14, 2020 •

edited

Loading

ettersi commented Jul 16, 2020 •

edited

Loading