Math/test/testApproxMath.cpp

0001 #include "DataFormats/Math/interface/approx_exp.h"
0002 #include "DataFormats/Math/interface/approx_log.h"
0003 #include "DataFormats/Math/interface/approx_erf.h"
0004 #include <cstdio>
0005 #include <cstdlib>
0006 #include <iostream>
0007
0008 void printEXP(float x) {
0009   printf("x= %a val= %a %a %a\n", x, std::exp(double(x)), unsafe_expf<6>(x), unsafe_expf<4>(x));
0010   printf("x= %f val= %f %f %f\n", x, std::exp(double(x)), unsafe_expf<6>(x), unsafe_expf<4>(x));
0011   printf("x= %f val= %f %f %f\n\n", x, std::exp(x), approx_expf<6>(x), approx_expf<4>(x));
0012 }
0013
0014 void printLOG(float x) {
0015   printf("x= %a val= %a %a %a\n", x, std::log(double(x)), unsafe_logf<8>(x), unsafe_logf<4>(x));
0016   printf("x= %f val= %f %f %f\n", x, std::log(double(x)), unsafe_logf<8>(x), unsafe_logf<4>(x));
0017   printf("x= %f val= %f %f %f\n\n", x, std::log(x), approx_logf<6>(x), approx_logf<4>(x));
0018 }
0019
0020 template <typename PRINT>
0021 void printLim(PRINT printit) {
0022   constexpr float zero_threshold_ftz = -float(0x5.75628p4);
0023   constexpr float inf_threshold = float(0x5.8b90cp4);
0024   std::cout << "\nexp valid for " << zero_threshold_ftz << " < x < " << inf_threshold << std::endl;
0025   printit(zero_threshold_ftz);
0026   printit(zero_threshold_ftz - 1);
0027   printit(zero_threshold_ftz + 1);
0028   printit(4 * zero_threshold_ftz);
0029   printit(0);
0030   printit(1.);
0031   printit(-1.);
0032   printit(inf_threshold);
0033   printit(inf_threshold + 1);
0034   printit(inf_threshold - 1);
0035   printit(4 * inf_threshold);
0036   printit(std::sqrt(-1));
0037   printit(std::numeric_limits<float>::infinity());
0038   printit(-std::numeric_limits<float>::infinity());
0039
0040   std::cout << "\n" << std::endl;
0041 }
0042
0043 namespace justcomp {
0044   constexpr int NN = 1024 * 1024;
0045   float a[NN], b[NN];
0046   template <int DEGREE>
0047   void bar() {
0048     for (int i = 0; i != NN; ++i)
0049       b[i] = approx_expf<DEGREE>(a[i]);
0050   }
0051   template <int DEGREE>
0052   void foo() {
0053     for (int i = 0; i != NN; ++i)
0054       b[i] = unsafe_expf<DEGREE>(a[i]);
0055   }
0056   template <int DEGREE>
0057   void lar() {
0058     for (int i = 0; i != NN; ++i)
0059       b[i] = approx_logf<DEGREE>(a[i]);
0060   }
0061   template <int DEGREE>
0062   void loo() {
0063     for (int i = 0; i != NN; ++i)
0064       b[i] = unsafe_logf<DEGREE>(a[i]);
0065   }
0066
0067 }  // namespace justcomp
0068
0069 template <typename STD, typename APPROX>
0070 void accTest(STD stdf, APPROX approx, int degree) {
0071   using namespace approx_math;
0072   std::cout << std::endl << "launching  exhaustive test for degree " << degree << std::endl;
0073   binary32 x, r, ref;
0074   int maxdiff = 0;
0075   int n127 = 0;
0076   int n16393 = 0;
0077   x.ui32 = 0;  // should be 0 but
0078   while (x.ui32 < 0xffffffff) {
0079     x.ui32++;
0080     // remove nans..
0081     if ((x.ui32 & 0x7f80000) && x.ui32 & 0x7FFFFF)
0082       continue;
0083     r.f = approx(x.f);
0084     ref.f = stdf(double(x.f));  // double-prec one  (no hope with -fno-math-errno)
0085     int d = abs(r.i32 - ref.i32);
0086     if (d > maxdiff) {
0087       // std::cout << "new maxdiff for x=" << x.f << " : " << d << std::endl;
0088       maxdiff = d;
0089     }
0090     if (d > 127)
0091       ++n127;
0092     if (d > 16393)
0093       ++n16393;
0094   }
0095   std::cout << "maxdiff / diff >127 / diff >16393 " << maxdiff << " / " << n127 << " / " << n16393 << std::endl;
0096 }
0097
0098 template <typename STD, typename APPROX>
0099 void accuTest(STD stdf,
0100               APPROX approx,
0101               const char* name,
0102               float mm = std::numeric_limits<float>::min(),
0103               float mx = std::numeric_limits<float>::max()) {
0104   using namespace approx_math;
0105   std::cout << std::endl << "launching  exhaustive test for " << name << std::endl;
0106   binary32 x, pend, r, ref;
0107   int maxdiff = 0;
0108   int n127 = 0;
0109   int n16393 = 0;
0110   float ad = 0., rd = 0;
0111   x.f = mm;
0112   x.ui32++;
0113   pend.f = mx;
0114   pend.ui32--;
0115   std::cout << "limits " << x.f << " " << pend.f << " " << pend.ui32 - x.ui32 << std::endl;
0116   while (x.ui32 < pend.ui32) {
0117     x.ui32++;
0118     r.f = approx(x.f);
0119     ref.f = stdf(x.f);  // double-prec one  (no hope with -fno-math-errno)
0120     ad = std::max(ad, std::abs(r.f - ref.f));
0121     rd = std::max(rd, std::abs((r.f / ref.f) - 1.f));
0122     int d = abs(r.i32 - ref.i32);
0123     if (d > maxdiff) {
0124       // std::cout << "new maxdiff for x=" << x.f << " : " << d << std::endl;
0125       maxdiff = d;
0126     }
0127     if (d > 127)
0128       ++n127;
0129     if (d > 16393)
0130       ++n16393;
0131   }
0132   std::cout << "absdiff / reldeff/ maxdiff / diff >127 / diff >16393 :  " << ad << " / " << rd << " / " << maxdiff
0133             << " / " << n127 << " / " << n16393 << std::endl;
0134 }
0135
0136 // performance test
0137 #include "rdtscp.h"
0138
0139 template <int DEGREE, int WHAT>
0140 struct Measure {
0141   inline void operator()(unsigned long long& t) const;
0142 };
0143
0144 template <int DEGREE>
0145 struct Measure<DEGREE, 0> {
0146   inline void operator()(unsigned long long& t) const {
0147     t -= rdtsc();
0148     justcomp::foo<DEGREE>();
0149     t += rdtsc();
0150   }
0151 };
0152
0153 template <int DEGREE>
0154 struct Measure<DEGREE, 1> {
0155   inline void operator()(unsigned long long& t) const {
0156     t -= rdtsc();
0157     justcomp::bar<DEGREE>();
0158     t += rdtsc();
0159   }
0160 };
0161
0162 template <int DEGREE>
0163 struct Measure<DEGREE, 2> {
0164   inline void operator()(unsigned long long& t) const {
0165     t -= rdtsc();
0166     justcomp::loo<DEGREE>();
0167     t += rdtsc();
0168   }
0169 };
0170
0171 template <int DEGREE>
0172 struct Measure<DEGREE, 3> {
0173   inline void operator()(unsigned long long& t) const {
0174     t -= rdtsc();
0175     justcomp::lar<DEGREE>();
0176     t += rdtsc();
0177   }
0178 };
0179
0180 template <int DEGREE, int WHAT = 1>
0181 void perf() {
0182   if (!has_rdtscp())
0183     return;
0184   Measure<DEGREE, WHAT> measure;
0185   using namespace approx_math;
0186   unsigned long long t = 0;
0187   binary32 x;
0188   float sum = 0;
0189   long long ntot = 0;
0190   x.f = 1.0;          // should be 0 but
0191   while (x.f < 32) {  // this is 5*2^23 tests
0192     ++ntot;
0193     int i = 0;
0194     while (i < justcomp::NN) {
0195       x.ui32++;
0196       justcomp::a[i++] = x.f;
0197       justcomp::a[i++] = (WHAT < 2) ? -x.f : 1.f / x.f;
0198     }
0199     measure(t);
0200     // binary32 r;
0201     //  r.f=approx_expf<6>(x.f);// time 0m1.180s
0202     // r.f=expf(x.f);   // time 0m4.372s
0203     // r.f=exp(x.f);  // time   0m1.789s
0204     for (int i = 0; i != justcomp::NN; ++i)
0205       sum += justcomp::b[i];
0206   }
0207   const char* what[] = {"exp unsafe", "exp apporx", "log unsafe", "log approx"};
0208   std::cout << "time for " << what[WHAT] << " degree " << DEGREE << " is " << double(t) / double(justcomp::NN * ntot)
0209             << std::endl;
0210   std::cout << "sum= " << sum << " to prevent compiler optim." << std::endl;
0211   ;
0212 }
0213
0214 int main() {
0215   printLim(printEXP);
0216   printLim(printLOG);
0217
0218   perf<2, 0>();
0219   perf<3, 0>();
0220   perf<3>();
0221   perf<4>();
0222   perf<6>();
0223
0224   perf<2, 2>();
0225   perf<2, 3>();
0226   perf<4, 2>();
0227   perf<4, 3>();
0228   perf<8, 2>();
0229   perf<8, 3>();
0230
0231   accTest(::exp, approx_expf<2>, 2);
0232   accTest(::exp, approx_expf<3>, 3);
0233   accTest(::exp, approx_expf<4>, 4);
0234   accTest(::exp, approx_expf<6>, 6);
0235
0236   accTest(::log, approx_logf<2>, 2);
0237   accTest(::log, approx_logf<4>, 4);
0238   accTest(::log, approx_logf<8>, 8);
0239
0240   accuTest(::erf, approx_erf, "erf", .01, 8);
0241
0242   return 0;
0243 }