boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-01-30 20:02:11 +00:00
Code Issues Projects Releases Wiki Activity

[formaulas][util] Define coefficient containers for computing series expansions

Browse Source 
The coefficient containers are defined as structs in
series_expansion.hpp file. They allow the caller to
compute expansions without specifying the size
for the output array.
This commit is contained in:
Adeel Ahmad
2018-06-25 17:19:44 +05:00
parent 1972bcda3e
commit dedccdbdae
2 changed files with 103 additions and 78 deletions
Download Patch File Download Diff File Expand all files Collapse all files

35
include/boost/geometry/formulas/karney_direct.hpp
View File

@@ -131,8 +131,7 @@ public:
= series_expansion::evaluate_series_A1<CT, SeriesOrder>(epsilon);
// Index zero element of coeffs_C1 is unused.
boost::array<CT, SeriesOrder + 1> coeffs_C1;
series_expansion::evaluate_coeffs_C1(epsilon, coeffs_C1);
series_expansion::coeffs_C1<CT, SeriesOrder> const coeffs_C1(epsilon);
// Tau is an integration variable.
CT const tau12 = distance / (b * (c1 + expansion_A1));
@@ -148,7 +147,7 @@ public:
math::normalize_values<CT>(sin_sigma1, cos_sigma1);
CT const B11 =
series_expansion::sin_cos_series<CT, SeriesOrder + 1>(sin_sigma1, cos_sigma1, coeffs_C1);
series_expansion::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C1);
CT const sin_B11 = sin(B11);
CT const cos_B11 = cos(B11);
@@ -158,11 +157,10 @@ public:
= cos_sigma1 * cos_B11 - sin_sigma1 * sin_B11;
// Index zero element of coeffs_C1p is unused.
boost::array<CT, SeriesOrder + 1> coeffs_C1p;
series_expansion::evaluate_coeffs_C1p(epsilon, coeffs_C1p);
series_expansion::coeffs_C1p<CT, SeriesOrder> const coeffs_C1p(epsilon);
CT const B12 =
- series_expansion::sin_cos_series<CT, SeriesOrder + 1>
- series_expansion::sin_cos_series
(sin_tau1 * cos_tau12 + cos_tau1 * sin_tau12,
cos_tau1 * cos_tau12 - sin_tau1 * sin_tau12,
coeffs_C1p);
@@ -207,26 +205,18 @@ public:
CT const omega12 = atan2(sin_omega2 * cos_omega1 - cos_omega2 * sin_omega1,
cos_omega2 * cos_omega1 + sin_omega2 * sin_omega1);
CT coeffs_A3[SeriesOrder];
series_expansion::evaluate_coeffs_A3<CT, SeriesOrder>(n, coeffs_A3);
series_expansion::coeffs_A3<CT, SeriesOrder> coeffs_A3(n);
CT const A3 = math::horner_evaluate(epsilon, coeffs_A3, coeffs_A3 + SeriesOrder);
CT const A3 = math::horner_evaluate(epsilon, coeffs_A3.begin(), coeffs_A3.end());
CT const A3c = -f * sin_alpha0 * A3;
// Compute the size of coefficient array.
size_t const coeffs_C3_size = (SeriesOrder * (SeriesOrder - 1)) / 2;
CT coeffs_C3x[coeffs_C3_size];
series_expansion::evaluate_coeffs_C3x<CT, SeriesOrder>(n, coeffs_C3x);
// Evaluate C3 coefficients.
CT coeffs_C3[SeriesOrder];
series_expansion::evaluate_coeffs_C3<CT, SeriesOrder>(epsilon, coeffs_C3, coeffs_C3x);
series_expansion::coeffs_C3<CT, SeriesOrder> coeffs_C3(n, epsilon);
CT const B31 =
series_expansion::sin_cos_series<CT, SeriesOrder>(sin_sigma1, cos_sigma1, coeffs_C3);
series_expansion::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C3);
CT const lam12 = omega12 + A3c *
(sigma12 + (series_expansion::sin_cos_series<CT, SeriesOrder>
(sigma12 + (series_expansion::sin_cos_series
(sin_sigma2,
cos_sigma2,
coeffs_C3) - B31));
@@ -248,13 +238,12 @@ public:
{
// Evaluate the coefficients for C2.
// Index zero element of coeffs_C2 is unused.
boost::array<CT, SeriesOrder + 1> coeffs_C2;
series_expansion::evaluate_coeffs_C2(epsilon, coeffs_C2);
series_expansion::coeffs_C2<CT, SeriesOrder> coeffs_C2(epsilon);
CT const B21 =
series_expansion::sin_cos_series<CT, SeriesOrder + 1>(sin_sigma1, cos_sigma1, coeffs_C2);
series_expansion::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C2);
CT const B22 =
series_expansion::sin_cos_series<CT, SeriesOrder + 1>(sin_sigma2, cos_sigma2, coeffs_C2);
series_expansion::sin_cos_series(sin_sigma2, cos_sigma2, coeffs_C2);
// Find the coefficients for A2 by computing the
// series expansion using Horner scehme.

146
include/boost/geometry/util/series_expansion.hpp
View File

@@ -28,6 +28,66 @@
namespace boost { namespace geometry { namespace series_expansion {
/*
The coefficient containers for the series expansions.
These structs allow the caller to only know the series order.
*/
template <typename CT, size_t SeriesOrder>
struct coeffs_C1 : boost::array<CT, SeriesOrder + 1>
{
coeffs_C1(CT const& epsilon)
{
evaluate_coeffs_C1(*this, epsilon);
}
};
template <typename CT, size_t SeriesOrder>
struct coeffs_C1p : boost::array<CT, SeriesOrder + 1>
{
coeffs_C1p(CT const& epsilon)
{
evaluate_coeffs_C1p(*this, epsilon);
}
};
template <typename CT, size_t SeriesOrder>
struct coeffs_C2 : boost::array<CT, SeriesOrder + 1>
{
coeffs_C2(CT const& epsilon)
{
evaluate_coeffs_C2(*this, epsilon);
}
};
template <typename CT, size_t SeriesOrder>
struct coeffs_C3x : boost::array<CT, (SeriesOrder * (SeriesOrder - 1)) / 2>
{
coeffs_C3x(CT const& n)
{
evaluate_coeffs_C3x(*this, SeriesOrder, n);
}
};
template <typename T, size_t SeriesOrder>
struct coeffs_C3 : boost::array<T, SeriesOrder>
{
coeffs_C3(T const& n, T const& epsilon)
{
coeffs_C3x<T, SeriesOrder> coeffs_C3x(n);
evaluate_coeffs_C3(*this, coeffs_C3x, epsilon);
}
};
template <typename CT, size_t SeriesOrder>
struct coeffs_A3 : boost::array<CT, SeriesOrder>
{
coeffs_A3(CT const& n)
{
evaluate_coeffs_A3(*this, n);
}
};
/*
Generate and evaluate the series expansion of the following integral
@@ -53,7 +113,7 @@ namespace boost { namespace geometry { namespace series_expansion {
sed -e 's/[0-9]\+/CT(&)/g; s/\[CT/\[/g; s/)\]/\]/g;
s/case\sCT(/case /g; s/):/:/g; s/epsCT(2)/eps2/g;'
*/
template <typename CT, std::size_t SeriesOrder>
template <typename CT, size_t SeriesOrder>
static inline CT evaluate_series_A1(CT eps)
{
CT eps2 = math::sqr(eps);
@@ -98,7 +158,7 @@ namespace boost { namespace geometry { namespace series_expansion {
generated by the following Maxima script:
geometry/doc/other/maxima/geod.mac
*/
template <typename CT, std::size_t SeriesOrder>
template <typename CT, size_t SeriesOrder>
CT evaluate_series_A2(CT const& eps)
{
CT const eps2 = math::sqr(eps);
@@ -116,7 +176,7 @@ namespace boost { namespace geometry { namespace series_expansion {
case 3:
t = eps2*((-CT(11)*eps2-CT(28))*eps2-CT(192))/CT(256);
break;
case 4:
default:
t = eps2*(eps2*((-CT(375)*eps2-CT(704))*eps2-CT(1792))-CT(12288))/CT(16384);
break;
}
@@ -143,11 +203,10 @@ namespace boost { namespace geometry { namespace series_expansion {
generated by the following Maxima script:
geometry/doc/other/maxima/geod.mac
*/
// TODO: this produces different results that geographiclib
template <typename CT, std::size_t SeriesOrder>
static inline void evaluate_coeffs_A3(CT const& n, CT c[])
template <typename Coeffs, typename CT>
static inline void evaluate_coeffs_A3(Coeffs &c, CT const& n)
{
switch (SeriesOrder) {
switch (int(Coeffs::static_size)) {
case 0:
break;
case 1:
@@ -192,7 +251,7 @@ namespace boost { namespace geometry { namespace series_expansion {
c[5] = (-CT(5)*n-CT(3))/CT(128);
c[6] = -CT(5)/CT(256);
break;
case 8:
default:
c[0] = CT(1);
c[1] = (n-CT(1))/CT(2);
c[2] = (n*(CT(3)*n-CT(1))-CT(2))/CT(8);
@@ -213,12 +272,12 @@ namespace boost { namespace geometry { namespace series_expansion {
generated by the following Maxima script:
geometry/doc/other/maxima/geod.mac
*/
template <typename CT, std::size_t SeriesOrder>
static inline void evaluate_coeffs_C1(CT eps, boost::array<CT, SeriesOrder>& c)
template <typename Coeffs, typename CT>
static inline void evaluate_coeffs_C1(Coeffs &c, CT const& eps)
{
CT eps2 = math::sqr(eps);
CT d = eps;
switch (c.size() - 1) {
switch (int(Coeffs::static_size) - 1) {
case 0:
break;
case 1:
@@ -284,7 +343,7 @@ namespace boost { namespace geometry { namespace series_expansion {
d *= eps;
c[7] = -CT(33)*d/CT(14336);
break;
case 8:
default:
c[1] = d*(eps2*(eps2*(CT(19)*eps2-CT(64))+CT(384))-CT(1024))/CT(2048);
d *= eps;
c[2] = d*(eps2*(eps2*(CT(7)*eps2-CT(18))+CT(128))-CT(256))/CT(4096);
@@ -304,6 +363,7 @@ namespace boost { namespace geometry { namespace series_expansion {
}
}
/*
The coefficients C1p[l] in the Fourier expansion of B1p.
@@ -312,12 +372,12 @@ namespace boost { namespace geometry { namespace series_expansion {
generated by the following Maxima script:
geometry/doc/other/maxima/geod.mac
*/
template <typename CT, std::size_t SeriesOrder>
static inline void evaluate_coeffs_C1p(CT eps, boost::array<CT, SeriesOrder>& c)
template <typename Coeffs, typename CT>
static inline void evaluate_coeffs_C1p(Coeffs& c, CT const& eps)
{
CT const eps2 = math::sqr(eps);
CT d = eps;
switch (c.size() - 1) {
switch (int(Coeffs::static_size) - 1) {
case 0:
break;
case 1:
@@ -383,7 +443,7 @@ namespace boost { namespace geometry { namespace series_expansion {
d *= eps;
c[7] = CT(459485)*d/CT(516096);
break;
case 8:
default:
c[1] = d*(eps2*((CT(9840)-CT(4879)*eps2)*eps2-CT(20736))+CT(36864))/CT(73728);
d *= eps;
c[2] = d*(eps2*((CT(120150)-CT(86171)*eps2)*eps2-CT(142080))+CT(115200))/CT(368640);
@@ -411,12 +471,12 @@ namespace boost { namespace geometry { namespace series_expansion {
generated by the following Maxima script:
geometry/doc/other/maxima/geod.mac
*/
template <typename CT, std::size_t SeriesOrder>
static inline void evaluate_coeffs_C2(CT const& eps, boost::array<CT, SeriesOrder>& c)
template <typename Coeffs, typename CT>
static inline void evaluate_coeffs_C2(Coeffs& c, CT const& eps)
{
CT const eps2 = math::sqr(eps);
CT d = eps;
switch (c.size() - 1) {
switch (int(Coeffs::static_size) - 1) {
case 0:
break;
case 1:
@@ -482,7 +542,7 @@ namespace boost { namespace geometry { namespace series_expansion {
d *= eps;
c[7] = CT(429)*d/CT(14336);
break;
case 8:
default:
c[1] = d*(eps2*(eps2*(CT(41)*eps2+CT(64))+CT(128))+CT(1024))/CT(2048);
d *= eps;
c[2] = d*(eps2*(eps2*(CT(47)*eps2+CT(70))+CT(128))+CT(768))/CT(4096);
@@ -510,8 +570,8 @@ namespace boost { namespace geometry { namespace series_expansion {
generated by the following Maxima script:
geometry/doc/other/maxima/geod.mac
*/
template <typename CT, size_t SeriesOrder>
static inline void evaluate_coeffs_C3x(CT const& n, CT c[]) {
template <typename Coeffs, typename CT>
static inline void evaluate_coeffs_C3x(Coeffs &c, size_t const& SeriesOrder, CT const& n) {
const CT n2 = math::sqr(n);
switch (SeriesOrder) {
case 0:
@@ -586,7 +646,7 @@ namespace boost { namespace geometry { namespace series_expansion {
c[19] = (CT(9)-CT(15)*n)/CT(1024);
c[20] = (CT(44)-CT(99)*n)/CT(8192);
break;
case 8:
default:
c[0] = (CT(1)-n)/CT(4);
c[1] = (CT(1)-n2)/CT(8);
c[2] = (n*((-CT(5)*n-CT(1))*n+CT(3))+CT(3))/CT(64);
@@ -625,19 +685,19 @@ namespace boost { namespace geometry { namespace series_expansion {
Elements coeffs1[1] through coeffs1[SeriesOrder - 1] are set.
*/
template <typename CT, size_t SeriesOrder>
static inline void evaluate_coeffs_C3(CT eps, CT coeffs1[], CT coeffs2[])
template <typename Coeffs1, typename Coeffs2, typename CT>
static inline void evaluate_coeffs_C3(Coeffs1 &coeffs1, Coeffs2 &coeffs2, CT const& eps)
{
CT mult = 1;
int offset = 0;
// l is the index of C3[l].
for (size_t l = 1; l < SeriesOrder; ++l)
for (size_t l = 1; l < Coeffs1::static_size; ++l)
{
// Order of polynomial in eps.
int m = SeriesOrder - l - 1;
int m = Coeffs1::static_size - l - 1;
mult *= eps;
coeffs1[l] = mult * math::polyval(m, coeffs2 + offset, eps);
coeffs1[l] = mult * math::polyval(m, coeffs2.begin() + offset, eps);
offset += m + 1;
}
// Post condition: offset == coeffs_C3_size
@@ -650,33 +710,10 @@ namespace boost { namespace geometry { namespace series_expansion {
using Clenshaw summation.
*/
template <typename CT, size_t SeriesOrder>
static inline CT sin_cos_series(CT sinx, CT cosx, const CT coeffs[])
template <typename CT, typename Coeffs>
static inline CT sin_cos_series(CT const& sinx, CT const& cosx, Coeffs const& coeffs)
{
size_t n = SeriesOrder;
// Point to one beyond last element.
coeffs += (n + 1);
CT ar = 2 * (cosx - sinx) * (cosx + sinx);
CT k0 = n & 1 ? *--coeffs : 0;
CT k1 = 0;
// Make n even.
n /= 2;
while (n--) {
// Unroll loop x 2, so accumulators return to their original role.
k1 = ar * k0 - k1 + *--coeffs;
k0 = ar * k1 - k0 + *--coeffs;
}
return 2 * sinx * cosx * k0;
}
template <typename CT, size_t SeriesOrder>
static inline CT sin_cos_series(CT sinx, CT cosx, const boost::array<CT, SeriesOrder> coeffs)
{
size_t n = SeriesOrder - 1;
size_t n = Coeffs::static_size - 1;
size_t index = 0;
// Point to one beyond last element.
@@ -697,7 +734,6 @@ namespace boost { namespace geometry { namespace series_expansion {
return 2 * sinx * cosx * k0;
}
}}} // namespace boost::geometry::series_expansion
#endif // BOOST_GEOMETRY_UTIL_SERIES_EXPANSION_HPP