Math/test/testApproximations.cpp

0001 #include <iostream>
0002 #include <sstream>
0003 #include <fstream>
0004 #include <string>
0005 #include <vector>
0006 #include <cstdlib>
0007 #include <cstdio>
0008 #include <cmath>
0009
0010 #include "DataFormats/Math/interface/approx_exp.h"
0011 #include "DataFormats/Math/interface/approx_log.h"
0012
0013 inline int diff(float a, float b) {
0014   approx_math::binary32 ba(a);
0015   approx_math::binary32 bb(b);
0016   return ba.i32 - bb.i32;
0017 }
0018
0019 inline int bits(int a) {
0020   unsigned int aa = abs(a);
0021   int b = 0;
0022   if (a == 0)
0023     return 0;
0024   while ((aa /= 2) > 0)
0025     ++b;
0026   return (a > 0) ? b : -b;
0027 }
0028
0029 void testIt() {
0030   float y[10];
0031   int a[10]{0}, h[10]{0}, l[10]{9999999};
0032   const int N = 1000000;
0033   for (int i = 0; i != N; ++i) {
0034     float x = 1.e-9 + 1.e9 * drand48();
0035     y[0] = logf(x);
0036     y[2] = unsafe_logf<2>(x);
0037     y[3] = unsafe_logf<3>(x);
0038     y[4] = unsafe_logf<4>(x);
0039     y[5] = unsafe_logf<5>(x);
0040     y[6] = unsafe_logf<6>(x);
0041     y[7] = unsafe_logf<7>(x);
0042     y[8] = unsafe_logf<8>(x);
0043     for (int k = 2; k != 9; k++) {
0044       a[k] += diff(y[0], y[k]);
0045       h[k] = std::max(h[k], diff(y[0], y[k]));
0046       l[k] = std::min(l[k], diff(y[0], y[k]));
0047     }
0048   }
0049   for (int k = 2; k != 9; k++) {
0050     std::cout << k << ": ave/min/max " << double(a[k]) / double(N) << " " << l[k] << " " << h[k] << std::endl;
0051   }
0052 }
0053
0054 inline float ms(float radLen, float m2, float p2) {
0055   constexpr float amscon = 1.8496e-4;  // (13.6MeV)**2
0056   float e2 = p2 + m2;
0057
0058   float fact = 1.f + 0.038f * log(radLen);
0059   fact /= p2;
0060   fact *= fact;
0061   float a = e2 * fact;
0062   return amscon * radLen * a;
0063 }
0064
0065 inline float msf(float radLen, float m2, float p2) {
0066   constexpr float amscon = 1.8496e-4;  // (13.6MeV)**2
0067   float e2 = p2 + m2;
0068
0069   float fact = 1.f + 0.038f * unsafe_logf<2>(radLen);
0070   fact /= p2;
0071   fact *= fact;
0072   float a = e2 * fact;
0073   return amscon * radLen * a;
0074 }
0075
0076 inline float ms2(float radLen, float m2, float p2) {
0077   constexpr float amscon = 1.8496e-4;  // (13.6MeV)**2
0078   float e2 = p2 + m2;
0079   float beta2 = p2 / e2;
0080   float fact = 1.f + 0.038f * log(radLen);
0081   fact *= fact;
0082   float a = fact / (beta2 * p2);
0083   return amscon * radLen * a;
0084 }
0085
0086 template <typename T>
0087 inline float bb2(float xi, float m2, float p2) {
0088   const T emass = 0.511e-3;
0089   const T poti = 16.e-9 * 10.75;                     // = 16 eV * Z**0.9, for Si Z=14
0090   const T eplasma = 28.816e-9 * sqrt(2.33 * 0.498);  // 28.816 eV * sqrt(rho*(Z/A)) for Si
0091   const T delta0 = 2 * log(eplasma / poti) - 1.;
0092
0093   // calculate general physics things
0094   T p = sqrt(p2);
0095   T m = sqrt(m2);
0096   T e = sqrt(p2 + m2);
0097   T beta = p / e;
0098   T gamma = e / m;
0099   T eta2 = beta * gamma;
0100   eta2 *= eta2;
0101   T ratio = emass / m;
0102   T emax = 2. * emass * eta2 / (1. + 2. * ratio * gamma + ratio * ratio);
0103
0104   xi /= (beta * beta);
0105
0106   return xi * (log(2. * emass * emax / (poti * poti)) - 2. * (beta * beta) - delta0);
0107 }
0108
0109 template <typename T>
0110 inline float bb(float xi, float m2, float p2) {
0111   const T emass = 0.511e-3;
0112   const T poti = 16.e-9 * 10.75;                     // = 16 eV * Z**0.9, for Si Z=14
0113   const T eplasma = 28.816e-9 * sqrt(2.33 * 0.498);  // 28.816 eV * sqrt(rho*(Z/A)) for Si
0114   const T delta0 = 2 * log(eplasma / poti) - 1.;
0115
0116   // calculate general physics things
0117   T im2 = T(1.) / m2;
0118   T e2 = p2 + m2;
0119   T e = sqrt(e2);
0120   T beta2 = p2 / e2;
0121   T eta2 = p2 * im2;
0122   T ratio2 = (emass * emass) * im2;
0123   T emax = T(2.) * emass * eta2 / (T(1.) + T(2.) * emass * e * im2 + ratio2);
0124
0125   xi /= beta2;
0126
0127   return xi * (log(T(2.) * emass * emax / (poti * poti)) - T(2.) * (beta2)-delta0);
0128 }
0129
0130 template <typename T>
0131 inline float bbf(float xi, float m2, float p2) {
0132   const T emass = 0.511e-3;
0133   const T poti = 16.e-9 * 10.75;                     // = 16 eV * Z**0.9, for Si Z=14
0134   const T eplasma = 28.816e-9 * sqrt(2.33 * 0.498);  // 28.816 eV * sqrt(rho*(Z/A)) for Si
0135   const T delta0 = 2 * log(eplasma / poti) - 1.;
0136
0137   // calculate general physics things
0138   T im2 = T(1.) / m2;
0139   T e2 = p2 + m2;
0140   T e = sqrt(e2);
0141   T beta2 = p2 / e2;
0142   T eta2 = p2 * im2;
0143   T ratio2 = (emass * emass) * im2;
0144   T emax = T(2.) * emass * eta2 / (T(1.) + T(2.) * emass * e * im2 + ratio2);
0145
0146   xi /= beta2;
0147
0148   return xi * (unsafe_logf<2>(T(2.) * emass * emax / (poti * poti)) - T(2.) * (beta2)-delta0);
0149 }
0150
0151 template <typename T>
0152 inline float bbf2(float xi, float m2, float p2) {
0153   const T emass = 0.511e-3;
0154   const T poti = 16.e-9 * 10.75;                     // = 16 eV * Z**0.9, for Si Z=14
0155   const T eplasma = 28.816e-9 * sqrt(2.33 * 0.498);  // 28.816 eV * sqrt(rho*(Z/A)) for Si
0156   const T delta0 = 2 * log(eplasma / poti) - 1.;
0157
0158   // calculate general physics things
0159   T im2 = T(1.) / m2;
0160   T e2 = p2;  //  + m2;
0161   T e = sqrt(e2);
0162   T beta2 = T(1);  //  p2/e2;
0163   T eta2 = p2 * im2;
0164   T ratio2 = (emass * emass) * im2;
0165   T emax = T(2.) * emass * eta2 / (T(1.) + T(2.) * emass * e * im2 + ratio2);
0166
0167   xi /= beta2;
0168
0169   return xi * (unsafe_logf<2>(T(2.) * emass * emax / (poti * poti)) - T(2.) * (beta2)-delta0);
0170 }
0171
0172 template <typename Fun>
0173 void compare(Fun F, Fun F2, Fun Fapx) {
0174   std::cout << std::endl;
0175
0176   float m2 = 0.138;
0177   m2 *= m2;
0178
0179   int d1 = 0, d2 = 0, d3 = 0;
0180   int c1 = 99999999, c2 = c1, c3 = c1;
0181   int dm = 99999999;
0182
0183   float p2 = 0.01;
0184   for (int i = 0; i != 6; ++i) {
0185     p2 *= 10;
0186     float rl = 0.001;
0187     for (int j = 0; j != 4; ++j) {
0188       rl *= 10;
0189       float ref = F(rl, m2, p2);
0190       float rp = F(rl * 1.001, m2, p2);
0191       float rm = F(rl * 0.999, m2, p2);
0192       float apx = Fapx(rl, m2, p2);
0193
0194       int dd = std::min(abs(diff(rm, ref)), abs(diff(rp, ref)));
0195       dd -= abs(diff(apx, ref));  // negative if apx-ref is bigger than the uncer-interval
0196       dm = std::min(dm, dd);
0197
0198       d1 = std::max(d1, abs(diff(F2(rl, m2, p2), ref)));
0199       d2 = std::max(d2, abs(diff(apx, ref)));
0200       d3 = std::max(d3, abs(diff(rp, ref)));
0201       d3 = std::max(d3, abs(diff(rm, ref)));
0202
0203       c1 = std::min(c1, abs(diff(F2(rl, m2, p2), ref)));
0204       c2 = std::min(c2, abs(diff(apx, ref)));
0205       c3 = std::min(c3, abs(diff(rp, ref)));
0206       c3 = std::min(c3, abs(diff(rm, ref)));
0207
0208       // std::cout << diff(ms2(rl,m2,p2),ref) << std::endl;
0209       // std::cout << diff(msf(rl,m2,p2),ref) << std::endl;
0210       // std::cout << diff(ms(1.001*rl,m2,p2),ref) << std::endl;
0211       // std::cout << diff(ms(0.999*rl,m2,p2),ref) << std::endl;
0212     }
0213   }
0214
0215   std::cout << dm << "," << bits(dm) << std::endl;
0216
0217   std::cout << d1 << "," << bits(d1) << " " << d2 << "," << bits(d2) << " " << d3 << "," << bits(d3) << " "
0218             << std::endl;
0219
0220   std::cout << c1 << "," << bits(c1) << " " << c2 << "," << bits(c2) << " " << c3 << "," << bits(c3) << " "
0221             << std::endl;
0222 }
0223
0224 int main() {
0225   std::cout << bits(0) << " " << bits(1) << " " << bits(2) << " " << bits(-31) << " " << bits(32) << std::endl;
0226
0227   testIt();
0228   compare(ms, ms2, msf);
0229   compare(bb<float>, bb2<float>, bbf<float>);
0230   compare(bb<float>, bb2<float>, bbf2<float>);
0231   compare(bb<double>, bb2<double>, bbf<double>);
0232
0233   return 0;
0234 }