boost/math/interpolators/detail/cardinal_trigonometric_detail.hpp

   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:
  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:
  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
  87   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
 113   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
 137   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
 161   float period() const
 162   {
 163     return m_T;
 164   }
 165
 166   float integrate() const
 167   {
 168     return m_T*m_gamma[0][0];
 169   }
 170
 171   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
 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:
 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
 248   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
 274   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
 298   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
 322   double period() const
 323   {
 324     return m_T;
 325   }
 326
 327   double integrate() const
 328   {
 329     return m_T*m_gamma[0][0];
 330   }
 331
 332   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
 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:
 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
 405   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
 430   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
 454   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
 478   long double period() const
 479   {
 480     return m_T;
 481   }
 482
 483   long double integrate() const
 484   {
 485     return m_T*m_gamma[0][0];
 486   }
 487
 488   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
 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:
 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
 561   __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
 586   __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
 610   __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
 634   __float128 period() const
 635   {
 636     return m_T;
 637   }
 638
 639   __float128 integrate() const
 640   {
 641     return m_T*m_gamma[0][0];
 642   }
 643
 644   __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
 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