Math/test/testAtan2.cpp

0001 #include "DataFormats/Math/interface/approx_atan2.h"
0002 #include "DataFormats/Math/interface/deltaPhi.h"
0003
0004 #include <cstdio>
0005 #include <cstdlib>
0006 #include <iostream>
0007 #include <cassert>
0008
0009 #include "rdtscp.h"
0010
0011 namespace {
0012   template <typename T>
0013   inline T toPhi(T phi) {
0014     return reco::reducePhiRange(phi);
0015   }
0016
0017   template <typename T>
0018   inline T deltaPhi(T phi1, T phi2) {
0019     return reco::reducePhiRange(phi1 - phi2);
0020   }
0021
0022   inline bool phiLess(float x, float y) {
0023     auto ix = phi2int(toPhi(x));
0024     auto iy = phi2int(toPhi(y));
0025
0026     return (ix - iy) < 0;
0027   }
0028
0029 }  // namespace
0030
0031 template <>
0032 float unsafe_atan2f<0>(float, float) {
0033   return 1.f;
0034 }
0035 template <>
0036 float unsafe_atan2f<99>(float y, float x) {
0037   return std::atan2(y, x);
0038 }
0039 template <>
0040 int unsafe_atan2i<0>(float, float) {
0041   return 1;
0042 }
0043 template <>
0044 int unsafe_atan2i<99>(float y, float x) {
0045   return phi2int(std::atan2(y, x));
0046 }
0047
0048 template <int DEGREE>
0049 void diff() {
0050   float mdiff = 0;
0051   int idiff = 0;
0052   constexpr float xmin = -100.001;  // avoid 0
0053   constexpr float incr = 0.04;
0054   constexpr int N = 2. * std::abs(xmin) / incr;
0055   auto y = xmin;
0056   for (int i = 0; i < N; ++i) {
0057     auto x = xmin;
0058     for (int i = 0; i < N; ++i) {
0059       auto approx = unsafe_atan2f<DEGREE>(y, x);
0060       auto iapprox = unsafe_atan2i<DEGREE>(y, x);
0061       auto std = std::atan2(y, x);
0062       mdiff = std::max(std::abs(std - approx), mdiff);
0063       idiff = std::max(std::abs(phi2int(std) - iapprox), idiff);
0064       x += incr;
0065     }
0066     y += incr;
0067   }
0068
0069   std::cout << "for degree " << DEGREE << " max diff is " << mdiff << ' ' << idiff << ' ' << int2phi(idiff)
0070             << std::endl;
0071 }
0072
0073 template <int DEGREE>
0074 void speedf() {
0075   if (!has_rdtscp())
0076     return;
0077   long long t = -rdtsc();
0078   float approx = 0;
0079   constexpr float xmin = -400.001;  // avoid 0
0080   constexpr float incr = 0.02;
0081   constexpr int N = 2. * std::abs(xmin) / incr;
0082   auto y = xmin;
0083   for (int i = 0; i < N; ++i) {
0084     auto x = xmin;
0085     for (int i = 0; i < N; ++i) {
0086       approx += unsafe_atan2f<DEGREE>(y, x);
0087       x += incr;
0088     }
0089     y += incr;
0090   }
0091   t += rdtsc();
0092
0093   std::cout << "f for degree " << DEGREE << " clock is " << t << " " << approx << std::endl;
0094 }
0095
0096 template <int DEGREE>
0097 void speeds() {
0098   if (!has_rdtscp())
0099     return;
0100   long long t = -rdtsc();
0101   float approx = 0;
0102   constexpr float xmin = -400.001;  // avoid 0
0103   constexpr float incr = 0.02;
0104   constexpr int N = 2. * std::abs(xmin) / incr;
0105   auto y = xmin;
0106   for (int i = 0; i < N; ++i) {
0107     auto x = xmin;
0108     for (int i = 0; i < N; ++i) {
0109       approx += safe_atan2f<DEGREE>(y, x);
0110       x += incr;
0111     }
0112     y += incr;
0113   }
0114   t += rdtsc();
0115
0116   std::cout << "s for degree " << DEGREE << " clock is " << t << " " << approx << std::endl;
0117 }
0118
0119 template <int DEGREE>
0120 void speedi() {
0121   if (!has_rdtscp())
0122     return;
0123   long long t = -rdtsc();
0124   int approx = 0;
0125   constexpr float xmin = -400.001;  // avoid 0
0126   constexpr float incr = 0.02;
0127   constexpr int N = 2. * std::abs(xmin) / incr;
0128   auto y = xmin;
0129   for (int i = 0; i < N; ++i) {
0130     auto x = xmin;
0131     for (int i = 0; i < N; ++i) {
0132       approx += unsafe_atan2i<DEGREE>(y, x);
0133       x += incr;
0134     }
0135     y += incr;
0136   }
0137   t += rdtsc();
0138
0139   std::cout << "i for degree " << DEGREE << " clock is " << t << " " << approx << std::endl;
0140 }
0141
0142 void testIntPhi() {
0143   constexpr long long maxint = (long long)(std::numeric_limits<int>::max()) + 1LL;
0144   constexpr int pi2 = int(maxint / 2LL);
0145   constexpr int pi4 = int(maxint / 4LL);
0146   constexpr int pi34 = int(3LL * maxint / 4LL);
0147
0148   std::cout << "pi,  pi2,  pi4, p34 " << maxint << ' ' << pi2 << ' ' << pi4 << ' ' << pi34 << ' ' << pi2 + pi4 << '\n';
0149   std::cout << "Maximum value for int: " << std::numeric_limits<int>::max() << '\n';
0150   std::cout << "Maximum value for int+1 as LL: " << (long long)(std::numeric_limits<int>::max()) + 1LL << std::endl;
0151
0152   std::cout << "Maximum value for short: " << std::numeric_limits<short>::max() << '\n';
0153   std::cout << "Maximum value for short+1 as int: " << (int)(std::numeric_limits<short>::max()) + 1 << std::endl;
0154
0155   auto d = float(M_PI) - std::nextafter(float(M_PI), 0.f);
0156   std::cout << "abs res at pi for float " << d << ' ' << phi2int(d) << std::endl;
0157   std::cout << "abs res at for int " << int2dphi(1) << std::endl;
0158   std::cout << "abs res at for short " << short2phi(1) << std::endl;
0159
0160   assert(phiLess(0.f, 2.f));
0161   assert(phiLess(6.f, 0.f));
0162   assert(phiLess(3.2f, 0.f));
0163   assert(phiLess(3.0f, 3.2f));
0164
0165   assert(phiLess(-0.3f, 0.f));
0166   assert(phiLess(-0.3f, 0.1f));
0167   assert(phiLess(-3.0f, 0.f));
0168   assert(phiLess(3.0f, -3.0f));
0169   assert(phiLess(0.f, -3.4f));
0170
0171   // go around the clock
0172   constexpr float eps = 1.e-5;
0173   auto ok = [](float a, float b) { assert(std::abs(a - b) < eps); };
0174   float phi1 = -7.;
0175   while (phi1 < 8) {
0176     auto p1 = toPhi(phi1);
0177     auto ip1 = phi2int(p1);
0178     std::cout << "phi1 " << phi1 << ' ' << p1 << ' ' << ip1 << ' ' << int2phi(ip1) << std::endl;
0179     ok(p1, int2phi(ip1));
0180     float phi2 = -7.2;
0181     while (phi2 < 8) {
0182       auto p2 = toPhi(phi2);
0183       auto ip2 = phi2int(p2);
0184       std::cout << "phi2 " << phi2 << ' ' << deltaPhi(phi1, phi2) << ' ' << deltaPhi(phi2, phi1) << ' '
0185                 << int2phi(ip1 - ip2) << ' ' << int2phi(ip2 - ip1) << ' ' << toPhi(phi2 + phi1) << ' '
0186                 << int2phi(ip1 + ip2) << std::endl;
0187       ok(deltaPhi(phi1, phi2), int2phi(ip1 - ip2));
0188       ok(deltaPhi(phi2, phi1), int2phi(ip2 - ip1));
0189       ok(toPhi(phi2 + phi1), int2phi(ip1 + ip2));
0190       phi2 += 1;
0191     }
0192
0193     phi1 += 1;
0194   }
0195 }
0196
0197 int main() {
0198   constexpr float p2i = ((int)(std::numeric_limits<short>::max()) + 1) / M_PI;
0199   constexpr float i2p = M_PI / ((int)(std::numeric_limits<short>::max()) + 1);
0200   std::cout << std::hexfloat << "p2i i2p " << p2i << " " << i2p << std::defaultfloat << std::endl;
0201
0202   std::cout << unsafe_atan2f<5>(0.f, 0.f) << " " << std::atan2(0., 0.) << std::endl;
0203   std::cout << unsafe_atan2f<5>(0.5f, 0.5f) << " " << std::atan2(0.5, 0.5) << std::endl;
0204   std::cout << unsafe_atan2f<5>(0.5f, -0.5f) << " " << std::atan2(0.5, -0.5) << std::endl;
0205   std::cout << unsafe_atan2f<5>(-0.5f, -0.5f) << " " << std::atan2(-0.5, -0.5) << std::endl;
0206   std::cout << unsafe_atan2f<5>(-0.5f, 0.5f) << " " << std::atan2(-0.5, 0.5) << std::endl;
0207
0208   std::cout << safe_atan2f<15>(0.f, 0.f) << " " << std::atan2(0., 0.) << std::endl;
0209   std::cout << safe_atan2f<15>(0.5f, 0.5f) << " " << std::atan2(0.5, 0.5) << std::endl;
0210   std::cout << safe_atan2f<15>(0.5f, -0.5f) << " " << std::atan2(0.5, -0.5) << std::endl;
0211   std::cout << safe_atan2f<15>(-0.5f, -0.5f) << " " << std::atan2(-0.5, -0.5) << std::endl;
0212   std::cout << safe_atan2f<15>(-0.5f, 0.5f) << " " << std::atan2(-0.5, 0.5) << std::endl;
0213
0214   testIntPhi();
0215
0216   diff<3>();
0217   diff<5>();
0218   diff<7>();
0219   diff<9>();
0220   diff<11>();
0221   diff<13>();
0222   diff<15>();
0223   diff<99>();
0224
0225   speedf<0>();
0226   speedf<3>();
0227   speedf<5>();
0228   speedf<7>();
0229   speedf<9>();
0230   speedf<11>();
0231   speedf<13>();
0232   speedf<15>();
0233   speedf<99>();
0234
0235   speeds<5>();
0236   speeds<11>();
0237   speeds<15>();
0238
0239   speedi<0>();
0240   speedi<3>();
0241   speedi<5>();
0242   speedi<7>();
0243   speedi<9>();
0244   speedi<11>();
0245   speedi<13>();
0246   speedi<99>();
0247
0248   return 0;
0249 }