Math/test/eta_t.cpp

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 }