faster elliptic

[Fil]

Implementation by @toja — https://observablehq.com/d/9c104687d45ef00e

[jrus]

@toja I suspect there may be faster possible implementations than either of these. I haven’t looked very closely but do you know of https://dlmf.nist.gov/19.36 ?

[toja]

@jrus Ah no, I wasn't aware of it. But while I was working on the topic I also implemented Carlson’s elliptic integral of the first kind, R_F(x, y, z) after Numerical Recipes, p. 312 for comparison, but it was slower.

I also came across T. Fukushima, Fast computation of incomplete elliptic integral of first kind by half argument transformation. Numer. Math. 116 (4), pp. 687–719, that is mentioned at https://dlmf.nist.gov/19.36. I want to give it a try at some point, but it requires a lot more code than Bulirsch's algorithm.

[jrus]

I’m pretty confused about the code:

var k_ = (sqrt2 - 1) / (sqrt2 + 1),
    k = sqrt(1 - k_ * k_),
    K = ellipticF(halfPi, k * k),

Why do we take the square root and then immediately square it, never using the square rooted version again? Why aren’t these constants (k_, k*k, K) all precomputed?

[jrus]

Here’s one paper, but I’m unconvinced it actually has the fastest possible methods for this. https://doi.org/10.1080/2151237X.2009.10129277

---

The conformal sphere (with 4 slits from one pole to the equator) -> square projection can probably be most cleanly thought of starting from the stereographic projection. From there it is a map from the quarter of the complex plane with positive real and imaginary parts to the square where the 4 corners come from 0, 1, ∞, i. [With the other 4 quarters of the projection found by copying the original signs.]

Alternately we can think of it as the conformal map from unit disk to to a square with vertices at [1, i, –1, –i ], or from a quarter disk to a triangle with vertices at [0, 1, i], with parts outside the disk filled in by symmetry [invert the magnitude of any points outside the disk, then reflect across the line x + y = 1 at the end]

I should experiment a bit to figure out precisely which of these maps are fastest to approximate. I think we will be able to make a complex rational approximation without needing to call any additional special functions at runtime.

For example, using the technology from
https://arxiv.org/pdf/1911.03696
https://arxiv.org/pdf/1804.08127

Is there any standard in d3-geo-projection for how accurate we want projections to be? i.e. how many digits of precision we are looking for? I expect if we can get away with something accurate to 8 digits, that should be somewhat faster to compute than if we need 13 or 15 digits.

[Fil]

we take the square root and then immediately square it

Good catch!

Is there any standard in d3-geo-projection for how accurate we want projections to be?

I'm not aware of any, and we might want to do a survey. I think it varies a lot between projections, for good and bad reasons. I don't think accuracy is very important in "fancy" projections, which are never going to be used in life or death situations, and speed matters more (esp. for animations/interactions). There is value in having a certain regularity of the derivatives—an approximation which "jumps" might create some Moiré or jagged images.

[toja]

@jrus If you want to revisit the Guyou / Peirce quincuncial you might want to take a look here:
https://observablehq.com/d/fc1d2db223f7ed90

I've started to implement them following Snyder's Album of map projections. Since these version use a static k² = 0.5 it's possible to use a Chebyshev series to compute the elliptic integral. But there're some problems at the poles and also the inverse needs a lot optimization

[jrus]

I haven’t read this and am headed to sleep, but it looks promising:

http://www.acsel-lab.com/arithmetic/arith22/papers/8664a050.pdf
https://scholar.google.com/scholar?cluster=17190421556177330303

[jrus]

I have been trying to make rational-function approximations to the quincuncial projection using the AAA algorithm (https://people.maths.ox.ac.uk/trefethen/AAAfinal.pdf), with pretty good success. If I reduce the domain to half an octant of the (stereographically projected) globe, and then approximate a map returning (1 - p(z))^2, where p is the quincuncial projection, a degree 10 rational approximation is accurate to about 14.5 digits throughout the domain. It takes a degree 9 approximant for the inverse map accurate to about 14.5 digits. [I think these could be made accurate to another digit or two, but would have to be constructed using higher-precision arithmetic, which I can’t easily do from Matlab.]

So these functions are going to require overall a bit of logic and possibly a couple conditional branches for domain reduction, then about 12–15 divisions, and 3 square roots (to perform the single complex square root; if this is a bottleneck it might be possible to do something more optimized than the naive implementation), plus maybe 40 additions and multiplications.

I think they should turn out quite snappy, but I still have to implement some barycentric complex rational function evaluation code in Javascript before I can show you something concrete.

[jrus]

Basically, the idea is if we map between these two shapes, there are no square root singularities nearby, making it much easier for a rational function to handle. The black dots are the nodes used for the barycentric formulation of the approximation, and the red dots are the poles of the approximation.

https://i.imgur.com/4cicDa9.png

(I rotated the output shape by 180° in this picture so it’s easier to see which vertex goes to which.)

The maximum errors around the boundary are ~3e-15 in the forward direction and ~2e-15 in the inverse direction; should be a comparable level of error in the interior.

[jrus]

Still working on this (e.g. I'm going to switch to nodes that are slightly outside the reduced domain so that there’s no accidental chance of a NaN from evaluating Infinity / Infinity), but feel free to play with it:

https://observablehq.com/d/89028e12f98e34aa

I was a bit optimistic with my estimate of the number of multiplications/additions. More like 300 of them, plus the 12–14 divisions and 3 square roots, per projected point. (If Javascript could do SIMD, we could use all the SIMD lanes and cut that down by a factor of 4 or whatever, since most of the arithmetic is independent and can be done in parallel.) Still should be pretty quick though I expect.

I’ll try to test the speed out when I get a chance. Edit: it can handle about 9 million projected points per second on my desktop computer.

[jrus]

@toja here’s a fork of your notebook where I only replaced the forward direction but not the inverse

https://observablehq.com/d/47ceb9539492ff76

also ping @Fil

[jrus]

@toja @Fil okay, here’s the inverse as well https://observablehq.com/@jrus/fast-peirce-quincuncial

[Fil]

I'm testing @jrus's implementation in d3-geo-projection.

The (trivial) difference is that the current projection has a default rotate of [-90, -90, 45], which creates a weird default aspect (shown in the difference (black&white) map below). For a nicer aspect, one has to configure it as, e.g. d3.geoPeirceQuincuncial().rotate([-90 + 25, -90, 45]).

With Jacob's implementation the aspect given by the suggested rotate([25, 0, 0]) is perfect, as in this yellow map:

Other than that, the difference shown below (with the new implementation at [-90, 0, 0] to match the original aspect) is minimal—I suppose it comes from the new projection being more precise :)

Notebook: https://observablehq.com/d/6bbd20a7462e68f5

One question if we want to use the new implementation, is to know if we can replace the Guyou projection too, or if we're going to have a diverging code base?

[Fil]

Added fast PeirceQuincuncial in #199

[jrus]

Can definitely use this to implement Guyou as well. I might not have the bandwidth to look at it in the next few days though.