Cardinal trigonometric interpolation: Implement squared l2 norm.

doc/interpolators/cardinal_trigonometric.qbk

@@ -25,6 +25,8 @@ public:
     Real period() const;
 
     Real integrate() const;
+
+    Real squared_l2() const;
 };
 }}}
 ```
@@ -54,7 +56,7 @@ The period is always given by `v.size()*h`.
 Off-by-one errors are common in programming, and hence if you wonder what this interpolator believes the period to be, you can query it with the `.period()` member function.
 
 In addition, the transform into the trigonometric basis gives a trivial way to compute the integral of the function over a period; this is done via the `.integrate()` member function.
-
+Evaluation of the square of the L[super 2] norm is trivial in this basis; it is computed by the `.squared_l2()` member function.
 
 [heading Caveats]
 

include/boost/math/interpolators/detail/cardinal_trigonometric_detail.hpp

@@ -86,6 +86,7 @@ public:
       float s = m_gamma[0][0];
       float x = two_pi<float>()*(t - m_t0)/m_T;
       fftwf_complex z;
+      // boost::math::cos_pi with a redefinition of x? Not now . . .
       z[0] = cos(x);
       z[1] = sin(x);
       fftwf_complex b{0, 0};
@@ -107,10 +108,24 @@ public:
     return m_T;
   }
 
-  float integrate() const {
+  float integrate() const
+  {
       return m_T*m_gamma[0][0];
   }
 
+  float squared_l2() const
+  {
+    float s = 0;
+    // Always add smallest to largest for accuracy.
+    for (size_t i = m_complex_vector_size - 1; i >= 1; --i)
+    {
+        s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]);
+    }
+    s *= 2;
+    s += m_gamma[0][0]*m_gamma[0][0];
+    return s*m_T;
+  }
+
   ~cardinal_trigonometric_detail()
   {
     if (m_gamma)
@@ -205,6 +220,18 @@ public:
     return m_T*m_gamma[0][0];
   }
 
+  double squared_l2() const
+  {
+    double s = 0;
+    for (size_t i = m_complex_vector_size - 1; i >= 1; --i)
+    {
+        s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]);
+    }
+    s *= 2;
+    s += m_gamma[0][0]*m_gamma[0][0];
+    return s*m_T;
+  }
+
   ~cardinal_trigonometric_detail()
   {
     if (m_gamma)
@@ -292,11 +319,23 @@ public:
     return m_T;
   }
 
-  double integrate() const
+  long double integrate() const
   {
     return m_T*m_gamma[0][0];
   }
 
+  long double squared_l2() const
+  {
+    long double s = 0;
+    for (size_t i = m_complex_vector_size - 1; i >= 1; --i)
+    {
+        s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]);
+    }
+    s *= 2;
+    s += m_gamma[0][0]*m_gamma[0][0];
+    return s*m_T;
+  }
+
   ~cardinal_trigonometric_detail()
   {
     if (m_gamma)
@@ -389,6 +428,17 @@ public:
     return m_T*m_gamma[0][0];
   }
 
+  __float128 squared_l2() const
+  {
+    __float128 s = 0;
+    for (size_t i = m_complex_vector_size - 1; i >= 1; --i)
+    {
+        s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]);
+    }
+    s *= 2;
+    s += m_gamma[0][0]*m_gamma[0][0];
+    return s*m_T;
+  }
 
   ~cardinal_trigonometric_detail()
   {

test/cardinal_trigonometric_test.cpp

@@ -31,6 +31,8 @@ void test_constant()
       auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
       CHECK_ULP_CLOSE(c, ct(0.3), 3);
       CHECK_ULP_CLOSE(c*h*n, ct.integrate(), 3);
+
+      CHECK_ULP_CLOSE(c*c*h*n, ct.squared_l2(), 3);
     }
 }
 
@@ -119,7 +121,12 @@ void test_bump()
 
   // Wolfram Alpha:
   // NIntegrate[Exp[-1/(1-x*x)],{x,-1,1}]
-  CHECK_ULP_CLOSE(Real(0.4439938161680794378), ct.integrate(), 3);
+  CHECK_ULP_CLOSE(Real(0.443993816168079437823L), ct.integrate(), 3);
+
+  // NIntegrate[Exp[-2/(1-x*x)],{x,-1,1}]
+  CHECK_ULP_CLOSE(Real(0.1330861208449942715569473279553285713625791551628130055345002588895389L), ct.squared_l2(), 1);
+
+
 }
 
 
@@ -132,6 +139,7 @@ int main()
     test_bump<float>();
 #endif
 
+
 #ifdef TEST2
     test_constant<double>();
     test_sampled_sine<double>();