Improve numerical precision of inverse spherical mercator

[kbevers]

Supersedes #931

[cffk]

Much better would be to replace

pi/2 - 2*atan(exp(x))

by

-atan(sinh(x))

This is accurate everywhere, is demonstrably odd, and doesn't require branches or C99 functions.

[kbevers]

This is accurate everywhere, is demonstrably odd, and doesn't require branches or C99 functions.

I knew you would have a better way to do this. Updated the code and it is indeed much simpler. Thanks for the suggestion!

[rouault]

lots of negation there. I guess atan(sinh(xy.y / P->k0)) would do ?

[cffk]

Yes, the two negations cancel.

[kbevers]

You guessed right! Nice catch.

[rouault]

Perhaps there should be a new test point in gie tests to demonstrate the increased accuracy ?

[kbevers]

Perhaps there should be a new test point in gie tests to demonstrate the increased accuracy ?

Yeah. The original PR didn't have on and for some reason I can't come up with one that demonstrates that this approach is better. Perhaps it isn't more precise than the original code?

[rouault]

The difference between both approaches only show up at extremely small absolute values

>>> math.pi / 2 - 2 * math.atan(math.exp(-1e-15))
8.881784197001252e-16
>>> math.atan(math.sinh(1e-15))
1e-15

At 1e-10 and above, the relative error is below 1e-7
At 1e-5 and above, the relative error is below 1e-10
Perhaps the advantage of the new formula would be that it involves less floating point operations. Although it is not clear if sinh() implementations are equivalently optimized as the common exp()

[kbevers]

Yeah, and my tests has been using extremely small inputs. From the below I see a change, but not exactly an improvement. I would expect to the input to be repeated in the output in this example:

# This PR
$ echo 0 1e-15 | ./bin/proj -I -f %.9e +proj=merc +R=1
0.000000000e+00 5.729577951e-14

# PROJ 5.1
$ echo 0 1e-15 | proj -I -f %.9e +proj=merc +R=1
0.000000000e+00 5.088887490e-14

[kbevers]

Disregard my conclusion from the previous post. With the changes in this PR the precision is in fact improved. I forgot to account for the RAD2DEG conversion that proj does to the output. The output matches the input when that is taken into account.