Math/test/FastMath_t.cpp

0001 #include "DataFormats/Math/interface/FastMath.h"
0002
0003 #include <iostream>
0004 // #include <pmmintrin.h>
0005 #include <typeinfo>
0006
0007 #include "FWCore/Utilities/interface/HRRealTime.h"
0008
0009 namespace {
0010
0011   template <typename T>
0012   std::pair<T, T> stdatan2r(T x, T y) {
0013     return std::pair<T, T>(std::atan2(x, y), std::sqrt(x * x + y * y));
0014   }
0015
0016   template <typename T>
0017   struct Stat {
0018     std::string name;
0019     size_t n;
0020     size_t npos;
0021     T bias;
0022     T ave;
0023     T rms;
0024     T amax;
0025     Stat(std::string in) : name(in), n(0), npos(0), bias(0), ave(0), rms(0), amax(0) {}
0026     void operator()(T x, T ref);
0027
0028     ~Stat() {
0029       std::cout << name << " " << n << " " << npos << " " << bias / n << " " << ave / n << " "
0030                 << (n * rms - bias * bias) / (n * (n - 1)) << " " << amax << std::endl;
0031     }
0032   };
0033
0034   template <typename T>
0035   void Stat<T>::operator()(T x, T ref) {
0036     n++;
0037     if (x > ref)
0038       npos++;
0039     T d = (x - ref) / std::abs(ref);
0040     bias += d;
0041     ave += std::abs(d);
0042     rms += d * d;
0043     amax = std::max(amax, std::abs(d));
0044   }
0045
0046   double dummy;
0047   template <typename T>
0048   inline T eta(T x, T y, T z) {
0049     x = z / std::sqrt(x * x + y * y);
0050     return std::log(x + std::sqrt(x * x + T(1)));
0051   }
0052
0053   template <typename T>
0054   void sampleSquare() {
0055     edm::HRTimeType tf = 0;
0056     edm::HRTimeType ts = 0;
0057     Stat<T> stata("atan2");
0058     Stat<T> statr("r");
0059     T fac[8] = {-8, -5., -2., -1., 1., 2., 5., 8.};
0060     for (int k = 0; k < 100; k++)
0061       for (T x = 1e-15; x < 1.1e+15; x *= 10)
0062         for (T y = 1e-15; y < 1.1e+15; y *= 10)
0063           for (int i = 0; i != 8; ++i)
0064             for (int j = 0; j != 8; ++j) {
0065               T xx = x * fac[i];
0066               T yy = y * fac[j];
0067               edm::HRTimeType sf = edm::hrRealTime();
0068               std::pair<T, T> res = fastmath::atan2r(xx, yy);
0069               tf += (edm::hrRealTime() - sf);
0070               for (int l = 0; l < i + j; ++l)
0071                 dummy += yy;  // add a bit of random instruction
0072               edm::HRTimeType ss = edm::hrRealTime();
0073               std::pair<T, T> ref = stdatan2r(xx, yy);
0074               ts += (edm::hrRealTime() - ss);
0075               stata(res.first, ref.first);
0076               statr(res.second, ref.second);
0077             }
0078     std::cout << typeid(T).name() << " times " << tf << " " << ts << std::endl;
0079   }
0080
0081 }  // namespace
0082
0083 int main() {
0084   // _mm_setcsr (_mm_getcsr () | 0x8040);    // on Intel, treat denormals as zero for full speed
0085
0086   {
0087     std::pair<double, double> res = fastmath::atan2r(-3., 4.);
0088     std::cout << res.first << " " << res.second << std::endl;
0089     std::cout << atan2(-3., 4.) << std::endl;
0090   }
0091   {
0092     std::pair<double, double> res = fastmath::etaphi(-3., 4., 5.);
0093     std::cout << res.first << " " << res.second << std::endl;
0094     std::cout << eta(-3., 4., 5.) << " " << std::atan2(4., -3.) << std::endl;
0095   }
0096
0097   {
0098     std::pair<float, float> res = fastmath::atan2r(3.f, -4.f);
0099     std::cout << res.first << " " << res.second << std::endl;
0100     std::cout << atan2f(3.f, -4.f) << std::endl;
0101   }
0102   {
0103     std::pair<double, double> res = fastmath::etaphi(3.f, -4.f, -5.f);
0104     std::cout << res.first << " " << res.second << std::endl;
0105     std::cout << eta(3.f, -4.f, -5.f) << " " << std::atan2(-4.f, 3.f) << std::endl;
0106   }
0107
0108   sampleSquare<float>();
0109   sampleSquare<double>();
0110   return 0;
0111 }