File indexing completed on 2023-10-25 09:39:31
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::reduceRange(phi);
0015 }
0016
0017 template <typename T>
0018 inline T deltaPhi(T phi1, T phi2) {
0019 return reco::reduceRange(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 }
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;
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;
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;
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;
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
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 }