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