ceph/src/boost/boost/geometry/util/series_expansion.hpp

   1 // Boost.Geometry
   2
   3 // Copyright (c) 2018 Adeel Ahmad, Islamabad, Pakistan.
   4
   5 // Contributed and/or modified by Adeel Ahmad, as part of Google Summer of Code 2018 program.
   6
   7 // This file was modified by Oracle on 2019.
   8 // Modifications copyright (c) 2019 Oracle and/or its affiliates.
   9
  10 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  11
  12 // Use, modification and distribution is subject to the Boost Software License,
  13 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  14 // http://www.boost.org/LICENSE_1_0.txt)
  15
  16 // This file is converted from GeographicLib, https://geographiclib.sourceforge.io
  17 // GeographicLib is originally written by Charles Karney.
  18
  19 // Author: Charles Karney (2008-2017)
  20
  21 // Last updated version of GeographicLib: 1.49
  22
  23 // Original copyright notice:
  24
  25 // Copyright (c) Charles Karney (2008-2017) <charles@karney.com> and licensed
  26 // under the MIT/X11 License. For more information, see
  27 // https://geographiclib.sourceforge.io
  28
  29 #ifndef BOOST_GEOMETRY_UTIL_SERIES_EXPANSION_HPP
  30 #define BOOST_GEOMETRY_UTIL_SERIES_EXPANSION_HPP
  31
  32 #include <boost/geometry/core/assert.hpp>
  33 #include <boost/geometry/util/math.hpp>
  34
  35 namespace boost { namespace geometry { namespace series_expansion {
  36
  37     /*
  38      Generate and evaluate the series expansion of the following integral
  39
  40      I1 = integrate( sqrt(1+k2*sin(sigma1)^2), sigma1, 0, sigma )
  41
  42      which is valid for k2 small. We substitute k2 = 4 * eps / (1 - eps)^2
  43      and expand (1 - eps) * I1 retaining terms up to order eps^maxpow
  44      in A1 and C1[l].
  45
  46      The resulting series is of the form
  47
  48      A1 * ( sigma + sum(C1[l] * sin(2*l*sigma), l, 1, maxpow) ).
  49
  50      The scale factor A1-1 = mean value of (d/dsigma)I1 - 1
  51
  52      The expansion above is performed in Maxima, a Computer Algebra System.
  53      The C++ code (that yields the function evaluate_A1 below) is
  54      generated by the following Maxima script:
  55      geometry/doc/other/maxima/geod.mac
  56
  57      To replace each number x by CT(x) the following
  58      script can be used:
  59        sed -e 's/[0-9]\+/CT(&)/g; s/\[CT/\[/g; s/)\]/\]/g;
  60                s/case\sCT(/case /g; s/):/:/g; s/epsCT(2)/eps2/g;'
  61     */
  62     template <size_t SeriesOrder, typename CT>
  63     inline CT evaluate_A1(CT eps)
  64     {
  65         CT eps2 = math::sqr(eps);
  66         CT t;
  67         switch (SeriesOrder/2)
  68         {
  69         case 0:
  70             t = CT(0);
  71             break;
  72         case 1:
  73             t = eps2/CT(4);
  74             break;
  75         case 2:
  76             t = eps2*(eps2+CT(16))/CT(64);
  77             break;
  78         case 3:
  79             t = eps2*(eps2*(eps2+CT(4))+CT(64))/CT(256);
  80             break;
  81         case 4:
  82             t = eps2*(eps2*(eps2*(CT(25)*eps2+CT(64))+CT(256))+CT(4096))/CT(16384);
  83             break;
  84         }
  85         return (t + eps) / (CT(1) - eps);
  86     }
  87
  88     /*
  89      Generate and evaluate the series expansion of the following integral
  90
  91      I2 = integrate( 1/sqrt(1+k2*sin(sigma1)^2), sigma1, 0, sigma )
  92
  93      which is valid for k2 small. We substitute k2 = 4 * eps / (1 - eps)^2
  94      and expand (1 - eps) * I2 retaining terms up to order eps^maxpow
  95      in A2 and C2[l].
  96
  97      The resulting series is of the form
  98
  99      A2 * ( sigma + sum(C2[l] * sin(2*l*sigma), l, 1, maxpow) )
 100
 101      The scale factor A2-1 = mean value of (d/dsigma)2 - 1
 102
 103      The expansion above is performed in Maxima, a Computer Algebra System.
 104      The C++ code (that yields the function evaluate_A2 below) is
 105      generated by the following Maxima script:
 106      geometry/doc/other/maxima/geod.mac
 107     */
 108     template <size_t SeriesOrder, typename CT>
 109     inline CT evaluate_A2(CT const& eps)
 110     {
 111         CT const eps2 = math::sqr(eps);
 112         CT t;
 113         switch (SeriesOrder/2)
 114         {
 115         case 0:
 116             t = CT(0);
 117             break;
 118         case 1:
 119             t = -CT(3)*eps2/CT(4);
 120             break;
 121         case 2:
 122             t = (-CT(7)*eps2-CT(48))*eps2/CT(64);
 123             break;
 124         case 3:
 125             t = eps2*((-CT(11)*eps2-CT(28))*eps2-CT(192))/CT(256);
 126             break;
 127         default:
 128             t = eps2*(eps2*((-CT(375)*eps2-CT(704))*eps2-CT(1792))-CT(12288))/CT(16384);
 129             break;
 130         }
 131         return (t - eps) / (CT(1) + eps);
 132     }
 133
 134     /*
 135      Express
 136
 137         I3 = integrate( (2-f)/(1+(1-f)*sqrt(1+k2*sin(sigma1)^2)), sigma1, 0, sigma )
 138
 139      as a series
 140
 141         A3 * ( sigma + sum(C3[l] * sin(2*l*sigma), l, 1, maxpow-1) )
 142
 143      valid for f and k2 small.  It is convenient to write k2 = 4 * eps / (1 -
 144      eps)^2 and f = 2*n/(1+n) and expand in eps and n.  This procedure leads
 145      to a series where the coefficients of eps^j are terminating series in n.
 146
 147      The scale factor A3 = mean value of (d/dsigma)I3
 148
 149      The expansion above is performed in Maxima, a Computer Algebra System.
 150      The C++ code (that yields the function evaluate_coeffs_A3 below) is
 151      generated by the following Maxima script:
 152      geometry/doc/other/maxima/geod.mac
 153     */
 154     template <typename Coeffs, typename CT>
 155     inline void evaluate_coeffs_A3(Coeffs &c, CT const& n)
 156     {
 157         switch (int(Coeffs::static_size))
 158         {
 159         case 0:
 160             break;
 161         case 1:
 162             c[0] = CT(1);
 163             break;
 164         case 2:
 165             c[0] = CT(1);
 166             c[1] = -CT(1)/CT(2);
 167             break;
 168         case 3:
 169             c[0] = CT(1);
 170             c[1] = (n-CT(1))/CT(2);
 171             c[2] = -CT(1)/CT(4);
 172             break;
 173         case 4:
 174             c[0] = CT(1);
 175             c[1] = (n-CT(1))/CT(2);
 176             c[2] = (-n-CT(2))/CT(8);
 177             c[3] = -CT(1)/CT(16);
 178             break;
 179         case 5:
 180             c[0] = CT(1);
 181             c[1] = (n-CT(1))/CT(2);
 182             c[2] = (n*(CT(3)*n-CT(1))-CT(2))/CT(8);
 183             c[3] = (-CT(3)*n-CT(1))/CT(16);
 184             c[4] = -CT(3)/CT(64);
 185             break;
 186         case 6:
 187             c[0] = CT(1);
 188             c[1] = (n-CT(1))/CT(2);
 189             c[2] = (n*(CT(3)*n-CT(1))-CT(2))/CT(8);
 190             c[3] = ((-n-CT(3))*n-CT(1))/CT(16);
 191             c[4] = (-CT(2)*n-CT(3))/CT(64);
 192             c[5] = -CT(3)/CT(128);
 193             break;
 194         case 7:
 195             c[0] = CT(1);
 196             c[1] = (n-CT(1))/CT(2);
 197             c[2] = (n*(CT(3)*n-CT(1))-CT(2))/CT(8);
 198             c[3] = (n*(n*(CT(5)*n-CT(1))-CT(3))-CT(1))/CT(16);
 199             c[4] = ((-CT(10)*n-CT(2))*n-CT(3))/CT(64);
 200             c[5] = (-CT(5)*n-CT(3))/CT(128);
 201             c[6] = -CT(5)/CT(256);
 202             break;
 203         default:
 204             c[0] = CT(1);
 205             c[1] = (n-CT(1))/CT(2);
 206             c[2] = (n*(CT(3)*n-CT(1))-CT(2))/CT(8);
 207             c[3] = (n*(n*(CT(5)*n-CT(1))-CT(3))-CT(1))/CT(16);
 208             c[4] = (n*((-CT(5)*n-CT(20))*n-CT(4))-CT(6))/CT(128);
 209             c[5] = ((-CT(5)*n-CT(10))*n-CT(6))/CT(256);
 210             c[6] = (-CT(15)*n-CT(20))/CT(1024);
 211             c[7] = -CT(25)/CT(2048);
 212             break;
 213         }
 214     }
 215
 216     /*
 217      The coefficients C1[l] in the Fourier expansion of B1.
 218
 219      The expansion below is performed in Maxima, a Computer Algebra System.
 220      The C++ code (that yields the function evaluate_coeffs_C1 below) is
 221      generated by the following Maxima script:
 222      geometry/doc/other/maxima/geod.mac
 223     */
 224     template <typename Coeffs, typename CT>
 225     inline void evaluate_coeffs_C1(Coeffs &c, CT const& eps)
 226     {
 227         CT eps2 = math::sqr(eps);
 228         CT d = eps;
 229         switch (int(Coeffs::static_size) - 1)
 230         {
 231         case 0:
 232             break;
 233         case 1:
 234             c[1] = -d/CT(2);
 235             break;
 236         case 2:
 237             c[1] = -d/CT(2);
 238             d *= eps;
 239             c[2] = -d/CT(16);
 240             break;
 241         case 3:
 242             c[1] = d*(CT(3)*eps2-CT(8))/CT(16);
 243             d *= eps;
 244             c[2] = -d/CT(16);
 245             d *= eps;
 246             c[3] = -d/CT(48);
 247             break;
 248         case 4:
 249             c[1] = d*(CT(3)*eps2-CT(8))/CT(16);
 250             d *= eps;
 251             c[2] = d*(eps2-CT(2))/CT(32);
 252             d *= eps;
 253             c[3] = -d/CT(48);
 254             d *= eps;
 255             c[4] = -CT(5)*d/CT(512);
 256             break;
 257         case 5:
 258             c[1] = d*((CT(6)-eps2)*eps2-CT(16))/CT(32);
 259             d *= eps;
 260             c[2] = d*(eps2-CT(2))/CT(32);
 261             d *= eps;
 262             c[3] = d*(CT(9)*eps2-CT(16))/CT(768);
 263             d *= eps;
 264             c[4] = -CT(5)*d/CT(512);
 265             d *= eps;
 266             c[5] = -CT(7)*d/CT(1280);
 267             break;
 268         case 6:
 269             c[1] = d*((CT(6)-eps2)*eps2-CT(16))/CT(32);
 270             d *= eps;
 271             c[2] = d*((CT(64)-CT(9)*eps2)*eps2-CT(128))/CT(2048);
 272             d *= eps;
 273             c[3] = d*(CT(9)*eps2-CT(16))/CT(768);
 274             d *= eps;
 275             c[4] = d*(CT(3)*eps2-CT(5))/CT(512);
 276             d *= eps;
 277             c[5] = -CT(7)*d/CT(1280);
 278             d *= eps;
 279             c[6] = -CT(7)*d/CT(2048);
 280             break;
 281         case 7:
 282             c[1] = d*(eps2*(eps2*(CT(19)*eps2-CT(64))+CT(384))-CT(1024))/CT(2048);
 283             d *= eps;
 284             c[2] = d*((CT(64)-CT(9)*eps2)*eps2-CT(128))/CT(2048);
 285             d *= eps;
 286             c[3] = d*((CT(72)-CT(9)*eps2)*eps2-CT(128))/CT(6144);
 287             d *= eps;
 288             c[4] = d*(CT(3)*eps2-CT(5))/CT(512);
 289             d *= eps;
 290             c[5] = d*(CT(35)*eps2-CT(56))/CT(10240);
 291             d *= eps;
 292             c[6] = -CT(7)*d/CT(2048);
 293             d *= eps;
 294             c[7] = -CT(33)*d/CT(14336);
 295             break;
 296         default:
 297             c[1] = d*(eps2*(eps2*(CT(19)*eps2-CT(64))+CT(384))-CT(1024))/CT(2048);
 298             d *= eps;
 299             c[2] = d*(eps2*(eps2*(CT(7)*eps2-CT(18))+CT(128))-CT(256))/CT(4096);
 300             d *= eps;
 301             c[3] = d*((CT(72)-CT(9)*eps2)*eps2-CT(128))/CT(6144);
 302             d *= eps;
 303             c[4] = d*((CT(96)-CT(11)*eps2)*eps2-CT(160))/CT(16384);
 304             d *= eps;
 305             c[5] = d*(CT(35)*eps2-CT(56))/CT(10240);
 306             d *= eps;
 307             c[6] = d*(CT(9)*eps2-CT(14))/CT(4096);
 308             d *= eps;
 309             c[7] = -CT(33)*d/CT(14336);
 310             d *= eps;
 311             c[8] = -CT(429)*d/CT(262144);
 312             break;
 313         }
 314     }
 315
 316     /*
 317      The coefficients C1p[l] in the Fourier expansion of B1p.
 318
 319      The expansion below is performed in Maxima, a Computer Algebra System.
 320      The C++ code (that yields the function evaluate_coeffs_C1p below) is
 321      generated by the following Maxima script:
 322      geometry/doc/other/maxima/geod.mac
 323     */
 324     template <typename Coeffs, typename CT>
 325     inline void evaluate_coeffs_C1p(Coeffs& c, CT const& eps)
 326     {
 327         CT const eps2 = math::sqr(eps);
 328         CT d = eps;
 329         switch (int(Coeffs::static_size) - 1)
 330         {
 331         case 0:
 332             break;
 333         case 1:
 334             c[1] = d/CT(2);
 335             break;
 336         case 2:
 337             c[1] = d/CT(2);
 338             d *= eps;
 339             c[2] = CT(5)*d/CT(16);
 340             break;
 341         case 3:
 342             c[1] = d*(CT(16)-CT(9)*eps2)/CT(32);
 343             d *= eps;
 344             c[2] = CT(5)*d/CT(16);
 345             d *= eps;
 346             c[3] = CT(29)*d/CT(96);
 347             break;
 348         case 4:
 349             c[1] = d*(CT(16)-CT(9)*eps2)/CT(32);
 350             d *= eps;
 351             c[2] = d*(CT(30)-CT(37)*eps2)/CT(96);
 352             d *= eps;
 353             c[3] = CT(29)*d/CT(96);
 354             d *= eps;
 355             c[4] = CT(539)*d/CT(1536);
 356             break;
 357         case 5:
 358             c[1] = d*(eps2*(CT(205)*eps2-CT(432))+CT(768))/CT(1536);
 359             d *= eps;
 360             c[2] = d*(CT(30)-CT(37)*eps2)/CT(96);
 361             d *= eps;
 362             c[3] = d*(CT(116)-CT(225)*eps2)/CT(384);
 363             d *= eps;
 364             c[4] = CT(539)*d/CT(1536);
 365             d *= eps;
 366             c[5] = CT(3467)*d/CT(7680);
 367             break;
 368         case 6:
 369             c[1] = d*(eps2*(CT(205)*eps2-CT(432))+CT(768))/CT(1536);
 370             d *= eps;
 371             c[2] = d*(eps2*(CT(4005)*eps2-CT(4736))+CT(3840))/CT(12288);
 372             d *= eps;
 373             c[3] = d*(CT(116)-CT(225)*eps2)/CT(384);
 374             d *= eps;
 375             c[4] = d*(CT(2695)-CT(7173)*eps2)/CT(7680);
 376             d *= eps;
 377             c[5] = CT(3467)*d/CT(7680);
 378             d *= eps;
 379             c[6] = CT(38081)*d/CT(61440);
 380             break;
 381         case 7:
 382             c[1] = d*(eps2*((CT(9840)-CT(4879)*eps2)*eps2-CT(20736))+CT(36864))/CT(73728);
 383             d *= eps;
 384             c[2] = d*(eps2*(CT(4005)*eps2-CT(4736))+CT(3840))/CT(12288);
 385             d *= eps;
 386             c[3] = d*(eps2*(CT(8703)*eps2-CT(7200))+CT(3712))/CT(12288);
 387             d *= eps;
 388             c[4] = d*(CT(2695)-CT(7173)*eps2)/CT(7680);
 389             d *= eps;
 390             c[5] = d*(CT(41604)-CT(141115)*eps2)/CT(92160);
 391             d *= eps;
 392             c[6] = CT(38081)*d/CT(61440);
 393             d *= eps;
 394             c[7] = CT(459485)*d/CT(516096);
 395             break;
 396         default:
 397             c[1] = d*(eps2*((CT(9840)-CT(4879)*eps2)*eps2-CT(20736))+CT(36864))/CT(73728);
 398             d *= eps;
 399             c[2] = d*(eps2*((CT(120150)-CT(86171)*eps2)*eps2-CT(142080))+CT(115200))/CT(368640);
 400             d *= eps;
 401             c[3] = d*(eps2*(CT(8703)*eps2-CT(7200))+CT(3712))/CT(12288);
 402             d *= eps;
 403             c[4] = d*(eps2*(CT(1082857)*eps2-CT(688608))+CT(258720))/CT(737280);
 404             d *= eps;
 405             c[5] = d*(CT(41604)-CT(141115)*eps2)/CT(92160);
 406             d *= eps;
 407             c[6] = d*(CT(533134)-CT(2200311)*eps2)/CT(860160);
 408             d *= eps;
 409             c[7] = CT(459485)*d/CT(516096);
 410             d *= eps;
 411             c[8] = CT(109167851)*d/CT(82575360);
 412             break;
 413         }
 414     }
 415
 416     /*
 417      The coefficients C2[l] in the Fourier expansion of B2.
 418
 419      The expansion below is performed in Maxima, a Computer Algebra System.
 420      The C++ code (that yields the function evaluate_coeffs_C2 below) is
 421      generated by the following Maxima script:
 422      geometry/doc/other/maxima/geod.mac
 423     */
 424     template <typename Coeffs, typename CT>
 425     inline void evaluate_coeffs_C2(Coeffs& c, CT const& eps)
 426     {
 427         CT const eps2 = math::sqr(eps);
 428         CT d = eps;
 429         switch (int(Coeffs::static_size) - 1)
 430         {
 431         case 0:
 432             break;
 433         case 1:
 434             c[1] = d/CT(2);
 435             break;
 436         case 2:
 437             c[1] = d/CT(2);
 438             d *= eps;
 439             c[2] = CT(3)*d/CT(16);
 440             break;
 441         case 3:
 442             c[1] = d*(eps2+CT(8))/CT(16);
 443             d *= eps;
 444             c[2] = CT(3)*d/CT(16);
 445             d *= eps;
 446             c[3] = CT(5)*d/CT(48);
 447             break;
 448         case 4:
 449             c[1] = d*(eps2+CT(8))/CT(16);
 450             d *= eps;
 451             c[2] = d*(eps2+CT(6))/CT(32);
 452             d *= eps;
 453             c[3] = CT(5)*d/CT(48);
 454             d *= eps;
 455             c[4] = CT(35)*d/CT(512);
 456             break;
 457         case 5:
 458             c[1] = d*(eps2*(eps2+CT(2))+CT(16))/CT(32);
 459             d *= eps;
 460             c[2] = d*(eps2+CT(6))/CT(32);
 461             d *= eps;
 462             c[3] = d*(CT(15)*eps2+CT(80))/CT(768);
 463             d *= eps;
 464             c[4] = CT(35)*d/CT(512);
 465             d *= eps;
 466             c[5] = CT(63)*d/CT(1280);
 467             break;
 468         case 6:
 469             c[1] = d*(eps2*(eps2+CT(2))+CT(16))/CT(32);
 470             d *= eps;
 471             c[2] = d*(eps2*(CT(35)*eps2+CT(64))+CT(384))/CT(2048);
 472             d *= eps;
 473             c[3] = d*(CT(15)*eps2+CT(80))/CT(768);
 474             d *= eps;
 475             c[4] = d*(CT(7)*eps2+CT(35))/CT(512);
 476             d *= eps;
 477             c[5] = CT(63)*d/CT(1280);
 478             d *= eps;
 479             c[6] = CT(77)*d/CT(2048);
 480             break;
 481         case 7:
 482             c[1] = d*(eps2*(eps2*(CT(41)*eps2+CT(64))+CT(128))+CT(1024))/CT(2048);
 483             d *= eps;
 484             c[2] = d*(eps2*(CT(35)*eps2+CT(64))+CT(384))/CT(2048);
 485             d *= eps;
 486             c[3] = d*(eps2*(CT(69)*eps2+CT(120))+CT(640))/CT(6144);
 487             d *= eps;
 488             c[4] = d*(CT(7)*eps2+CT(35))/CT(512);
 489             d *= eps;
 490             c[5] = d*(CT(105)*eps2+CT(504))/CT(10240);
 491             d *= eps;
 492             c[6] = CT(77)*d/CT(2048);
 493             d *= eps;
 494             c[7] = CT(429)*d/CT(14336);
 495             break;
 496         default:
 497             c[1] = d*(eps2*(eps2*(CT(41)*eps2+CT(64))+CT(128))+CT(1024))/CT(2048);
 498             d *= eps;
 499             c[2] = d*(eps2*(eps2*(CT(47)*eps2+CT(70))+CT(128))+CT(768))/CT(4096);
 500             d *= eps;
 501             c[3] = d*(eps2*(CT(69)*eps2+CT(120))+CT(640))/CT(6144);
 502             d *= eps;
 503             c[4] = d*(eps2*(CT(133)*eps2+CT(224))+CT(1120))/CT(16384);
 504             d *= eps;
 505             c[5] = d*(CT(105)*eps2+CT(504))/CT(10240);
 506             d *= eps;
 507             c[6] = d*(CT(33)*eps2+CT(154))/CT(4096);
 508             d *= eps;
 509             c[7] = CT(429)*d/CT(14336);
 510             d *= eps;
 511             c[8] = CT(6435)*d/CT(262144);
 512             break;
 513         }
 514     }
 515
 516     /*
 517      The coefficients C3[l] in the Fourier expansion of B3.
 518
 519      The expansion below is performed in Maxima, a Computer Algebra System.
 520      The C++ code (that yields the function evaluate_coeffs_C3 below) is
 521      generated by the following Maxima script:
 522      geometry/doc/other/maxima/geod.mac
 523     */
 524     template <size_t SeriesOrder, typename Coeffs, typename CT>
 525     inline void evaluate_coeffs_C3x(Coeffs &c, CT const& n) {
 526         BOOST_GEOMETRY_ASSERT((Coeffs::static_size == (SeriesOrder * (SeriesOrder - 1)) / 2));
 527
 528         CT const n2 = math::sqr(n);
 529         switch (SeriesOrder)
 530         {
 531         case 0:
 532             break;
 533         case 1:
 534             break;
 535         case 2:
 536             c[0] = (CT(1)-n)/CT(4);
 537             break;
 538         case 3:
 539             c[0] = (CT(1)-n)/CT(4);
 540             c[1] = (CT(1)-n2)/CT(8);
 541             c[2] = ((n-CT(3))*n+CT(2))/CT(32);
 542             break;
 543         case 4:
 544             c[0] = (CT(1)-n)/CT(4);
 545             c[1] = (CT(1)-n2)/CT(8);
 546             c[2] = (n*((-CT(5)*n-CT(1))*n+CT(3))+CT(3))/CT(64);
 547             c[3] = ((n-CT(3))*n+CT(2))/CT(32);
 548             c[4] = (n*(n*(CT(2)*n-CT(3))-CT(2))+CT(3))/CT(64);
 549             c[5] = (n*((CT(5)-n)*n-CT(9))+CT(5))/CT(192);
 550             break;
 551         case 5:
 552             c[0] = (CT(1)-n)/CT(4);
 553             c[1] = (CT(1)-n2)/CT(8);
 554             c[2] = (n*((-CT(5)*n-CT(1))*n+CT(3))+CT(3))/CT(64);
 555             c[3] = (n*((CT(2)-CT(2)*n)*n+CT(2))+CT(5))/CT(128);
 556             c[4] = ((n-CT(3))*n+CT(2))/CT(32);
 557             c[5] = (n*(n*(CT(2)*n-CT(3))-CT(2))+CT(3))/CT(64);
 558             c[6] = (n*((-CT(6)*n-CT(9))*n+CT(2))+CT(6))/CT(256);
 559             c[7] = (n*((CT(5)-n)*n-CT(9))+CT(5))/CT(192);
 560             c[8] = (n*(n*(CT(10)*n-CT(6))-CT(10))+CT(9))/CT(384);
 561             c[9] = (n*((CT(20)-CT(7)*n)*n-CT(28))+CT(14))/CT(1024);
 562             break;
 563         case 6:
 564             c[0] = (CT(1)-n)/CT(4);
 565             c[1] = (CT(1)-n2)/CT(8);
 566             c[2] = (n*((-CT(5)*n-CT(1))*n+CT(3))+CT(3))/CT(64);
 567             c[3] = (n*((CT(2)-CT(2)*n)*n+CT(2))+CT(5))/CT(128);
 568             c[4] = (n*(CT(3)*n+CT(11))+CT(12))/CT(512);
 569             c[5] = ((n-CT(3))*n+CT(2))/CT(32);
 570             c[6] = (n*(n*(CT(2)*n-CT(3))-CT(2))+CT(3))/CT(64);
 571             c[7] = (n*((-CT(6)*n-CT(9))*n+CT(2))+CT(6))/CT(256);
 572             c[8] = ((CT(1)-CT(2)*n)*n+CT(5))/CT(256);
 573             c[9] = (n*((CT(5)-n)*n-CT(9))+CT(5))/CT(192);
 574             c[10] = (n*(n*(CT(10)*n-CT(6))-CT(10))+CT(9))/CT(384);
 575             c[11] = ((-CT(77)*n-CT(8))*n+CT(42))/CT(3072);
 576             c[12] = (n*((CT(20)-CT(7)*n)*n-CT(28))+CT(14))/CT(1024);
 577             c[13] = ((-CT(7)*n-CT(40))*n+CT(28))/CT(2048);
 578             c[14] = (n*(CT(75)*n-CT(90))+CT(42))/CT(5120);
 579             break;
 580         case 7:
 581             c[0] = (CT(1)-n)/CT(4);
 582             c[1] = (CT(1)-n2)/CT(8);
 583             c[2] = (n*((-CT(5)*n-CT(1))*n+CT(3))+CT(3))/CT(64);
 584             c[3] = (n*((CT(2)-CT(2)*n)*n+CT(2))+CT(5))/CT(128);
 585             c[4] = (n*(CT(3)*n+CT(11))+CT(12))/CT(512);
 586             c[5] = (CT(10)*n+CT(21))/CT(1024);
 587             c[6] = ((n-CT(3))*n+CT(2))/CT(32);
 588             c[7] = (n*(n*(CT(2)*n-CT(3))-CT(2))+CT(3))/CT(64);
 589             c[8] = (n*((-CT(6)*n-CT(9))*n+CT(2))+CT(6))/CT(256);
 590             c[9] = ((CT(1)-CT(2)*n)*n+CT(5))/CT(256);
 591             c[10] = (CT(69)*n+CT(108))/CT(8192);
 592             c[11] = (n*((CT(5)-n)*n-CT(9))+CT(5))/CT(192);
 593             c[12] = (n*(n*(CT(10)*n-CT(6))-CT(10))+CT(9))/CT(384);
 594             c[13] = ((-CT(77)*n-CT(8))*n+CT(42))/CT(3072);
 595             c[14] = (CT(12)-n)/CT(1024);
 596             c[15] = (n*((CT(20)-CT(7)*n)*n-CT(28))+CT(14))/CT(1024);
 597             c[16] = ((-CT(7)*n-CT(40))*n+CT(28))/CT(2048);
 598             c[17] = (CT(72)-CT(43)*n)/CT(8192);
 599             c[18] = (n*(CT(75)*n-CT(90))+CT(42))/CT(5120);
 600             c[19] = (CT(9)-CT(15)*n)/CT(1024);
 601             c[20] = (CT(44)-CT(99)*n)/CT(8192);
 602             break;
 603         default:
 604             c[0] = (CT(1)-n)/CT(4);
 605             c[1] = (CT(1)-n2)/CT(8);
 606             c[2] = (n*((-CT(5)*n-CT(1))*n+CT(3))+CT(3))/CT(64);
 607             c[3] = (n*((CT(2)-CT(2)*n)*n+CT(2))+CT(5))/CT(128);
 608             c[4] = (n*(CT(3)*n+CT(11))+CT(12))/CT(512);
 609             c[5] = (CT(10)*n+CT(21))/CT(1024);
 610             c[6] = CT(243)/CT(16384);
 611             c[7] = ((n-CT(3))*n+CT(2))/CT(32);
 612             c[8] = (n*(n*(CT(2)*n-CT(3))-CT(2))+CT(3))/CT(64);
 613             c[9] = (n*((-CT(6)*n-CT(9))*n+CT(2))+CT(6))/CT(256);
 614             c[10] = ((CT(1)-CT(2)*n)*n+CT(5))/CT(256);
 615             c[11] = (CT(69)*n+CT(108))/CT(8192);
 616             c[12] = CT(187)/CT(16384);
 617             c[13] = (n*((CT(5)-n)*n-CT(9))+CT(5))/CT(192);
 618             c[14] = (n*(n*(CT(10)*n-CT(6))-CT(10))+CT(9))/CT(384);
 619             c[15] = ((-CT(77)*n-CT(8))*n+CT(42))/CT(3072);
 620             c[16] = (CT(12)-n)/CT(1024);
 621             c[17] = CT(139)/CT(16384);
 622             c[18] = (n*((CT(20)-CT(7)*n)*n-CT(28))+CT(14))/CT(1024);
 623             c[19] = ((-CT(7)*n-CT(40))*n+CT(28))/CT(2048);
 624             c[20] = (CT(72)-CT(43)*n)/CT(8192);
 625             c[21] = CT(127)/CT(16384);
 626             c[22] = (n*(CT(75)*n-CT(90))+CT(42))/CT(5120);
 627             c[23] = (CT(9)-CT(15)*n)/CT(1024);
 628             c[24] = CT(99)/CT(16384);
 629             c[25] = (CT(44)-CT(99)*n)/CT(8192);
 630             c[26] = CT(99)/CT(16384);
 631             c[27] = CT(429)/CT(114688);
 632             break;
 633         }
 634     }
 635
 636     /*
 637     \brief Given the set of coefficients coeffs2[] evaluate on
 638       C3 and return the set of coefficients coeffs1[].
 639
 640       Elements coeffs1[1] through coeffs1[SeriesOrder - 1] are set.
 641     */
 642     template <typename Coeffs1, typename Coeffs2, typename CT>
 643     inline void evaluate_coeffs_C3(Coeffs1 &coeffs1, Coeffs2 &coeffs2, CT const& eps)
 644     {
 645         CT mult = 1;
 646         size_t offset = 0;
 647
 648         // i is the index of C3[i].
 649         for (size_t i = 1; i < Coeffs1::static_size; ++i)
 650         {
 651             // Order of polynomial in eps.
 652             size_t m = Coeffs1::static_size - i;
 653             mult *= eps;
 654
 655             coeffs1[i] = mult * math::horner_evaluate(eps, coeffs2.begin() + offset,
 656                                                            coeffs2.begin() + offset + m);
 657
 658             offset += m;
 659         }
 660         // Post condition: offset == coeffs_C3_size
 661     }
 662
 663     /*
 664     \brief Evaluate the following:
 665
 666      y = sum(c[i] * sin(2*i * x), i, 1, n)
 667
 668      using Clenshaw summation.
 669     */
 670     template <typename CT, typename Coeffs>
 671     inline CT sin_cos_series(CT const& sinx, CT const& cosx, Coeffs const& coeffs)
 672     {
 673         size_t n = Coeffs::static_size - 1;
 674         size_t index = 0;
 675
 676         // Point to one beyond last element.
 677         index += (n + 1);
 678         CT ar = 2 * (cosx - sinx) * (cosx + sinx);
 679
 680         // If n is odd, get the last element.
 681         CT k0 = n & 1 ? coeffs[--index] : 0;
 682         CT k1 = 0;
 683
 684         // Make n even.
 685         n /= 2;
 686         while (n--) {
 687             // Unroll loop x 2, so accumulators return to their original role.
 688             k1 = ar * k0 - k1 + coeffs[--index];
 689             k0 = ar * k1 - k0 + coeffs[--index];
 690         }
 691
 692         return 2 * sinx * cosx * k0;
 693     }
 694
 695     /*
 696      The coefficient containers for the series expansions.
 697      These structs allow the caller to only know the series order.
 698     */
 699     template <size_t SeriesOrder, typename CT>
 700     struct coeffs_C1 : boost::array<CT, SeriesOrder + 1>
 701     {
 702         coeffs_C1(CT const& epsilon)
 703         {
 704             evaluate_coeffs_C1(*this, epsilon);
 705         }
 706     };
 707
 708     template <size_t SeriesOrder, typename CT>
 709     struct coeffs_C1p : boost::array<CT, SeriesOrder + 1>
 710     {
 711         coeffs_C1p(CT const& epsilon)
 712         {
 713             evaluate_coeffs_C1p(*this, epsilon);
 714         }
 715     };
 716
 717     template <size_t SeriesOrder, typename CT>
 718     struct coeffs_C2 : boost::array<CT, SeriesOrder + 1>
 719     {
 720         coeffs_C2(CT const& epsilon)
 721         {
 722             evaluate_coeffs_C2(*this, epsilon);
 723         }
 724     };
 725
 726     template <size_t SeriesOrder, typename CT>
 727     struct coeffs_C3x : boost::array<CT, (SeriesOrder * (SeriesOrder - 1)) / 2>
 728     {
 729         coeffs_C3x(CT const& n)
 730         {
 731             evaluate_coeffs_C3x<SeriesOrder>(*this, n);
 732         }
 733     };
 734
 735     template <size_t SeriesOrder, typename CT>
 736     struct coeffs_C3 : boost::array<CT, SeriesOrder>
 737     {
 738         coeffs_C3(CT const& n, CT const& epsilon)
 739         {
 740             coeffs_C3x<SeriesOrder, CT> coeffs_C3x(n);
 741
 742             evaluate_coeffs_C3(*this, coeffs_C3x, epsilon);
 743         }
 744     };
 745
 746     template <size_t SeriesOrder, typename CT>
 747     struct coeffs_A3 : boost::array<CT, SeriesOrder>
 748     {
 749         coeffs_A3(CT const& n)
 750         {
 751             evaluate_coeffs_A3(*this, n);
 752         }
 753     };
 754
 755 }}} // namespace boost::geometry::series_expansion
 756
 757 #endif // BOOST_GEOMETRY_UTIL_SERIES_EXPANSION_HPP