diff --git a/include/boost/geometry/formulas/karney_direct.hpp b/include/boost/geometry/formulas/karney_direct.hpp index c844f977c..df9aaabe4 100644 --- a/include/boost/geometry/formulas/karney_direct.hpp +++ b/include/boost/geometry/formulas/karney_direct.hpp @@ -131,8 +131,7 @@ public: = series_expansion::evaluate_series_A1(epsilon); // Index zero element of coeffs_C1 is unused. - boost::array coeffs_C1; - series_expansion::evaluate_coeffs_C1(epsilon, coeffs_C1); + series_expansion::coeffs_C1 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(sin_sigma1, cos_sigma1); CT const B11 = - series_expansion::sin_cos_series(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 coeffs_C1p; - series_expansion::evaluate_coeffs_C1p(epsilon, coeffs_C1p); + series_expansion::coeffs_C1p const coeffs_C1p(epsilon); CT const B12 = - - series_expansion::sin_cos_series + - 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(n, coeffs_A3); + series_expansion::coeffs_A3 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(n, coeffs_C3x); - - // Evaluate C3 coefficients. - CT coeffs_C3[SeriesOrder]; - series_expansion::evaluate_coeffs_C3(epsilon, coeffs_C3, coeffs_C3x); + series_expansion::coeffs_C3 coeffs_C3(n, epsilon); CT const B31 = - series_expansion::sin_cos_series(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 + (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 coeffs_C2; - series_expansion::evaluate_coeffs_C2(epsilon, coeffs_C2); + series_expansion::coeffs_C2 coeffs_C2(epsilon); CT const B21 = - series_expansion::sin_cos_series(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(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. diff --git a/include/boost/geometry/util/series_expansion.hpp b/include/boost/geometry/util/series_expansion.hpp index 492393e19..f9a063c5c 100644 --- a/include/boost/geometry/util/series_expansion.hpp +++ b/include/boost/geometry/util/series_expansion.hpp @@ -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 + struct coeffs_C1 : boost::array + { + coeffs_C1(CT const& epsilon) + { + evaluate_coeffs_C1(*this, epsilon); + } + }; + + template + struct coeffs_C1p : boost::array + { + coeffs_C1p(CT const& epsilon) + { + evaluate_coeffs_C1p(*this, epsilon); + } + }; + + template + struct coeffs_C2 : boost::array + { + coeffs_C2(CT const& epsilon) + { + evaluate_coeffs_C2(*this, epsilon); + } + }; + + template + struct coeffs_C3x : boost::array + { + coeffs_C3x(CT const& n) + { + evaluate_coeffs_C3x(*this, SeriesOrder, n); + } + }; + + template + struct coeffs_C3 : boost::array + { + coeffs_C3(T const& n, T const& epsilon) + { + coeffs_C3x coeffs_C3x(n); + + evaluate_coeffs_C3(*this, coeffs_C3x, epsilon); + } + }; + + template + struct coeffs_A3 : boost::array + { + 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 + template 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 + template 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 - static inline void evaluate_coeffs_A3(CT const& n, CT c[]) + template + 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 - static inline void evaluate_coeffs_C1(CT eps, boost::array& c) + template + 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 - static inline void evaluate_coeffs_C1p(CT eps, boost::array& c) + template + 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 - static inline void evaluate_coeffs_C2(CT const& eps, boost::array& c) + template + 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 - static inline void evaluate_coeffs_C3x(CT const& n, CT c[]) { + template + 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 - static inline void evaluate_coeffs_C3(CT eps, CT coeffs1[], CT coeffs2[]) + template + 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 - static inline CT sin_cos_series(CT sinx, CT cosx, const CT coeffs[]) + template + 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 - static inline CT sin_cos_series(CT sinx, CT cosx, const boost::array 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