Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:54:04

0001 #pragma GCC diagnostic ignored "-Wformat"
0002 #include <iostream>
0003 #include <cmath>
0004 #include <cstdlib>
0005 #include <cstdio>
0006 
0007 namespace almostEqualDetail {
0008   union fasi {
0009     int i;
0010     float f;
0011   };
0012   
0013   union dasi {
0014     long long i;
0015     double f;
0016   };
0017   
0018 }
0019 
0020 inline int intDiff(float a, float b)
0021 {
0022   using namespace almostEqualDetail;
0023   // Make sure maxUlps is non-negative and small enough that the
0024   // default NAN won't compare as equal to anything.
0025   fasi fa; fa.f = a;
0026   // Make aInt lexicographically ordered as a twos-complement int
0027   if (fa.i < 0) fa.i = 0x80000000 - fa.i;
0028   // Make bInt lexicographically ordered as a twos-complement int
0029   fasi fb; fb.f = b;
0030   if (fb.i < 0) fb.i = 0x80000000 - fb.i;
0031   return std::abs(fa.i - fb.i);
0032 }
0033 
0034 
0035 inline long long intDiff(double a, double b)
0036 {
0037   using namespace almostEqualDetail;
0038   dasi fa; fa.f = a;
0039   // Make aInt lexicographically ordered as a twos-complement int
0040   if (fa.i < 0) fa.i = 0x8000000000000000LL - fa.i;
0041   // Make bInt lexicographically ordered as a twos-complement int
0042   dasi fb; fb.f = b;
0043   if (fb.i < 0) fb.i = 0x8000000000000000LL - fb.i;
0044   return std::abs(fa.i - fb.i);
0045 }
0046 
0047 template<typename T>
0048 inline bool almostEqual(T a, T b, int maxUlps) {
0049   // Make sure maxUlps is non-negative and small enough that the
0050   // default NAN won't compare as equal to anything.
0051   // assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024);
0052   return intDiff(a,b) <= maxUlps;
0053 }
0054 
0055 
0056 namespace {
0057   template<typename T>
0058   inline T eta(T x, T y, T z) { T t(z/std::sqrt(x*x+y*y)); return std::log(t+std::sqrt(t*t+T(1)));} 
0059   template<typename T>
0060   inline T eta2(T x, T y, T z) { T t = (z*z)/(x*x+y*y); return copysign(std::log(std::sqrt(t)+std::sqrt(t+T(1))), z); }
0061 
0062   inline float eta3(float x, float y, float z) { float t(z/std::sqrt(x*x+y*y)); return ::asinhf(t);} 
0063   inline double eta3(double x, double y, double z) { double t(z/std::sqrt(x*x+y*y)); return ::asinh(t);} 
0064 
0065   
0066   void look(float x) {
0067     int e;
0068     float r = ::frexpf(x,&e);
0069     std::cout << x << " exp " << e << " res " << r << std::endl;
0070     
0071     union {
0072       float val;
0073       int bin;
0074     } f;
0075     
0076     f.val = x;
0077     printf("%e %a %x\n",  f.val,  f.val,  f.bin);
0078     // printf("%e %x\n", f.val,  f.bin);
0079     int log_2 = ((f.bin >> 23) & 255) - 127;  //exponent
0080     f.bin &= 0x7FFFFF;                               //mantissa (aka significand)
0081     
0082     std::cout << "exp " << log_2 << " mant in binary "  << std::hex << f.bin 
0083           << " mant as float " <<  std::dec << (f.bin|0x800000)*::pow(2.,-23)
0084           << std::endl << std::endl;
0085   }
0086 
0087   void look(double x) {
0088     // int e;
0089     // float r = ::frexpf(x,&e);
0090     // std::cout << x << " exp " << e << " res " << r << std::endl;
0091     
0092     union {
0093       double val;
0094       long long bin;
0095     } f;
0096     
0097     f.val = x;
0098 
0099     printf("%e %a %x\n",  f.val,  f.val,  f.bin);
0100     // printf("%e %x\n", f.val,  f.bin);
0101 
0102   }
0103 
0104 }
0105 
0106 void peta() {
0107   std::cout << "T t(z/std::sqrt(x*x+y*y)); return std::log(t+std::sqrt(t*t+T(1)));" << std::endl;
0108   {
0109     float 
0110       xn = 122.436f, yn = 10.7118f, zn = -1115.f;
0111     float etan = eta(xn,yn,zn);
0112  
0113     std::cout << etan << std::endl;
0114     look(etan);
0115     std::cout << -etan << std::endl;
0116     look(-etan);
0117 
0118    float 
0119       xp = 122.436f, yp = 10.7118f, zp = 1115.f;
0120     float etap = eta(xp,yp,zp);
0121  
0122 
0123     std::cout << etap << std::endl;
0124     look(etap);
0125 
0126     std::cout << intDiff(etap,-etan) << std::endl << std::endl;
0127 
0128   }
0129 
0130  {
0131     double 
0132       xn = 122.436, yn = 10.7118, zn = -1115.;
0133     double etan = eta(xn,yn,zn);
0134  
0135     std::cout << etan << std::endl;
0136     look(etan);
0137     std::cout << -etan << std::endl;
0138     look(-etan);
0139 
0140    double 
0141       xp = 122.436, yp = 10.7118, zp = 1115.;
0142     double etap = eta(xp,yp,zp);
0143  
0144 
0145     std::cout << etap << std::endl;
0146     look(etap);
0147 
0148     std::cout << intDiff(etap,-etan) << std::endl << std::endl;
0149 
0150   }
0151 
0152 }
0153 
0154 
0155 void peta2() {
0156   std::cout << "T t = (z*z)/(x*x+y*y); return copysign(std::log(std::sqrt(t)+std::sqrt(t+T(1))), z);"<< std::endl;
0157   {
0158     float 
0159       xn = 122.436f, yn = 10.7118f, zn = -1115.f;
0160     float etan = eta2(xn,yn,zn);
0161  
0162     std::cout << etan << std::endl;
0163     look(etan);
0164     std::cout << -etan << std::endl;
0165     look(-etan);
0166 
0167    float 
0168       xp = 122.436f, yp = 10.7118f, zp = 1115.f;
0169     float etap = eta2(xp,yp,zp);
0170  
0171 
0172     std::cout << etap << std::endl;
0173     look(etap);
0174 
0175     std::cout << intDiff(etap,-etan) << std::endl << std::endl;
0176 
0177   }
0178 
0179  {
0180     double 
0181       xn = 122.436, yn = 10.7118, zn = -1115.;
0182     double etan = eta2(xn,yn,zn);
0183  
0184     std::cout << etan << std::endl;
0185     look(etan);
0186     std::cout << -etan << std::endl;
0187     look(-etan);
0188 
0189    double 
0190       xp = 122.436, yp = 10.7118, zp = 1115.;
0191     double etap = eta2(xp,yp,zp);
0192  
0193 
0194     std::cout << etap << std::endl;
0195     look(etap);
0196 
0197     std::cout << intDiff(etap,-etan) << std::endl << std::endl;
0198 
0199   }
0200 
0201 }
0202 
0203 void peta3() {
0204   std::cout << "t(z/std::sqrt(x*x+y*y)); return ::asinh(t);" << std::endl;
0205   {
0206     float 
0207       xn = 122.436f, yn = 10.7118f, zn = -1115.f;
0208     float etan = eta3(xn,yn,zn);
0209  
0210     std::cout << etan << std::endl;
0211     look(etan);
0212     std::cout << -etan << std::endl;
0213     look(-etan);
0214 
0215    float 
0216       xp = 122.436f, yp = 10.7118f, zp = 1115.f;
0217     float etap = eta3(xp,yp,zp);
0218  
0219 
0220     std::cout << etap << std::endl;
0221     look(etap);
0222 
0223     std::cout << intDiff(etap,-etan) << std::endl << std::endl;
0224 
0225   }
0226 
0227  {
0228     double 
0229       xn = 122.436, yn = 10.7118, zn = -1115.;
0230     double etan = eta3(xn,yn,zn);
0231  
0232     std::cout << etan << std::endl;
0233     look(etan);
0234     std::cout << -etan << std::endl;
0235     look(-etan);
0236 
0237    double 
0238       xp = 122.436, yp = 10.7118, zp = 1115.;
0239     double etap = eta3(xp,yp,zp);
0240  
0241 
0242     std::cout << etap << std::endl;
0243     look(etap);
0244 
0245     std::cout << intDiff(etap,-etan) << std::endl << std::endl;
0246 
0247   }
0248 
0249 }
0250 
0251 int main() {
0252 
0253   peta();
0254   peta2();
0255   peta3();
0256 
0257     return 0;
0258 
0259 }