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