<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Dear Even, Charles, Thomas, All,</p>
<p>Please find below a couple revised implementations of the
authalic ==> geodetic conversion using Horner's method and
Clenshaw summation algorithm, both sharing this table of
coefficients from A20:<br>
</p>
<p>
<blockquote type="cite"><font size="4" face="monospace">#define
AUTH_ORDER 6<br>
</font><br>
<font size="4" face="monospace">static const double Cphimu[21] =
// Cφξ (A20) - coefficients to convert authalic latitude to
geodetic latitude <br>
{ <br>
4 / 3.0, 4 / 45.0, -16/35.0, -2582 /14175.0, 60136
/467775.0, 28112932/ 212837625.0, <br>
46 / 45.0, 152/945.0, -11966 /14175.0, -21016 /
51975.0, 251310128/ 638512875.0, <br>
3044/2835.0, 3802 /14175.0, -94388 /
66825.0, -8797648/ 10945935.0, <br>
6059 / 4725.0, 41072 /
93555.0, -1472637812/ 638512875.0, <br>
768272
/467775.0, -455935736/ 638512875.0, <br>
4210684958/1915538625.0 <br>
}; <br>
</font></blockquote>
</p>
<p>This first one is using the existing functions from <i>mlfn.cpp</i>
(untouched other than possibly different formatting here):<br>
</p>
<p>
<blockquote type="cite"><font size="4" face="monospace">//
Evaluate sum(p[i] * x^i, i, 0, N) via Horner's method (p is of
length N+1) <br>
static inline double polyval(double x, const double p[], int
N) <br>
{ <br>
double y = N < 0 ? 0 : p[N]; <br>
while(N > 0) <br>
y = y * x + p[--N]; <br>
return y; <br>
} <br>
<br>
// Evaluate y = sum(c[k] * sin((2*k+2) * zeta), k, 0, K-1) <br>
static inline double clenshaw(double szeta, double czeta,
const double c[], int K) <br>
{ <br>
// Approx operation count = (K + 5) mult and (2 * K + 2)
add <br>
double u0 = 0, u1 = 0; // accumulators for sum <br>
double X = 2 * (czeta - szeta) * (czeta + szeta); // 2 *
cos(2*zeta) <br>
while(K > 0) <br>
{ <br>
double t = X * u0 - u1 + c[--K]; <br>
u1 = u0; <br>
u0 = t; <br>
} <br>
return 2 * szeta * czeta * u0; // sin(2*zeta) * u0 <br>
} <br>
<br>
// <a class="moz-txt-link-freetext" href="https://arxiv.org/pdf/2212.05818">https://arxiv.org/pdf/2212.05818</a> <br>
// ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1) -- (20) <br>
void pj_authset(double a, double b, double cp[AUTH_ORDER]) <br>
{ <br>
double n = (a - b) / (a + b); // Third flattening <br>
double d = n; <br>
int l, o; <br>
<br>
for(l = 0, o = 0; l < AUTH_ORDER; l++) <br>
{ <br>
int m = AUTH_ORDER - l - 1; <br>
<br>
cp[l] = d * polyval(n, Cphimu + o, m); <br>
d *= n; <br>
o += m + 1; <br>
} <br>
} <br>
<br>
double pj_auth2geodlat(const double * cp, double phi) <br>
{ <br>
return phi + clenshaw(sin(phi), cos(phi), cp, AUTH_ORDER);
<br>
} </font><br>
</blockquote>
</p>
<p>For this second implementation, I unrolled the loops to get rid
of the iterations (and associated counter incrementations) and
conditionals, which if the compiler is not optimizing out, could
potentially introduce some <a
href="https://en.algorithmica.org/hpc/pipelining/branching/">branching
costs</a>.<br>
This unrolled version remains quite compact (at least in this
particular formatting which the pre-commit hook will certainly
massacre). The sequence of operations is exactly the same, and
I've tested that the two are equivalent, and also equivalent with
the two earlier implementations that I shared which were not using
Horner and Clenshaw, and also equivalent to 8 decimals to the
existing <i>pj_authlat()</i> function in PROJ.<br>
</p>
<p>
<blockquote type="cite"><font size="4" face="monospace">//
<a class="moz-txt-link-freetext" href="https://arxiv.org/pdf/2212.05818">https://arxiv.org/pdf/2212.05818</a> <br>
// ∆η(ζ) = S^(L)(ζ) · Cηζ · P^(M) (n) + O(n^L+1) -- (20) <br>
void pj_authset(double a, double b, double cp[AUTH_ORDER]) <br>
{ <br>
// Precomputing coefficient based on Horner's method <br>
double n = (a - b) / (a + b); // Third flattening <br>
const double * C = Cphimu; <br>
double d = n; <br>
<br>
cp[0] = (((((C[ 5] * n + C[ 4]) * n + C[ 3]) * n + C[ 2]) *
n + C[ 1]) * n + C[ 0]) * d, d *= n; <br>
cp[1] = (((( C[10] * n + C[ 9]) * n + C[ 8]) *
n + C[ 7]) * n + C[ 6]) * d, d *= n; <br>
cp[2] = ((( C[14] * n + C[13]) *
n + C[12]) * n + C[11]) * d, d *= n; <br>
cp[3] = (( C[17] *
n + C[16]) * n + C[15]) * d, d *= n; <br>
cp[4] =
( C[19]
* n + C[18]) * d, d *= n; <br>
cp[5]
=
C[20] * d; <br>
} <br>
<br>
double pj_auth2geodlat(const double * cp, double phi) <br>
{ <br>
// Using Clenshaw summation algorithm (order 6) <br>
double szeta = sin(phi), czeta = cos(phi); <br>
// Approx operation count = (K + 5) mult and (2 * K + 2)
add <br>
double X = 2 * (czeta - szeta) * (czeta + szeta); // 2 *
cos(2*zeta) <br>
double u0 = 0, u1 = 0; // accumulators for sum <br>
double t; <br>
t = X * u0 - u1 + cp[5], u1 = u0, u0 = t; <br>
t = X * u0 - u1 + cp[4], u1 = u0, u0 = t; <br>
t = X * u0 - u1 + cp[3], u1 = u0, u0 = t; <br>
t = X * u0 - u1 + cp[2], u1 = u0, u0 = t; <br>
t = X * u0 - u1 + cp[1], u1 = u0, u0 = t; <br>
t = X * u0 - u1 + cp[0]; <br>
return phi + /* sin(2*zeta) * u0 */ 2 * szeta * czeta * t;
<br>
}</font><br>
</blockquote>
Note that the output of these two versions of <i>pj_authset()</i>
(the 6 constants precomputed from the authalic ==> geodetic A20
conversion matrix and the ellipsoid's third flattening) is exactly
the same as the previous version not using Horner's method, and I
believe also the same as the current output of <i>pj_autset() </i>except
that it currently uses only 3 constants for order 3 rather than 6
for order 6.</p>
<p>With both of these versions, we're down to only one <i>sin()</i>
and one <i>cos()</i> call, as per Even's suggestion, so I imagine
that the Clenshaw algorithm does take advantage of that
trigonometric identity trick.<br>
<br>
If we go with the separate <i>polyval()</i> and <i>clenshaw()</i>
functions, then I suggest we move these functions to a header file
so that we can share them between <i>mlfn.cpp</i> and <i>auth.cpp</i>
while allowing the compiler to hopefully efficiently inline them,
and also hopefully optimize the code close to or equivalent to the
unrolled version (we could always compare the disassembly to
verify whether this is the case or not, but I would leave that to
others).</p>
<p>My own preference would be for the unrolled version.</p>
<p>We could also make <i>C / Cphim</i>u a parameter to <i>pj_authset()</i>
(which could be named something else), since this could be used
for other conversions between auxiliary latitudes.<br>
Similarly, <i>pj_auth2geodlat() </i>could actually be used for
different conversions if passing it pre-computed coefficients for
other conversions, so perhaps it could have a more generic names.<br>
The rectifying latitude for <i>pj_enfn() </i>is a bit special
because it uses n^2 rather than n, which tripped me up for a
little while.</p>
<p>Thoughts / suggestions on how to move forward with this?<br>
<br>
As a next step I would prepare a Pull Request based on your
feedback, if you have a preference for the shared functions or the
unrolled loops approach.</p>
<p>Thank you very much for your help and guidance!</p>
<p>Kind regards,</p>
<p>-Jerome</p>
<div class="moz-cite-prefix">On 9/11/24 4:21 PM, Jérôme St-Louis
wrote:<br>
</div>
<blockquote type="cite"
cite="mid:859880f2-fba6-4e46-a615-ce8c2ac0cccf@ecere.com">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<p>So it seems like we already have an implementation of Horner
and Clenshaw in:</p>
<p> <a class="moz-txt-link-freetext"
href="https://github.com/OSGeo/PROJ/blob/master/src/mlfn.cpp"
moz-do-not-send="true">https://github.com/OSGeo/PROJ/blob/master/src/mlfn.cpp</a><br>
</p>
<div class="moz-signature">called <i>polyval()</i> and <i>clenshaw()</i>
just like in GeographicLib ( <a
href="https://github.com/geographiclib/geographiclib/blob/main/include/GeographicLib/Math.hpp#L280"
moz-do-not-send="true">polyval()</a> , <a
href="https://github.com/geographiclib/geographiclib/blob/main/src/AuxLatitude.cpp#L1319"
moz-do-not-send="true">Clenshaw()</a>).<br>
</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">It seems like Charles wrote or at least
updated that :)</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">That is using the Cµφ (C[mu phi]) (A5)
and Cφµ (C[phi mu]) (A6) from page 12 of the paper, where µ is
called the "rectifying latitude".<br>
I imagine that this is directly related to the "meridional
distance" ?</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">Perhaps we could re-organize this a bit
to share this <i>polyval()</i> and <i>clenshaw()</i> (they are
currently static functions local to this <i>mlfn.cpp</i>) for
use in <i>auth.cpp</i> ?</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">Thanks!</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">Kind regards,</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">-Jerome<br>
</div>
<div class="moz-signature"><br>
</div>
<div class="moz-cite-prefix">On 9/11/24 3:33 PM, Jérôme St-Louis
wrote:<br>
</div>
<blockquote type="cite"
cite="mid:9ab1b1a7-caae-43cc-97a5-ab5ddade3f26@ecere.com">
<meta http-equiv="Content-Type"
content="text/html; charset=UTF-8">
<p>Thanks a lot for the input Charles and Thomas,</p>
<p>I am not familiar with either <a
href="https://en.wikipedia.org/wiki/Horner%27s_method"
moz-do-not-send="true">Horner</a> or <a
href="https://en.wikipedia.org/wiki/Clenshaw_algorithm"
moz-do-not-send="true">Clenshaw</a>, but I do see the
mentions now on <i>Section 6 - Evaluating the series</i>
pages 6 and 7 of the papers.<br>
I implemented the simpler basic approach from section 3 / page
3, which also happened to more easily correspond to the
existing PROJ implementation.</p>
<p>I can definitely try to understand all this, with the help of
this Rust Geodesy code and the GeographicLib code, and have a
go at updating my proposed implementation for improved
accuracy and performance.</p>
<p>Kind regards,</p>
<p>-Jerome</p>
<div class="moz-cite-prefix">On 9/11/24 12:18 PM, Thomas Knudsen
wrote:<br>
</div>
<blockquote type="cite"
cite="mid:CAH0YoEMiYpNxR1LNYFxh9BqjM=NWWaiJKFZV_2yoe+kb__uy2g@mail.gmail.com">
<pre class="moz-quote-pre" wrap="">I totally agree with Charles regarding using Horner for polynomial
evaluation and Clenshaw for the trig series - for accuracy and speed.
I implemented all the material from Charles' preprint
<a class="moz-txt-link-freetext" href="https://arxiv.org/pdf/2212.05818"
moz-do-not-send="true">https://arxiv.org/pdf/2212.05818</a> for Rust Geodesy, when the preprint
appeared about 1½ years ago.
And although (being an experiment) my handling of the raw coefficients
is rather clumsy, at least it gave me a reason to revise my PROJ horner
and clenshaw implementations (which in turn were based on material from
Poder & Engsager: "Some Conformal Mappings...").
So Jérôme, perhaps take a look at the functions "taylor" and "fourier"
over at <a class="moz-txt-link-freetext"
href="https://github.com/busstoptaktik/geodesy/blob/main/src/math/series.rs"
moz-do-not-send="true">https://github.com/busstoptaktik/geodesy/blob/main/src/math/series.rs</a>
While written in Rust, translating to C++ should be rather trivial,
and they may be easier to follow than my decade-old versions already
in the PROJ code base.
</pre>
</blockquote>
</blockquote>
</blockquote>
</body>
</html>