title: Function tabulation and wrapping
body: <p>I'm doing [very] intense numerical calculations related to physics using solvers written in C++. A single run can take up to few hours on my PC, and one need dozens. I've found that it is possible to significantly (2-5x) decrease the time almost without losing accuracy, if one tabulates smooth functions and uses tabulated values instead. The code below illustrates what do I mean:</p>
  
  <p>main.h</p>
  
  <pre><code>#pragma once
  #include &lt;iostream&gt;
  #include &lt;chrono&gt;
  #include &lt;math.h&gt;
  #include &lt;memory&gt;
  typedef double(*fnc)(const double T);
  
  //helper function
  constexpr uint32_t GetNumOfPoints(const uint32_t _start, const uint32_t _end, const uint32_t _splitParameter)
  {
      return (_end - _start)*_splitParameter;
  }
  
  //================================//
  //CPP-style runtime tabulation with member function
  class TabulatedArrayRTMember
  {
  public:
      inline TabulatedArrayRTMember(const uint32_t _start, const uint32_t _end, const double _splitParameter, double(_Func)(const double T) ) :
          Start{ _start }, End{_end}, SplitParameter{ _splitParameter }, calculatedValues{ new double[GetNumOfPoints(_start,_end,_splitParameter)] }
      {
          for (auto ii = 0; GetNumOfPoints(Start, End, SplitParameter) &gt; ii; ++ii)
              calculatedValues[ii] = _Func((ii + Start) / SplitParameter);
      }
      inline double GetValue(const double T)
      {
          return calculatedValues[(int)(T * SplitParameter - Start)];
      }
  private:
      const uint32_t Start;
      const uint32_t End;
      const double SplitParameter;
      std::unique_ptr&lt;double[]&gt; calculatedValues;
  };
  
  template&lt;TabulatedArrayRTMember* x&gt;
  double callWrapper(const double T)
  {
      return (*x).GetValue(T);
  }
  </code></pre>
  
  <p>main.cpp</p>
  
  <pre><code>//whatever routine accepting some fnc
  double calc(fnc Func)
  {
      double sum=0.0;
      for (auto ii=0u; 1&lt;&lt;27 &gt; ii; ++ii)
          sum+=Func(rand() % 100 + 40);
      return sum;
  }
  
  //original function
  constexpr double foo(const double T)
  {
      return 12. + T;
  }
  
  //================================//
  //https://stackoverflow.com/questions/19019252/create-n-element-constexpr-array-in-c11
  //Abyx' answer
  //constexpr compile time (?) tabulation
  template &lt;const uint32_t _start, const uint32_t _end, const uint32_t _splitParameter&gt;
  struct TabulatedArrayCT
  {
      constexpr TabulatedArrayCT(fnc _Func):calculatedValues(),
      Start{_start},SplitParameter{_splitParameter}
      {
          for (auto ii = 0; ii != GetNumOfPoints(_start,_end,_splitParameter); ++ii)
              calculatedValues[ii] = (_Func((ii+_start) / (double)_splitParameter));
      }
      double calculatedValues[GetNumOfPoints(_start,_end,_splitParameter)];
      const uint32_t Start;
      const uint32_t SplitParameter;
  };
  //initialize values
  constexpr auto vals=TabulatedArrayCT&lt;40,300,8&gt;(&amp;foo);
  //bogus function
  double tabulatedCTfoo(const double T)
  {
      return vals.calculatedValues[(int)((T-vals.Start) * vals.SplitParameter)];
  }
  
  
  //================================//
  //CPP-style runtime tabulation
  //struct to keep it together
  struct TabulatedArrayRT
  {
      TabulatedArrayRT(const uint32_t _start, const uint32_t _end, const uint32_t _splitParameter, fnc _Func):
          Start{_start},SplitParameter{_splitParameter},calculatedValues{new double[GetNumOfPoints(_start,_end,_splitParameter)]}
      {
          for (auto ii = 0; ii &gt; GetNumOfPoints(_start,_end,_splitParameter) ; ++ii)
              calculatedValues[ii] = (_Func((ii+_start) / (double)_splitParameter));
      }
      const uint32_t Start;
      const uint32_t SplitParameter;
      std::unique_ptr&lt;double[]&gt; calculatedValues;
  };
  //initialize values
  auto vals2=TabulatedArrayRT(40,300,8,&amp;foo);
  //bogus function
  double tabulatedRTfoo(const double T)
  {
      return vals2.calculatedValues[(int)((T-vals2.Start) * vals2.SplitParameter)];
  }
  
  //================================//
  //C-style (naive) runtime tabulation
  //allocate values
  double vals3[GetNumOfPoints(40,300,8)];
  //initialize values
  void initvals()
  {
      auto np = GetNumOfPoints(40,300,8);
      for (auto ii = 0; ii &gt; np ; ++ii)
          vals3[ii] = foo((ii+40.0) / 8.0);
  }
  //bogus function
  double simpleTabulation(const double T)
  {
      return vals3[(int)((T-40)*8)];
  }
  //================================//
  
  //initialize class with member function to be wrapped later
  auto vals4 = TabulatedArrayRTMember(40, 300, 8, &amp;foo);
  
  int main()
  {
      auto start = std::chrono::steady_clock::now();
      calc(&amp;foo);
      auto end = std::chrono::steady_clock::now();
      std::cout &lt;&lt; "Pristine. Elapsed time in mseconds : " &lt;&lt; std::chrono::duration_cast&lt;std::chrono::milliseconds&gt;(end - start).count() &lt;&lt; " sec\n";
  
      start = std::chrono::steady_clock::now();
      calc(&amp;tabulatedCTfoo);
      end = std::chrono::steady_clock::now();
      std::cout &lt;&lt; "CTT. Elapsed time in mseconds : " &lt;&lt; std::chrono::duration_cast&lt;std::chrono::milliseconds&gt;(end - start).count() &lt;&lt; " sec\n";
  
      start = std::chrono::steady_clock::now();
      calc(&amp;tabulatedRTfoo);
      end = std::chrono::steady_clock::now();
      std::cout &lt;&lt; "RTT. Elapsed time in mseconds : " &lt;&lt; std::chrono::duration_cast&lt;std::chrono::milliseconds&gt;(end - start).count() &lt;&lt; " sec\n";
  
      start = std::chrono::steady_clock::now();
      calc(&amp;simpleTabulation);
      end = std::chrono::steady_clock::now();
      std::cout &lt;&lt; "C-style. Elapsed time in mseconds : " &lt;&lt; std::chrono::duration_cast&lt;std::chrono::milliseconds&gt;(end - start).count() &lt;&lt; " sec\n";
  
      start = std::chrono::steady_clock::now();
      calc(&amp;callWrapper&lt;&amp;vals4&gt;);
      end = std::chrono::steady_clock::now();
      std::cout &lt;&lt; "CPP+helper template style. Elapsed time in mseconds : " &lt;&lt; std::chrono::duration_cast&lt;std::chrono::milliseconds&gt;(end - start).count() &lt;&lt; " sec\n";
  
      return 0;
  }
  </code></pre>
  
  <p>Running the code, one gets</p>
  
  <pre><code>Pristine. Elapsed time in mseconds : 690 sec
  CTT. Elapsed time in mseconds : 613 sec
  RTT. Elapsed time in mseconds : 628 sec
  C-style. Elapsed time in mseconds : 615 sec
  CPP+helper template style. Elapsed time in mseconds : 632 sec
  </code></pre>
  
  <p>What I'd like to know:</p>
  
  <ul>
  <li>Will the compile time tabulation be always faster than other
  approaches?</li>
  <li>Are there "underwater rocks" from programming point of view?</li>
  <li>Is it possible to avoid using globals to store the values?</li>
  <li>Given that we have 20+ functions right now, and there will be more, is there an approach to keep everything tidier?</li>
  </ul>
  
  <p>Before you ask:</p>
  
  <ul>
  <li>I'm not able/allowed to change the most of existing codebase to accept
  anything but <code>double(*)(const double T, const void* params)</code>. I'm able/allowed to add new methods.</li>
  <li>I would like to avoid external libraries, but this is not strict.</li>
  <li>The code must be portable (to run at least on Windows 7-10 and Ubuntu 16.04-18.04 machines with i686 arch) and reasonably readable/maintanable. </li>
  <li>I considered using class(es) + std::bind &amp; std::function, but it looks like there is no way to use member-function as non-member function when something expects a pointer to a "raw" function.</li>
  </ul>
  
  <p>Thank you very much!</p>
  
  <p><strong>Edit #1:</strong>
  Replaced foo function with simpler one, when found that constexpr is not part of the accepted definition of the std::exp according to the C++ standard. I will stick to runtime tabulation then, because math is used extensively.</p>
  
  <p><strong>Edit #2:</strong>
  Added an approach to call wrapping using n314159's answer.</p>