1 // (C) Copyright Nick Thompson 2019. 2 // Use, modification and distribution are subject to the 3 // Boost Software License, Version 1.0. (See accompanying file 4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 5 6 #ifndef BOOST_MATH_INTERPOLATORS_DETAIL_CARDINAL_TRIGONOMETRIC_HPP 7 #define BOOST_MATH_INTERPOLATORS_DETAIL_CARDINAL_TRIGONOMETRIC_HPP 8 #include <cmath> 9 #include <stdexcept> 10 #include <fftw3.h> 11 #include <boost/math/constants/constants.hpp> 12 13 #ifdef BOOST_HAS_FLOAT128 14 #include <quadmath.h> 15 #endif 16 17 namespace boost { namespace math { namespace interpolators { namespace detail { 18 19 template<typename Real> 20 class cardinal_trigonometric_detail { 21 public: cardinal_trigonometric_detail(const Real * data,size_t length,Real t0,Real h)22 cardinal_trigonometric_detail(const Real* data, size_t length, Real t0, Real h) 23 { 24 m_data = data; 25 m_length = length; 26 m_t0 = t0; 27 m_h = h; 28 throw std::domain_error("Not implemented."); 29 } 30 private: 31 size_t m_length; 32 Real m_t0; 33 Real m_h; 34 Real* m_data; 35 }; 36 37 template<> 38 class cardinal_trigonometric_detail<float> { 39 public: cardinal_trigonometric_detail(const float * data,size_t length,float t0,float h)40 cardinal_trigonometric_detail(const float* data, size_t length, float t0, float h) : m_t0{t0}, m_h{h} 41 { 42 if (length == 0) 43 { 44 throw std::logic_error("At least one sample is required."); 45 } 46 if (h <= 0) 47 { 48 throw std::logic_error("The step size must be > 0"); 49 } 50 // The period sadly must be stored, since the complex vector has length that cannot be used to recover the period: 51 m_T = m_h*length; 52 m_complex_vector_size = length/2 + 1; 53 m_gamma = fftwf_alloc_complex(m_complex_vector_size); 54 // The const_cast is legitimate: FFTW does not change the data as long as FFTW_ESTIMATE is provided. 55 fftwf_plan plan = fftwf_plan_dft_r2c_1d(length, const_cast<float*>(data), m_gamma, FFTW_ESTIMATE); 56 // FFTW says a null plan is impossible with the basic interface we are using, and I have no reason to doubt them. 57 // But it just feels weird not to check this: 58 if (!plan) 59 { 60 throw std::logic_error("A null fftw plan was created."); 61 } 62 63 fftwf_execute(plan); 64 fftwf_destroy_plan(plan); 65 66 float denom = length; 67 for (size_t k = 0; k < m_complex_vector_size; ++k) 68 { 69 m_gamma[k][0] /= denom; 70 m_gamma[k][1] /= denom; 71 } 72 73 if (length % 2 == 0) 74 { 75 m_gamma[m_complex_vector_size -1][0] /= 2; 76 // numerically, m_gamma[m_complex_vector_size -1][1] should be zero . . . 77 // I believe, but need to check, that FFTW guarantees that it is identically zero. 78 } 79 } 80 81 cardinal_trigonometric_detail(const cardinal_trigonometric_detail& old) = delete; 82 83 cardinal_trigonometric_detail& operator=(const cardinal_trigonometric_detail&) = delete; 84 85 cardinal_trigonometric_detail(cardinal_trigonometric_detail &&) = delete; 86 operator ()(float t) const87 float operator()(float t) const 88 { 89 using std::sin; 90 using std::cos; 91 using boost::math::constants::two_pi; 92 using std::exp; 93 float s = m_gamma[0][0]; 94 float x = two_pi<float>()*(t - m_t0)/m_T; 95 fftwf_complex z; 96 // boost::math::cos_pi with a redefinition of x? Not now . . . 97 z[0] = cos(x); 98 z[1] = sin(x); 99 fftwf_complex b{0, 0}; 100 // u = b*z 101 fftwf_complex u; 102 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) { 103 u[0] = b[0]*z[0] - b[1]*z[1]; 104 u[1] = b[0]*z[1] + b[1]*z[0]; 105 b[0] = m_gamma[k][0] + u[0]; 106 b[1] = m_gamma[k][1] + u[1]; 107 } 108 109 s += 2*(b[0]*z[0] - b[1]*z[1]); 110 return s; 111 } 112 prime(float t) const113 float prime(float t) const 114 { 115 using std::sin; 116 using std::cos; 117 using boost::math::constants::two_pi; 118 using std::exp; 119 float x = two_pi<float>()*(t - m_t0)/m_T; 120 fftwf_complex z; 121 z[0] = cos(x); 122 z[1] = sin(x); 123 fftwf_complex b{0, 0}; 124 // u = b*z 125 fftwf_complex u; 126 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) 127 { 128 u[0] = b[0]*z[0] - b[1]*z[1]; 129 u[1] = b[0]*z[1] + b[1]*z[0]; 130 b[0] = k*m_gamma[k][0] + u[0]; 131 b[1] = k*m_gamma[k][1] + u[1]; 132 } 133 // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1]) 134 return -2*two_pi<float>()*(b[1]*z[0] + b[0]*z[1])/m_T; 135 } 136 double_prime(float t) const137 float double_prime(float t) const 138 { 139 using std::sin; 140 using std::cos; 141 using boost::math::constants::two_pi; 142 using std::exp; 143 float x = two_pi<float>()*(t - m_t0)/m_T; 144 fftwf_complex z; 145 z[0] = cos(x); 146 z[1] = sin(x); 147 fftwf_complex b{0, 0}; 148 // u = b*z 149 fftwf_complex u; 150 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) 151 { 152 u[0] = b[0]*z[0] - b[1]*z[1]; 153 u[1] = b[0]*z[1] + b[1]*z[0]; 154 b[0] = k*k*m_gamma[k][0] + u[0]; 155 b[1] = k*k*m_gamma[k][1] + u[1]; 156 } 157 // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1]) 158 return -2*two_pi<float>()*two_pi<float>()*(b[0]*z[0] - b[1]*z[1])/(m_T*m_T); 159 } 160 period() const161 float period() const 162 { 163 return m_T; 164 } 165 integrate() const166 float integrate() const 167 { 168 return m_T*m_gamma[0][0]; 169 } 170 squared_l2() const171 float squared_l2() const 172 { 173 float s = 0; 174 // Always add smallest to largest for accuracy. 175 for (size_t i = m_complex_vector_size - 1; i >= 1; --i) 176 { 177 s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]); 178 } 179 s *= 2; 180 s += m_gamma[0][0]*m_gamma[0][0]; 181 return s*m_T; 182 } 183 184 ~cardinal_trigonometric_detail()185 ~cardinal_trigonometric_detail() 186 { 187 if (m_gamma) 188 { 189 fftwf_free(m_gamma); 190 m_gamma = nullptr; 191 } 192 } 193 194 195 private: 196 float m_t0; 197 float m_h; 198 float m_T; 199 fftwf_complex* m_gamma; 200 size_t m_complex_vector_size; 201 }; 202 203 204 template<> 205 class cardinal_trigonometric_detail<double> { 206 public: cardinal_trigonometric_detail(const double * data,size_t length,double t0,double h)207 cardinal_trigonometric_detail(const double* data, size_t length, double t0, double h) : m_t0{t0}, m_h{h} 208 { 209 if (length == 0) 210 { 211 throw std::logic_error("At least one sample is required."); 212 } 213 if (h <= 0) 214 { 215 throw std::logic_error("The step size must be > 0"); 216 } 217 m_T = m_h*length; 218 m_complex_vector_size = length/2 + 1; 219 m_gamma = fftw_alloc_complex(m_complex_vector_size); 220 fftw_plan plan = fftw_plan_dft_r2c_1d(length, const_cast<double*>(data), m_gamma, FFTW_ESTIMATE); 221 if (!plan) 222 { 223 throw std::logic_error("A null fftw plan was created."); 224 } 225 226 fftw_execute(plan); 227 fftw_destroy_plan(plan); 228 229 double denom = length; 230 for (size_t k = 0; k < m_complex_vector_size; ++k) 231 { 232 m_gamma[k][0] /= denom; 233 m_gamma[k][1] /= denom; 234 } 235 236 if (length % 2 == 0) 237 { 238 m_gamma[m_complex_vector_size -1][0] /= 2; 239 } 240 } 241 242 cardinal_trigonometric_detail(const cardinal_trigonometric_detail& old) = delete; 243 244 cardinal_trigonometric_detail& operator=(const cardinal_trigonometric_detail&) = delete; 245 246 cardinal_trigonometric_detail(cardinal_trigonometric_detail &&) = delete; 247 operator ()(double t) const248 double operator()(double t) const 249 { 250 using std::sin; 251 using std::cos; 252 using boost::math::constants::two_pi; 253 using std::exp; 254 double s = m_gamma[0][0]; 255 double x = two_pi<double>()*(t - m_t0)/m_T; 256 fftw_complex z; 257 z[0] = cos(x); 258 z[1] = sin(x); 259 fftw_complex b{0, 0}; 260 // u = b*z 261 fftw_complex u; 262 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) 263 { 264 u[0] = b[0]*z[0] - b[1]*z[1]; 265 u[1] = b[0]*z[1] + b[1]*z[0]; 266 b[0] = m_gamma[k][0] + u[0]; 267 b[1] = m_gamma[k][1] + u[1]; 268 } 269 270 s += 2*(b[0]*z[0] - b[1]*z[1]); 271 return s; 272 } 273 prime(double t) const274 double prime(double t) const 275 { 276 using std::sin; 277 using std::cos; 278 using boost::math::constants::two_pi; 279 using std::exp; 280 double x = two_pi<double>()*(t - m_t0)/m_T; 281 fftw_complex z; 282 z[0] = cos(x); 283 z[1] = sin(x); 284 fftw_complex b{0, 0}; 285 // u = b*z 286 fftw_complex u; 287 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) 288 { 289 u[0] = b[0]*z[0] - b[1]*z[1]; 290 u[1] = b[0]*z[1] + b[1]*z[0]; 291 b[0] = k*m_gamma[k][0] + u[0]; 292 b[1] = k*m_gamma[k][1] + u[1]; 293 } 294 // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1]) 295 return -2*two_pi<double>()*(b[1]*z[0] + b[0]*z[1])/m_T; 296 } 297 double_prime(double t) const298 double double_prime(double t) const 299 { 300 using std::sin; 301 using std::cos; 302 using boost::math::constants::two_pi; 303 using std::exp; 304 double x = two_pi<double>()*(t - m_t0)/m_T; 305 fftw_complex z; 306 z[0] = cos(x); 307 z[1] = sin(x); 308 fftw_complex b{0, 0}; 309 // u = b*z 310 fftw_complex u; 311 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) 312 { 313 u[0] = b[0]*z[0] - b[1]*z[1]; 314 u[1] = b[0]*z[1] + b[1]*z[0]; 315 b[0] = k*k*m_gamma[k][0] + u[0]; 316 b[1] = k*k*m_gamma[k][1] + u[1]; 317 } 318 // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1]) 319 return -2*two_pi<double>()*two_pi<double>()*(b[0]*z[0] - b[1]*z[1])/(m_T*m_T); 320 } 321 period() const322 double period() const 323 { 324 return m_T; 325 } 326 integrate() const327 double integrate() const 328 { 329 return m_T*m_gamma[0][0]; 330 } 331 squared_l2() const332 double squared_l2() const 333 { 334 double s = 0; 335 for (size_t i = m_complex_vector_size - 1; i >= 1; --i) 336 { 337 s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]); 338 } 339 s *= 2; 340 s += m_gamma[0][0]*m_gamma[0][0]; 341 return s*m_T; 342 } 343 ~cardinal_trigonometric_detail()344 ~cardinal_trigonometric_detail() 345 { 346 if (m_gamma) 347 { 348 fftw_free(m_gamma); 349 m_gamma = nullptr; 350 } 351 } 352 353 private: 354 double m_t0; 355 double m_h; 356 double m_T; 357 fftw_complex* m_gamma; 358 size_t m_complex_vector_size; 359 }; 360 361 362 template<> 363 class cardinal_trigonometric_detail<long double> { 364 public: cardinal_trigonometric_detail(const long double * data,size_t length,long double t0,long double h)365 cardinal_trigonometric_detail(const long double* data, size_t length, long double t0, long double h) : m_t0{t0}, m_h{h} 366 { 367 if (length == 0) 368 { 369 throw std::logic_error("At least one sample is required."); 370 } 371 if (h <= 0) 372 { 373 throw std::logic_error("The step size must be > 0"); 374 } 375 m_T = m_h*length; 376 m_complex_vector_size = length/2 + 1; 377 m_gamma = fftwl_alloc_complex(m_complex_vector_size); 378 fftwl_plan plan = fftwl_plan_dft_r2c_1d(length, const_cast<long double*>(data), m_gamma, FFTW_ESTIMATE); 379 if (!plan) 380 { 381 throw std::logic_error("A null fftw plan was created."); 382 } 383 384 fftwl_execute(plan); 385 fftwl_destroy_plan(plan); 386 387 long double denom = length; 388 for (size_t k = 0; k < m_complex_vector_size; ++k) 389 { 390 m_gamma[k][0] /= denom; 391 m_gamma[k][1] /= denom; 392 } 393 394 if (length % 2 == 0) { 395 m_gamma[m_complex_vector_size -1][0] /= 2; 396 } 397 } 398 399 cardinal_trigonometric_detail(const cardinal_trigonometric_detail& old) = delete; 400 401 cardinal_trigonometric_detail& operator=(const cardinal_trigonometric_detail&) = delete; 402 403 cardinal_trigonometric_detail(cardinal_trigonometric_detail &&) = delete; 404 operator ()(long double t) const405 long double operator()(long double t) const 406 { 407 using std::sin; 408 using std::cos; 409 using boost::math::constants::two_pi; 410 using std::exp; 411 long double s = m_gamma[0][0]; 412 long double x = two_pi<long double>()*(t - m_t0)/m_T; 413 fftwl_complex z; 414 z[0] = cos(x); 415 z[1] = sin(x); 416 fftwl_complex b{0, 0}; 417 fftwl_complex u; 418 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) 419 { 420 u[0] = b[0]*z[0] - b[1]*z[1]; 421 u[1] = b[0]*z[1] + b[1]*z[0]; 422 b[0] = m_gamma[k][0] + u[0]; 423 b[1] = m_gamma[k][1] + u[1]; 424 } 425 426 s += 2*(b[0]*z[0] - b[1]*z[1]); 427 return s; 428 } 429 prime(long double t) const430 long double prime(long double t) const 431 { 432 using std::sin; 433 using std::cos; 434 using boost::math::constants::two_pi; 435 using std::exp; 436 long double x = two_pi<long double>()*(t - m_t0)/m_T; 437 fftwl_complex z; 438 z[0] = cos(x); 439 z[1] = sin(x); 440 fftwl_complex b{0, 0}; 441 // u = b*z 442 fftwl_complex u; 443 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) 444 { 445 u[0] = b[0]*z[0] - b[1]*z[1]; 446 u[1] = b[0]*z[1] + b[1]*z[0]; 447 b[0] = k*m_gamma[k][0] + u[0]; 448 b[1] = k*m_gamma[k][1] + u[1]; 449 } 450 // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1]) 451 return -2*two_pi<long double>()*(b[1]*z[0] + b[0]*z[1])/m_T; 452 } 453 double_prime(long double t) const454 long double double_prime(long double t) const 455 { 456 using std::sin; 457 using std::cos; 458 using boost::math::constants::two_pi; 459 using std::exp; 460 long double x = two_pi<long double>()*(t - m_t0)/m_T; 461 fftwl_complex z; 462 z[0] = cos(x); 463 z[1] = sin(x); 464 fftwl_complex b{0, 0}; 465 // u = b*z 466 fftwl_complex u; 467 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) 468 { 469 u[0] = b[0]*z[0] - b[1]*z[1]; 470 u[1] = b[0]*z[1] + b[1]*z[0]; 471 b[0] = k*k*m_gamma[k][0] + u[0]; 472 b[1] = k*k*m_gamma[k][1] + u[1]; 473 } 474 // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1]) 475 return -2*two_pi<long double>()*two_pi<long double>()*(b[0]*z[0] - b[1]*z[1])/(m_T*m_T); 476 } 477 period() const478 long double period() const 479 { 480 return m_T; 481 } 482 integrate() const483 long double integrate() const 484 { 485 return m_T*m_gamma[0][0]; 486 } 487 squared_l2() const488 long double squared_l2() const 489 { 490 long double s = 0; 491 for (size_t i = m_complex_vector_size - 1; i >= 1; --i) 492 { 493 s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]); 494 } 495 s *= 2; 496 s += m_gamma[0][0]*m_gamma[0][0]; 497 return s*m_T; 498 } 499 ~cardinal_trigonometric_detail()500 ~cardinal_trigonometric_detail() 501 { 502 if (m_gamma) 503 { 504 fftwl_free(m_gamma); 505 m_gamma = nullptr; 506 } 507 } 508 509 private: 510 long double m_t0; 511 long double m_h; 512 long double m_T; 513 fftwl_complex* m_gamma; 514 size_t m_complex_vector_size; 515 }; 516 517 #ifdef BOOST_HAS_FLOAT128 518 template<> 519 class cardinal_trigonometric_detail<__float128> { 520 public: cardinal_trigonometric_detail(const __float128 * data,size_t length,__float128 t0,__float128 h)521 cardinal_trigonometric_detail(const __float128* data, size_t length, __float128 t0, __float128 h) : m_t0{t0}, m_h{h} 522 { 523 if (length == 0) 524 { 525 throw std::logic_error("At least one sample is required."); 526 } 527 if (h <= 0) 528 { 529 throw std::logic_error("The step size must be > 0"); 530 } 531 m_T = m_h*length; 532 m_complex_vector_size = length/2 + 1; 533 m_gamma = fftwq_alloc_complex(m_complex_vector_size); 534 fftwq_plan plan = fftwq_plan_dft_r2c_1d(length, reinterpret_cast<__float128*>(const_cast<__float128*>(data)), m_gamma, FFTW_ESTIMATE); 535 if (!plan) 536 { 537 throw std::logic_error("A null fftw plan was created."); 538 } 539 540 fftwq_execute(plan); 541 fftwq_destroy_plan(plan); 542 543 __float128 denom = length; 544 for (size_t k = 0; k < m_complex_vector_size; ++k) 545 { 546 m_gamma[k][0] /= denom; 547 m_gamma[k][1] /= denom; 548 } 549 if (length % 2 == 0) 550 { 551 m_gamma[m_complex_vector_size -1][0] /= 2; 552 } 553 } 554 555 cardinal_trigonometric_detail(const cardinal_trigonometric_detail& old) = delete; 556 557 cardinal_trigonometric_detail& operator=(const cardinal_trigonometric_detail&) = delete; 558 559 cardinal_trigonometric_detail(cardinal_trigonometric_detail &&) = delete; 560 operator ()(__float128 t) const561 __float128 operator()(__float128 t) const 562 { 563 using std::sin; 564 using std::cos; 565 using boost::math::constants::two_pi; 566 using std::exp; 567 __float128 s = m_gamma[0][0]; 568 __float128 x = two_pi<__float128>()*(t - m_t0)/m_T; 569 fftwq_complex z; 570 z[0] = cosq(x); 571 z[1] = sinq(x); 572 fftwq_complex b{0, 0}; 573 fftwq_complex u; 574 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) 575 { 576 u[0] = b[0]*z[0] - b[1]*z[1]; 577 u[1] = b[0]*z[1] + b[1]*z[0]; 578 b[0] = m_gamma[k][0] + u[0]; 579 b[1] = m_gamma[k][1] + u[1]; 580 } 581 582 s += 2*(b[0]*z[0] - b[1]*z[1]); 583 return s; 584 } 585 prime(__float128 t) const586 __float128 prime(__float128 t) const 587 { 588 using std::sin; 589 using std::cos; 590 using boost::math::constants::two_pi; 591 using std::exp; 592 __float128 x = two_pi<__float128>()*(t - m_t0)/m_T; 593 fftwq_complex z; 594 z[0] = cosq(x); 595 z[1] = sinq(x); 596 fftwq_complex b{0, 0}; 597 // u = b*z 598 fftwq_complex u; 599 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) 600 { 601 u[0] = b[0]*z[0] - b[1]*z[1]; 602 u[1] = b[0]*z[1] + b[1]*z[0]; 603 b[0] = k*m_gamma[k][0] + u[0]; 604 b[1] = k*m_gamma[k][1] + u[1]; 605 } 606 // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1]) 607 return -2*two_pi<__float128>()*(b[1]*z[0] + b[0]*z[1])/m_T; 608 } 609 double_prime(__float128 t) const610 __float128 double_prime(__float128 t) const 611 { 612 using std::sin; 613 using std::cos; 614 using boost::math::constants::two_pi; 615 using std::exp; 616 __float128 x = two_pi<__float128>()*(t - m_t0)/m_T; 617 fftwq_complex z; 618 z[0] = cosq(x); 619 z[1] = sinq(x); 620 fftwq_complex b{0, 0}; 621 // u = b*z 622 fftwq_complex u; 623 for (size_t k = m_complex_vector_size - 1; k >= 1; --k) 624 { 625 u[0] = b[0]*z[0] - b[1]*z[1]; 626 u[1] = b[0]*z[1] + b[1]*z[0]; 627 b[0] = k*k*m_gamma[k][0] + u[0]; 628 b[1] = k*k*m_gamma[k][1] + u[1]; 629 } 630 // b*z = (b[0]*z[0] - b[1]*z[1]) + i(b[1]*z[0] + b[0]*z[1]) 631 return -2*two_pi<__float128>()*two_pi<__float128>()*(b[0]*z[0] - b[1]*z[1])/(m_T*m_T); 632 } 633 period() const634 __float128 period() const 635 { 636 return m_T; 637 } 638 integrate() const639 __float128 integrate() const 640 { 641 return m_T*m_gamma[0][0]; 642 } 643 squared_l2() const644 __float128 squared_l2() const 645 { 646 __float128 s = 0; 647 for (size_t i = m_complex_vector_size - 1; i >= 1; --i) 648 { 649 s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]); 650 } 651 s *= 2; 652 s += m_gamma[0][0]*m_gamma[0][0]; 653 return s*m_T; 654 } 655 ~cardinal_trigonometric_detail()656 ~cardinal_trigonometric_detail() 657 { 658 if (m_gamma) 659 { 660 fftwq_free(m_gamma); 661 m_gamma = nullptr; 662 } 663 } 664 665 666 private: 667 __float128 m_t0; 668 __float128 m_h; 669 __float128 m_T; 670 fftwq_complex* m_gamma; 671 size_t m_complex_vector_size; 672 }; 673 #endif 674 675 }}}} 676 #endif 677