Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:36:53 

0001 // Boost.Geometry
0002 
0003 // Copyright (c) 2018 Adeel Ahmad, Islamabad, Pakistan.
0004 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0005 
0006 // Contributed and/or modified by Adeel Ahmad, as part of Google Summer of Code 2018 program.
0007 
0008 // This file was modified by Oracle on 2019.
0009 // Modifications copyright (c) 2019 Oracle and/or its affiliates.
0010 
0011 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0012 
0013 // Use, modification and distribution is subject to the Boost Software License,
0014 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0015 // http://www.boost.org/LICENSE_1_0.txt)
0016 
0017 // This file is converted from GeographicLib, https://geographiclib.sourceforge.io
0018 // GeographicLib is originally written by Charles Karney.
0019 
0020 // Author: Charles Karney (2008-2017)
0021 
0022 // Last updated version of GeographicLib: 1.49
0023 
0024 // Original copyright notice:
0025 
0026 // Copyright (c) Charles Karney (2008-2017) <charles@karney.com> and licensed
0027 // under the MIT/X11 License. For more information, see
0028 // https://geographiclib.sourceforge.io
0029 
0030 #ifndef BOOST_GEOMETRY_UTIL_SERIES_EXPANSION_HPP
0031 #define BOOST_GEOMETRY_UTIL_SERIES_EXPANSION_HPP
0032 
0033 #include <boost/array.hpp>
0034 #include <boost/geometry/core/assert.hpp>
0035 #include <boost/geometry/util/math.hpp>
0036 
0037 namespace boost { namespace geometry { namespace series_expansion {
0038 
0039     /*
0040      Generate and evaluate the series expansion of the following integral
0041 
0042      I1 = integrate( sqrt(1+k2*sin(sigma1)^2), sigma1, 0, sigma )
0043 
0044      which is valid for k2 small. We substitute k2 = 4 * eps / (1 - eps)^2
0045      and expand (1 - eps) * I1 retaining terms up to order eps^maxpow
0046      in A1 and C1[l].
0047 
0048      The resulting series is of the form
0049 
0050      A1 * ( sigma + sum(C1[l] * sin(2*l*sigma), l, 1, maxpow) ).
0051 
0052      The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
0053 
0054      The expansion above is performed in Maxima, a Computer Algebra System.
0055      The C++ code (that yields the function evaluate_A1 below) is
0056      generated by the following Maxima script:
0057      geometry/doc/other/maxima/geod.mac
0058 
0059      To replace each number x by CT(x) the following
0060      script can be used:
0061        sed -e 's/[0-9]\+/CT(&)/g; s/\[CT/\[/g; s/)\]/\]/g;
0062                s/case\sCT(/case /g; s/):/:/g; s/epsCT(2)/eps2/g;'
0063     */
0064     template <size_t SeriesOrder, typename CT>
0065     inline CT evaluate_A1(CT const& eps)
0066     {
0067         CT const eps2 = math::sqr(eps);
0068         CT t;
0069         switch (SeriesOrder/2)
0070         {
0071         case 0:
0072             t = CT(0);
0073             break;
0074         case 1:
0075             t = eps2/CT(4);
0076             break;
0077         case 2:
0078             t = eps2*(eps2+CT(16))/CT(64);
0079             break;
0080         case 3:
0081             t = eps2*(eps2*(eps2+CT(4))+CT(64))/CT(256);
0082             break;
0083         default:
0084             t = eps2*(eps2*(eps2*(CT(25)*eps2+CT(64))+CT(256))+CT(4096))/CT(16384);
0085             break;
0086         }
0087         return (t + eps) / (CT(1) - eps);
0088     }
0089 
0090     /*
0091      Generate and evaluate the series expansion of the following integral
0092 
0093      I2 = integrate( 1/sqrt(1+k2*sin(sigma1)^2), sigma1, 0, sigma )
0094 
0095      which is valid for k2 small. We substitute k2 = 4 * eps / (1 - eps)^2
0096      and expand (1 - eps) * I2 retaining terms up to order eps^maxpow
0097      in A2 and C2[l].
0098 
0099      The resulting series is of the form
0100 
0101      A2 * ( sigma + sum(C2[l] * sin(2*l*sigma), l, 1, maxpow) )
0102 
0103      The scale factor A2-1 = mean value of (d/dsigma)2 - 1
0104 
0105      The expansion above is performed in Maxima, a Computer Algebra System.
0106      The C++ code (that yields the function evaluate_A2 below) is
0107      generated by the following Maxima script:
0108      geometry/doc/other/maxima/geod.mac
0109     */
0110     template <size_t SeriesOrder, typename CT>
0111     inline CT evaluate_A2(CT const& eps)
0112     {
0113         CT const eps2 = math::sqr(eps);
0114         CT t;
0115         switch (SeriesOrder/2)
0116         {
0117         case 0:
0118             t = CT(0);
0119             break;
0120         case 1:
0121             t = -CT(3)*eps2/CT(4);
0122             break;
0123         case 2:
0124             t = (-CT(7)*eps2-CT(48))*eps2/CT(64);
0125             break;
0126         case 3:
0127             t = eps2*((-CT(11)*eps2-CT(28))*eps2-CT(192))/CT(256);
0128             break;
0129         default:
0130             t = eps2*(eps2*((-CT(375)*eps2-CT(704))*eps2-CT(1792))-CT(12288))/CT(16384);
0131             break;
0132         }
0133         return (t - eps) / (CT(1) + eps);
0134     }
0135 
0136     /*
0137      Express
0138 
0139         I3 = integrate( (2-f)/(1+(1-f)*sqrt(1+k2*sin(sigma1)^2)), sigma1, 0, sigma )
0140 
0141      as a series
0142 
0143         A3 * ( sigma + sum(C3[l] * sin(2*l*sigma), l, 1, maxpow-1) )
0144 
0145      valid for f and k2 small.  It is convenient to write k2 = 4 * eps / (1 -
0146      eps)^2 and f = 2*n/(1+n) and expand in eps and n.  This procedure leads
0147      to a series where the coefficients of eps^j are terminating series in n.
0148 
0149      The scale factor A3 = mean value of (d/dsigma)I3
0150 
0151      The expansion above is performed in Maxima, a Computer Algebra System.
0152      The C++ code (that yields the function evaluate_coeffs_A3 below) is
0153      generated by the following Maxima script:
0154      geometry/doc/other/maxima/geod.mac
0155     */
0156     template <typename Coeffs, typename CT>
0157     inline void evaluate_coeffs_A3(Coeffs &c, CT const& n)
0158     {
0159         switch (int(Coeffs::static_size))
0160         {
0161         case 0:
0162             break;
0163         case 1:
0164             c[0] = CT(1);
0165             break;
0166         case 2:
0167             c[0] = CT(1);
0168             c[1] = -CT(1)/CT(2);
0169             break;
0170         case 3:
0171             c[0] = CT(1);
0172             c[1] = (n-CT(1))/CT(2);
0173             c[2] = -CT(1)/CT(4);
0174             break;
0175         case 4:
0176             c[0] = CT(1);
0177             c[1] = (n-CT(1))/CT(2);
0178             c[2] = (-n-CT(2))/CT(8);
0179             c[3] = -CT(1)/CT(16);
0180             break;
0181         case 5:
0182             c[0] = CT(1);
0183             c[1] = (n-CT(1))/CT(2);
0184             c[2] = (n*(CT(3)*n-CT(1))-CT(2))/CT(8);
0185             c[3] = (-CT(3)*n-CT(1))/CT(16);
0186             c[4] = -CT(3)/CT(64);
0187             break;
0188         case 6:
0189             c[0] = CT(1);
0190             c[1] = (n-CT(1))/CT(2);
0191             c[2] = (n*(CT(3)*n-CT(1))-CT(2))/CT(8);
0192             c[3] = ((-n-CT(3))*n-CT(1))/CT(16);
0193             c[4] = (-CT(2)*n-CT(3))/CT(64);
0194             c[5] = -CT(3)/CT(128);
0195             break;
0196         case 7:
0197             c[0] = CT(1);
0198             c[1] = (n-CT(1))/CT(2);
0199             c[2] = (n*(CT(3)*n-CT(1))-CT(2))/CT(8);
0200             c[3] = (n*(n*(CT(5)*n-CT(1))-CT(3))-CT(1))/CT(16);
0201             c[4] = ((-CT(10)*n-CT(2))*n-CT(3))/CT(64);
0202             c[5] = (-CT(5)*n-CT(3))/CT(128);
0203             c[6] = -CT(5)/CT(256);
0204             break;
0205         default:
0206             c[0] = CT(1);
0207             c[1] = (n-CT(1))/CT(2);
0208             c[2] = (n*(CT(3)*n-CT(1))-CT(2))/CT(8);
0209             c[3] = (n*(n*(CT(5)*n-CT(1))-CT(3))-CT(1))/CT(16);
0210             c[4] = (n*((-CT(5)*n-CT(20))*n-CT(4))-CT(6))/CT(128);
0211             c[5] = ((-CT(5)*n-CT(10))*n-CT(6))/CT(256);
0212             c[6] = (-CT(15)*n-CT(20))/CT(1024);
0213             c[7] = -CT(25)/CT(2048);
0214             break;
0215         }
0216     }
0217 
0218     /*
0219      The coefficients C1[l] in the Fourier expansion of B1.
0220 
0221      The expansion below is performed in Maxima, a Computer Algebra System.
0222      The C++ code (that yields the function evaluate_coeffs_C1 below) is
0223      generated by the following Maxima script:
0224      geometry/doc/other/maxima/geod.mac
0225     */
0226     template <typename Coeffs, typename CT>
0227     inline void evaluate_coeffs_C1(Coeffs &c, CT const& eps)
0228     {
0229         CT const eps2 = math::sqr(eps);
0230         CT d = eps;
0231         switch (int(Coeffs::static_size) - 1)
0232         {
0233         case 0:
0234             break;
0235         case 1:
0236             c[1] = -d/CT(2);
0237             break;
0238         case 2:
0239             c[1] = -d/CT(2);
0240             d *= eps;
0241             c[2] = -d/CT(16);
0242             break;
0243         case 3:
0244             c[1] = d*(CT(3)*eps2-CT(8))/CT(16);
0245             d *= eps;
0246             c[2] = -d/CT(16);
0247             d *= eps;
0248             c[3] = -d/CT(48);
0249             break;
0250         case 4:
0251             c[1] = d*(CT(3)*eps2-CT(8))/CT(16);
0252             d *= eps;
0253             c[2] = d*(eps2-CT(2))/CT(32);
0254             d *= eps;
0255             c[3] = -d/CT(48);
0256             d *= eps;
0257             c[4] = -CT(5)*d/CT(512);
0258             break;
0259         case 5:
0260             c[1] = d*((CT(6)-eps2)*eps2-CT(16))/CT(32);
0261             d *= eps;
0262             c[2] = d*(eps2-CT(2))/CT(32);
0263             d *= eps;
0264             c[3] = d*(CT(9)*eps2-CT(16))/CT(768);
0265             d *= eps;
0266             c[4] = -CT(5)*d/CT(512);
0267             d *= eps;
0268             c[5] = -CT(7)*d/CT(1280);
0269             break;
0270         case 6:
0271             c[1] = d*((CT(6)-eps2)*eps2-CT(16))/CT(32);
0272             d *= eps;
0273             c[2] = d*((CT(64)-CT(9)*eps2)*eps2-CT(128))/CT(2048);
0274             d *= eps;
0275             c[3] = d*(CT(9)*eps2-CT(16))/CT(768);
0276             d *= eps;
0277             c[4] = d*(CT(3)*eps2-CT(5))/CT(512);
0278             d *= eps;
0279             c[5] = -CT(7)*d/CT(1280);
0280             d *= eps;
0281             c[6] = -CT(7)*d/CT(2048);
0282             break;
0283         case 7:
0284             c[1] = d*(eps2*(eps2*(CT(19)*eps2-CT(64))+CT(384))-CT(1024))/CT(2048);
0285             d *= eps;
0286             c[2] = d*((CT(64)-CT(9)*eps2)*eps2-CT(128))/CT(2048);
0287             d *= eps;
0288             c[3] = d*((CT(72)-CT(9)*eps2)*eps2-CT(128))/CT(6144);
0289             d *= eps;
0290             c[4] = d*(CT(3)*eps2-CT(5))/CT(512);
0291             d *= eps;
0292             c[5] = d*(CT(35)*eps2-CT(56))/CT(10240);
0293             d *= eps;
0294             c[6] = -CT(7)*d/CT(2048);
0295             d *= eps;
0296             c[7] = -CT(33)*d/CT(14336);
0297             break;
0298         default:
0299             c[1] = d*(eps2*(eps2*(CT(19)*eps2-CT(64))+CT(384))-CT(1024))/CT(2048);
0300             d *= eps;
0301             c[2] = d*(eps2*(eps2*(CT(7)*eps2-CT(18))+CT(128))-CT(256))/CT(4096);
0302             d *= eps;
0303             c[3] = d*((CT(72)-CT(9)*eps2)*eps2-CT(128))/CT(6144);
0304             d *= eps;
0305             c[4] = d*((CT(96)-CT(11)*eps2)*eps2-CT(160))/CT(16384);
0306             d *= eps;
0307             c[5] = d*(CT(35)*eps2-CT(56))/CT(10240);
0308             d *= eps;
0309             c[6] = d*(CT(9)*eps2-CT(14))/CT(4096);
0310             d *= eps;
0311             c[7] = -CT(33)*d/CT(14336);
0312             d *= eps;
0313             c[8] = -CT(429)*d/CT(262144);
0314             break;
0315         }
0316     }
0317 
0318     /*
0319      The coefficients C1p[l] in the Fourier expansion of B1p.
0320 
0321      The expansion below is performed in Maxima, a Computer Algebra System.
0322      The C++ code (that yields the function evaluate_coeffs_C1p below) is
0323      generated by the following Maxima script:
0324      geometry/doc/other/maxima/geod.mac
0325     */
0326     template <typename Coeffs, typename CT>
0327     inline void evaluate_coeffs_C1p(Coeffs& c, CT const& eps)
0328     {
0329         CT const eps2 = math::sqr(eps);
0330         CT d = eps;
0331         switch (int(Coeffs::static_size) - 1)
0332         {
0333         case 0:
0334             break;
0335         case 1:
0336             c[1] = d/CT(2);
0337             break;
0338         case 2:
0339             c[1] = d/CT(2);
0340             d *= eps;
0341             c[2] = CT(5)*d/CT(16);
0342             break;
0343         case 3:
0344             c[1] = d*(CT(16)-CT(9)*eps2)/CT(32);
0345             d *= eps;
0346             c[2] = CT(5)*d/CT(16);
0347             d *= eps;
0348             c[3] = CT(29)*d/CT(96);
0349             break;
0350         case 4:
0351             c[1] = d*(CT(16)-CT(9)*eps2)/CT(32);
0352             d *= eps;
0353             c[2] = d*(CT(30)-CT(37)*eps2)/CT(96);
0354             d *= eps;
0355             c[3] = CT(29)*d/CT(96);
0356             d *= eps;
0357             c[4] = CT(539)*d/CT(1536);
0358             break;
0359         case 5:
0360             c[1] = d*(eps2*(CT(205)*eps2-CT(432))+CT(768))/CT(1536);
0361             d *= eps;
0362             c[2] = d*(CT(30)-CT(37)*eps2)/CT(96);
0363             d *= eps;
0364             c[3] = d*(CT(116)-CT(225)*eps2)/CT(384);
0365             d *= eps;
0366             c[4] = CT(539)*d/CT(1536);
0367             d *= eps;
0368             c[5] = CT(3467)*d/CT(7680);
0369             break;
0370         case 6:
0371             c[1] = d*(eps2*(CT(205)*eps2-CT(432))+CT(768))/CT(1536);
0372             d *= eps;
0373             c[2] = d*(eps2*(CT(4005)*eps2-CT(4736))+CT(3840))/CT(12288);
0374             d *= eps;
0375             c[3] = d*(CT(116)-CT(225)*eps2)/CT(384);
0376             d *= eps;
0377             c[4] = d*(CT(2695)-CT(7173)*eps2)/CT(7680);
0378             d *= eps;
0379             c[5] = CT(3467)*d/CT(7680);
0380             d *= eps;
0381             c[6] = CT(38081)*d/CT(61440);
0382             break;
0383         case 7:
0384             c[1] = d*(eps2*((CT(9840)-CT(4879)*eps2)*eps2-CT(20736))+CT(36864))/CT(73728);
0385             d *= eps;
0386             c[2] = d*(eps2*(CT(4005)*eps2-CT(4736))+CT(3840))/CT(12288);
0387             d *= eps;
0388             c[3] = d*(eps2*(CT(8703)*eps2-CT(7200))+CT(3712))/CT(12288);
0389             d *= eps;
0390             c[4] = d*(CT(2695)-CT(7173)*eps2)/CT(7680);
0391             d *= eps;
0392             c[5] = d*(CT(41604)-CT(141115)*eps2)/CT(92160);
0393             d *= eps;
0394             c[6] = CT(38081)*d/CT(61440);
0395             d *= eps;
0396             c[7] = CT(459485)*d/CT(516096);
0397             break;
0398         default:
0399             c[1] = d*(eps2*((CT(9840)-CT(4879)*eps2)*eps2-CT(20736))+CT(36864))/CT(73728);
0400             d *= eps;
0401             c[2] = d*(eps2*((CT(120150)-CT(86171)*eps2)*eps2-CT(142080))+CT(115200))/CT(368640);
0402             d *= eps;
0403             c[3] = d*(eps2*(CT(8703)*eps2-CT(7200))+CT(3712))/CT(12288);
0404             d *= eps;
0405             c[4] = d*(eps2*(CT(1082857)*eps2-CT(688608))+CT(258720))/CT(737280);
0406             d *= eps;
0407             c[5] = d*(CT(41604)-CT(141115)*eps2)/CT(92160);
0408             d *= eps;
0409             c[6] = d*(CT(533134)-CT(2200311)*eps2)/CT(860160);
0410             d *= eps;
0411             c[7] = CT(459485)*d/CT(516096);
0412             d *= eps;
0413             c[8] = CT(109167851)*d/CT(82575360);
0414             break;
0415         }
0416     }
0417 
0418     /*
0419      The coefficients C2[l] in the Fourier expansion of B2.
0420 
0421      The expansion below is performed in Maxima, a Computer Algebra System.
0422      The C++ code (that yields the function evaluate_coeffs_C2 below) is
0423      generated by the following Maxima script:
0424      geometry/doc/other/maxima/geod.mac
0425     */
0426     template <typename Coeffs, typename CT>
0427     inline void evaluate_coeffs_C2(Coeffs& c, CT const& eps)
0428     {
0429         CT const eps2 = math::sqr(eps);
0430         CT d = eps;
0431         switch (int(Coeffs::static_size) - 1)
0432         {
0433         case 0:
0434             break;
0435         case 1:
0436             c[1] = d/CT(2);
0437             break;
0438         case 2:
0439             c[1] = d/CT(2);
0440             d *= eps;
0441             c[2] = CT(3)*d/CT(16);
0442             break;
0443         case 3:
0444             c[1] = d*(eps2+CT(8))/CT(16);
0445             d *= eps;
0446             c[2] = CT(3)*d/CT(16);
0447             d *= eps;
0448             c[3] = CT(5)*d/CT(48);
0449             break;
0450         case 4:
0451             c[1] = d*(eps2+CT(8))/CT(16);
0452             d *= eps;
0453             c[2] = d*(eps2+CT(6))/CT(32);
0454             d *= eps;
0455             c[3] = CT(5)*d/CT(48);
0456             d *= eps;
0457             c[4] = CT(35)*d/CT(512);
0458             break;
0459         case 5:
0460             c[1] = d*(eps2*(eps2+CT(2))+CT(16))/CT(32);
0461             d *= eps;
0462             c[2] = d*(eps2+CT(6))/CT(32);
0463             d *= eps;
0464             c[3] = d*(CT(15)*eps2+CT(80))/CT(768);
0465             d *= eps;
0466             c[4] = CT(35)*d/CT(512);
0467             d *= eps;
0468             c[5] = CT(63)*d/CT(1280);
0469             break;
0470         case 6:
0471             c[1] = d*(eps2*(eps2+CT(2))+CT(16))/CT(32);
0472             d *= eps;
0473             c[2] = d*(eps2*(CT(35)*eps2+CT(64))+CT(384))/CT(2048);
0474             d *= eps;
0475             c[3] = d*(CT(15)*eps2+CT(80))/CT(768);
0476             d *= eps;
0477             c[4] = d*(CT(7)*eps2+CT(35))/CT(512);
0478             d *= eps;
0479             c[5] = CT(63)*d/CT(1280);
0480             d *= eps;
0481             c[6] = CT(77)*d/CT(2048);
0482             break;
0483         case 7:
0484             c[1] = d*(eps2*(eps2*(CT(41)*eps2+CT(64))+CT(128))+CT(1024))/CT(2048);
0485             d *= eps;
0486             c[2] = d*(eps2*(CT(35)*eps2+CT(64))+CT(384))/CT(2048);
0487             d *= eps;
0488             c[3] = d*(eps2*(CT(69)*eps2+CT(120))+CT(640))/CT(6144);
0489             d *= eps;
0490             c[4] = d*(CT(7)*eps2+CT(35))/CT(512);
0491             d *= eps;
0492             c[5] = d*(CT(105)*eps2+CT(504))/CT(10240);
0493             d *= eps;
0494             c[6] = CT(77)*d/CT(2048);
0495             d *= eps;
0496             c[7] = CT(429)*d/CT(14336);
0497             break;
0498         default:
0499             c[1] = d*(eps2*(eps2*(CT(41)*eps2+CT(64))+CT(128))+CT(1024))/CT(2048);
0500             d *= eps;
0501             c[2] = d*(eps2*(eps2*(CT(47)*eps2+CT(70))+CT(128))+CT(768))/CT(4096);
0502             d *= eps;
0503             c[3] = d*(eps2*(CT(69)*eps2+CT(120))+CT(640))/CT(6144);
0504             d *= eps;
0505             c[4] = d*(eps2*(CT(133)*eps2+CT(224))+CT(1120))/CT(16384);
0506             d *= eps;
0507             c[5] = d*(CT(105)*eps2+CT(504))/CT(10240);
0508             d *= eps;
0509             c[6] = d*(CT(33)*eps2+CT(154))/CT(4096);
0510             d *= eps;
0511             c[7] = CT(429)*d/CT(14336);
0512             d *= eps;
0513             c[8] = CT(6435)*d/CT(262144);
0514             break;
0515         }
0516     }
0517 
0518     /*
0519      The coefficients C3[l] in the Fourier expansion of B3.
0520 
0521      The expansion below is performed in Maxima, a Computer Algebra System.
0522      The C++ code (that yields the function evaluate_coeffs_C3 below) is
0523      generated by the following Maxima script:
0524      geometry/doc/other/maxima/geod.mac
0525     */
0526     template <size_t SeriesOrder, typename Coeffs, typename CT>
0527     inline void evaluate_coeffs_C3x(Coeffs &c, CT const& n) {
0528         BOOST_GEOMETRY_ASSERT((Coeffs::static_size == (SeriesOrder * (SeriesOrder - 1)) / 2));
0529 
0530         CT const n2 = math::sqr(n);
0531         switch (SeriesOrder)
0532         {
0533         case 0:
0534             break;
0535         case 1:
0536             break;
0537         case 2:
0538             c[0] = (CT(1)-n)/CT(4);
0539             break;
0540         case 3:
0541             c[0] = (CT(1)-n)/CT(4);
0542             c[1] = (CT(1)-n2)/CT(8);
0543             c[2] = ((n-CT(3))*n+CT(2))/CT(32);
0544             break;
0545         case 4:
0546             c[0] = (CT(1)-n)/CT(4);
0547             c[1] = (CT(1)-n2)/CT(8);
0548             c[2] = (n*((-CT(5)*n-CT(1))*n+CT(3))+CT(3))/CT(64);
0549             c[3] = ((n-CT(3))*n+CT(2))/CT(32);
0550             c[4] = (n*(n*(CT(2)*n-CT(3))-CT(2))+CT(3))/CT(64);
0551             c[5] = (n*((CT(5)-n)*n-CT(9))+CT(5))/CT(192);
0552             break;
0553         case 5:
0554             c[0] = (CT(1)-n)/CT(4);
0555             c[1] = (CT(1)-n2)/CT(8);
0556             c[2] = (n*((-CT(5)*n-CT(1))*n+CT(3))+CT(3))/CT(64);
0557             c[3] = (n*((CT(2)-CT(2)*n)*n+CT(2))+CT(5))/CT(128);
0558             c[4] = ((n-CT(3))*n+CT(2))/CT(32);
0559             c[5] = (n*(n*(CT(2)*n-CT(3))-CT(2))+CT(3))/CT(64);
0560             c[6] = (n*((-CT(6)*n-CT(9))*n+CT(2))+CT(6))/CT(256);
0561             c[7] = (n*((CT(5)-n)*n-CT(9))+CT(5))/CT(192);
0562             c[8] = (n*(n*(CT(10)*n-CT(6))-CT(10))+CT(9))/CT(384);
0563             c[9] = (n*((CT(20)-CT(7)*n)*n-CT(28))+CT(14))/CT(1024);
0564             break;
0565         case 6:
0566             c[0] = (CT(1)-n)/CT(4);
0567             c[1] = (CT(1)-n2)/CT(8);
0568             c[2] = (n*((-CT(5)*n-CT(1))*n+CT(3))+CT(3))/CT(64);
0569             c[3] = (n*((CT(2)-CT(2)*n)*n+CT(2))+CT(5))/CT(128);
0570             c[4] = (n*(CT(3)*n+CT(11))+CT(12))/CT(512);
0571             c[5] = ((n-CT(3))*n+CT(2))/CT(32);
0572             c[6] = (n*(n*(CT(2)*n-CT(3))-CT(2))+CT(3))/CT(64);
0573             c[7] = (n*((-CT(6)*n-CT(9))*n+CT(2))+CT(6))/CT(256);
0574             c[8] = ((CT(1)-CT(2)*n)*n+CT(5))/CT(256);
0575             c[9] = (n*((CT(5)-n)*n-CT(9))+CT(5))/CT(192);
0576             c[10] = (n*(n*(CT(10)*n-CT(6))-CT(10))+CT(9))/CT(384);
0577             c[11] = ((-CT(77)*n-CT(8))*n+CT(42))/CT(3072);
0578             c[12] = (n*((CT(20)-CT(7)*n)*n-CT(28))+CT(14))/CT(1024);
0579             c[13] = ((-CT(7)*n-CT(40))*n+CT(28))/CT(2048);
0580             c[14] = (n*(CT(75)*n-CT(90))+CT(42))/CT(5120);
0581             break;
0582         case 7:
0583             c[0] = (CT(1)-n)/CT(4);
0584             c[1] = (CT(1)-n2)/CT(8);
0585             c[2] = (n*((-CT(5)*n-CT(1))*n+CT(3))+CT(3))/CT(64);
0586             c[3] = (n*((CT(2)-CT(2)*n)*n+CT(2))+CT(5))/CT(128);
0587             c[4] = (n*(CT(3)*n+CT(11))+CT(12))/CT(512);
0588             c[5] = (CT(10)*n+CT(21))/CT(1024);
0589             c[6] = ((n-CT(3))*n+CT(2))/CT(32);
0590             c[7] = (n*(n*(CT(2)*n-CT(3))-CT(2))+CT(3))/CT(64);
0591             c[8] = (n*((-CT(6)*n-CT(9))*n+CT(2))+CT(6))/CT(256);
0592             c[9] = ((CT(1)-CT(2)*n)*n+CT(5))/CT(256);
0593             c[10] = (CT(69)*n+CT(108))/CT(8192);
0594             c[11] = (n*((CT(5)-n)*n-CT(9))+CT(5))/CT(192);
0595             c[12] = (n*(n*(CT(10)*n-CT(6))-CT(10))+CT(9))/CT(384);
0596             c[13] = ((-CT(77)*n-CT(8))*n+CT(42))/CT(3072);
0597             c[14] = (CT(12)-n)/CT(1024);
0598             c[15] = (n*((CT(20)-CT(7)*n)*n-CT(28))+CT(14))/CT(1024);
0599             c[16] = ((-CT(7)*n-CT(40))*n+CT(28))/CT(2048);
0600             c[17] = (CT(72)-CT(43)*n)/CT(8192);
0601             c[18] = (n*(CT(75)*n-CT(90))+CT(42))/CT(5120);
0602             c[19] = (CT(9)-CT(15)*n)/CT(1024);
0603             c[20] = (CT(44)-CT(99)*n)/CT(8192);
0604             break;
0605         default:
0606             c[0] = (CT(1)-n)/CT(4);
0607             c[1] = (CT(1)-n2)/CT(8);
0608             c[2] = (n*((-CT(5)*n-CT(1))*n+CT(3))+CT(3))/CT(64);
0609             c[3] = (n*((CT(2)-CT(2)*n)*n+CT(2))+CT(5))/CT(128);
0610             c[4] = (n*(CT(3)*n+CT(11))+CT(12))/CT(512);
0611             c[5] = (CT(10)*n+CT(21))/CT(1024);
0612             c[6] = CT(243)/CT(16384);
0613             c[7] = ((n-CT(3))*n+CT(2))/CT(32);
0614             c[8] = (n*(n*(CT(2)*n-CT(3))-CT(2))+CT(3))/CT(64);
0615             c[9] = (n*((-CT(6)*n-CT(9))*n+CT(2))+CT(6))/CT(256);
0616             c[10] = ((CT(1)-CT(2)*n)*n+CT(5))/CT(256);
0617             c[11] = (CT(69)*n+CT(108))/CT(8192);
0618             c[12] = CT(187)/CT(16384);
0619             c[13] = (n*((CT(5)-n)*n-CT(9))+CT(5))/CT(192);
0620             c[14] = (n*(n*(CT(10)*n-CT(6))-CT(10))+CT(9))/CT(384);
0621             c[15] = ((-CT(77)*n-CT(8))*n+CT(42))/CT(3072);
0622             c[16] = (CT(12)-n)/CT(1024);
0623             c[17] = CT(139)/CT(16384);
0624             c[18] = (n*((CT(20)-CT(7)*n)*n-CT(28))+CT(14))/CT(1024);
0625             c[19] = ((-CT(7)*n-CT(40))*n+CT(28))/CT(2048);
0626             c[20] = (CT(72)-CT(43)*n)/CT(8192);
0627             c[21] = CT(127)/CT(16384);
0628             c[22] = (n*(CT(75)*n-CT(90))+CT(42))/CT(5120);
0629             c[23] = (CT(9)-CT(15)*n)/CT(1024);
0630             c[24] = CT(99)/CT(16384);
0631             c[25] = (CT(44)-CT(99)*n)/CT(8192);
0632             c[26] = CT(99)/CT(16384);
0633             c[27] = CT(429)/CT(114688);
0634             break;
0635         }
0636     }
0637 
0638     /*
0639     \brief Given the set of coefficients coeffs2[] evaluate on
0640       C3 and return the set of coefficients coeffs1[].
0641 
0642       Elements coeffs1[1] through coeffs1[SeriesOrder - 1] are set.
0643     */
0644     template <typename Coeffs1, typename Coeffs2, typename CT>
0645     inline void evaluate_coeffs_C3(Coeffs1 &coeffs1, Coeffs2 &coeffs2, CT const& eps)
0646     {
0647         CT mult = 1;
0648         size_t offset = 0;
0649 
0650         // i is the index of C3[i].
0651         for (size_t i = 1; i < Coeffs1::static_size; ++i)
0652         {
0653             // Order of polynomial in eps.
0654             size_t m = Coeffs1::static_size - i;
0655             mult *= eps;
0656 
0657             coeffs1[i] = mult * math::horner_evaluate(eps, coeffs2.begin() + offset,
0658                                                            coeffs2.begin() + offset + m);
0659 
0660             offset += m;
0661         }
0662         // Post condition: offset == coeffs_C3_size
0663     }
0664 
0665     /*
0666     \brief Evaluate the following:
0667 
0668      y = sum(c[i] * sin(2*i * x), i, 1, n)
0669 
0670      using Clenshaw summation.
0671     */
0672     template <typename CT, typename Coeffs>
0673     inline CT sin_cos_series(CT const& sinx, CT const& cosx, Coeffs const& coeffs)
0674     {
0675         size_t n = Coeffs::static_size - 1;
0676         size_t index = 0;
0677 
0678         // Point to one beyond last element.
0679         index += (n + 1);
0680         CT ar = 2 * (cosx - sinx) * (cosx + sinx);
0681 
0682         // If n is odd, get the last element.
0683         CT k0 = n & 1 ? coeffs[--index] : 0;
0684         CT k1 = 0;
0685 
0686         // Make n even.
0687         n /= 2;
0688         while (n--) {
0689             // Unroll loop x 2, so accumulators return to their original role.
0690             k1 = ar * k0 - k1 + coeffs[--index];
0691             k0 = ar * k1 - k0 + coeffs[--index];
0692         }
0693 
0694         return 2 * sinx * cosx * k0;
0695     }
0696 
0697     /*
0698      The coefficient containers for the series expansions.
0699      These structs allow the caller to only know the series order.
0700     */
0701     template <size_t SeriesOrder, typename CT>
0702     struct coeffs_C1 : boost::array<CT, SeriesOrder + 1>
0703     {
0704         coeffs_C1(CT const& epsilon)
0705         {
0706             evaluate_coeffs_C1(*this, epsilon);
0707         }
0708     };
0709 
0710     template <size_t SeriesOrder, typename CT>
0711     struct coeffs_C1p : boost::array<CT, SeriesOrder + 1>
0712     {
0713         coeffs_C1p(CT const& epsilon)
0714         {
0715             evaluate_coeffs_C1p(*this, epsilon);
0716         }
0717     };
0718 
0719     template <size_t SeriesOrder, typename CT>
0720     struct coeffs_C2 : boost::array<CT, SeriesOrder + 1>
0721     {
0722         coeffs_C2(CT const& epsilon)
0723         {
0724             evaluate_coeffs_C2(*this, epsilon);
0725         }
0726     };
0727 
0728     template <size_t SeriesOrder, typename CT>
0729     struct coeffs_C3x : boost::array<CT, (SeriesOrder * (SeriesOrder - 1)) / 2>
0730     {
0731         coeffs_C3x(CT const& n)
0732         {
0733             evaluate_coeffs_C3x<SeriesOrder>(*this, n);
0734         }
0735     };
0736 
0737     template <size_t SeriesOrder, typename CT>
0738     struct coeffs_C3 : boost::array<CT, SeriesOrder>
0739     {
0740         coeffs_C3(CT const& n, CT const& epsilon)
0741         {
0742             coeffs_C3x<SeriesOrder, CT> coeffs_C3x(n);
0743 
0744             evaluate_coeffs_C3(*this, coeffs_C3x, epsilon);
0745         }
0746     };
0747 
0748     template <size_t SeriesOrder, typename CT>
0749     struct coeffs_A3 : boost::array<CT, SeriesOrder>
0750     {
0751         coeffs_A3(CT const& n)
0752         {
0753             evaluate_coeffs_A3(*this, n);
0754         }
0755     };
0756 
0757 }}} // namespace boost::geometry::series_expansion
0758 
0759 #endif // BOOST_GEOMETRY_UTIL_SERIES_EXPANSION_HPP
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!