Math/interface/approx_atan2.h

0001 #ifndef DataFormatsMathAPPROX_ATAN2_H
0002 #define DataFormatsMathAPPROX_ATAN2_H
0003
0004 /*
0005  * approximate atan2 evaluations
0006  *
0007  * Polynomials were obtained using Sollya scripts (in comments below)
0008  *
0009  *
0010 */
0011
0012 /*
0013 f= atan((1-x)/(1+x))-atan(1);
0014 I=[-1+10^(-4);1.0];
0015 filename="atan.txt";
0016 print("") > filename;
0017 for deg from 3 to 11 do begin
0018   p = fpminimax(f, deg,[|1,23...|],I, floating, absolute);
0019   display=decimal;
0020   acc=floor(-log2(sup(supnorm(p, f, I, absolute, 2^(-20)))));
0021   print( "   // degree = ", deg,
0022          "  => absolute accuracy is ",  acc, "bits" ) >> filename;
0023   print("template<> constexpr float approx_atan2f_P<", deg, ">(float x){") >> filename;
0024   display=hexadecimal;
0025   print(" return ", horner(p) , ";") >> filename;
0026   print("}") >> filename;
0027 end;
0028 */
0029
0030 #include <cstdint>
0031 #include <cmath>
0032 #include <limits>
0033 #include <algorithm>
0034
0035 // float
0036
0037 template <int DEGREE>
0038 constexpr float approx_atan2f_P(float x);
0039
0040 // degree =  3   => absolute accuracy is  7 bits
0041 template <>
0042 constexpr float approx_atan2f_P<3>(float x) {
0043   return x * (float(-0xf.8eed2p-4) + x * x * float(0x3.1238p-4));
0044 }
0045
0046 // degree =  5   => absolute accuracy is  10 bits
0047 template <>
0048 constexpr float approx_atan2f_P<5>(float x) {
0049   auto z = x * x;
0050   return x * (float(-0xf.ecfc8p-4) + z * (float(0x4.9e79dp-4) + z * float(-0x1.44f924p-4)));
0051 }
0052
0053 // degree =  7   => absolute accuracy is  13 bits
0054 template <>
0055 constexpr float approx_atan2f_P<7>(float x) {
0056   auto z = x * x;
0057   return x * (float(-0xf.fcc7ap-4) + z * (float(0x5.23886p-4) + z * (float(-0x2.571968p-4) + z * float(0x9.fb05p-8))));
0058 }
0059
0060 // degree =  9   => absolute accuracy is  16 bits
0061 template <>
0062 constexpr float approx_atan2f_P<9>(float x) {
0063   auto z = x * x;
0064   return x * (float(-0xf.ff73ep-4) +
0065               z * (float(0x5.48ee1p-4) +
0066                    z * (float(-0x2.e1efe8p-4) + z * (float(0x1.5cce54p-4) + z * float(-0x5.56245p-8)))));
0067 }
0068
0069 // degree =  11   => absolute accuracy is  19 bits
0070 template <>
0071 constexpr float approx_atan2f_P<11>(float x) {
0072   auto z = x * x;
0073   return x * (float(-0xf.ffe82p-4) +
0074               z * (float(0x5.526c8p-4) +
0075                    z * (float(-0x3.18bea8p-4) +
0076                         z * (float(0x1.dce3bcp-4) + z * (float(-0xd.7a64ap-8) + z * float(0x3.000eap-8))))));
0077 }
0078
0079 // degree =  13   => absolute accuracy is  21 bits
0080 template <>
0081 constexpr float approx_atan2f_P<13>(float x) {
0082   auto z = x * x;
0083   return x * (float(-0xf.fffbep-4) +
0084               z * (float(0x5.54adp-4) +
0085                    z * (float(-0x3.2b4df8p-4) +
0086                         z * (float(0x2.1df79p-4) +
0087                              z * (float(-0x1.46081p-4) + z * (float(0x8.99028p-8) + z * float(-0x1.be0bc4p-8)))))));
0088 }
0089
0090 // degree =  15   => absolute accuracy is  24 bits
0091 template <>
0092 constexpr float approx_atan2f_P<15>(float x) {
0093   auto z = x * x;
0094   return x * (float(-0xf.ffff4p-4) +
0095               z * (float(0x5.552f9p-4 + z * (float(-0x3.30f728p-4) +
0096                                              z * (float(0x2.39826p-4) +
0097                                                   z * (float(-0x1.8a880cp-4) +
0098                                                        z * (float(0xe.484d6p-8) +
0099                                                             z * (float(-0x5.93d5p-8) + z * float(0x1.0875dcp-8)))))))));
0100 }
0101
0102 template <int DEGREE>
0103 constexpr float unsafe_atan2f_impl(float y, float x) {
0104   constexpr float pi4f = 3.1415926535897932384626434 / 4;
0105   constexpr float pi34f = 3.1415926535897932384626434 * 3 / 4;
0106
0107   auto r = (std::abs(x) - std::abs(y)) / (std::abs(x) + std::abs(y));
0108   if (x < 0)
0109     r = -r;
0110
0111   auto angle = (x >= 0) ? pi4f : pi34f;
0112   angle += approx_atan2f_P<DEGREE>(r);
0113
0114   return ((y < 0)) ? -angle : angle;
0115 }
0116
0117 template <int DEGREE>
0118 constexpr float unsafe_atan2f(float y, float x) {
0119   return unsafe_atan2f_impl<DEGREE>(y, x);
0120 }
0121
0122 template <int DEGREE>
0123 constexpr float safe_atan2f(float y, float x) {
0124   return unsafe_atan2f_impl<DEGREE>(y, ((y == 0.f) & (x == 0.f)) ? 0.2f : x);
0125   // return (y==0.f)&(x==0.f) ? 0.f :  unsafe_atan2f_impl<DEGREE>( y, x);
0126 }
0127
0128 // integer...
0129 /*
0130   f= (2^31/pi)*(atan((1-x)/(1+x))-atan(1));
0131   I=[-1+10^(-4);1.0];
0132   p = fpminimax(f, [|1,3,5,7,9,11|],[|23...|],I, floating, absolute);
0133  */
0134
0135 template <int DEGREE>
0136 constexpr float approx_atan2i_P(float x);
0137
0138 // degree =  3   => absolute accuracy is  6*10^6
0139 template <>
0140 constexpr float approx_atan2i_P<3>(float x) {
0141   auto z = x * x;
0142   return x * (-664694912.f + z * 131209024.f);
0143 }
0144
0145 // degree =  5   => absolute accuracy is  4*10^5
0146 template <>
0147 constexpr float approx_atan2i_P<5>(float x) {
0148   auto z = x * x;
0149   return x * (-680392064.f + z * (197338400.f + z * (-54233256.f)));
0150 }
0151
0152 // degree =  7   => absolute accuracy is  6*10^4
0153 template <>
0154 constexpr float approx_atan2i_P<7>(float x) {
0155   auto z = x * x;
0156   return x * (-683027840.f + z * (219543904.f + z * (-99981040.f + z * 26649684.f)));
0157 }
0158
0159 // degree =  9   => absolute accuracy is  8000
0160 template <>
0161 constexpr float approx_atan2i_P<9>(float x) {
0162   auto z = x * x;
0163   return x * (-683473920.f + z * (225785056.f + z * (-123151184.f + z * (58210592.f + z * (-14249276.f)))));
0164 }
0165
0166 // degree =  11   => absolute accuracy is  1000
0167 template <>
0168 constexpr float approx_atan2i_P<11>(float x) {
0169   auto z = x * x;
0170   return x *
0171          (-683549696.f + z * (227369312.f + z * (-132297008.f + z * (79584144.f + z * (-35987016.f + z * 8010488.f)))));
0172 }
0173
0174 // degree =  13   => absolute accuracy is  163
0175 template <>
0176 constexpr float approx_atan2i_P<13>(float x) {
0177   auto z = x * x;
0178   return x * (-683562624.f +
0179               z * (227746080.f +
0180                    z * (-135400128.f + z * (90460848.f + z * (-54431464.f + z * (22973256.f + z * (-4657049.f)))))));
0181 }
0182
0183 template <>
0184 constexpr float approx_atan2i_P<15>(float x) {
0185   auto z = x * x;
0186   return x * (-683562624.f +
0187               z * (227746080.f +
0188                    z * (-135400128.f + z * (90460848.f + z * (-54431464.f + z * (22973256.f + z * (-4657049.f)))))));
0189 }
0190
0191 template <int DEGREE>
0192 constexpr int unsafe_atan2i_impl(float y, float x) {
0193   constexpr long long maxint = (long long)(std::numeric_limits<int>::max()) + 1LL;
0194   constexpr int pi4 = int(maxint / 4LL);
0195   constexpr int pi34 = int(3LL * maxint / 4LL);
0196
0197   auto r = (std::abs(x) - std::abs(y)) / (std::abs(x) + std::abs(y));
0198   if (x < 0)
0199     r = -r;
0200
0201   auto angle = (x >= 0) ? pi4 : pi34;
0202   angle += int(approx_atan2i_P<DEGREE>(r));
0203   // angle += int(std::round(approx_atan2i_P<DEGREE>(r)));
0204
0205   return (y < 0) ? -angle : angle;
0206 }
0207
0208 template <int DEGREE>
0209 constexpr int unsafe_atan2i(float y, float x) {
0210   return unsafe_atan2i_impl<DEGREE>(y, x);
0211 }
0212
0213 // short (16bits)
0214
0215 template <int DEGREE>
0216 constexpr float approx_atan2s_P(float x);
0217
0218 // degree =  3   => absolute accuracy is  53
0219 template <>
0220 constexpr float approx_atan2s_P<3>(float x) {
0221   auto z = x * x;
0222   return x * ((-10142.439453125f) + z * 2002.0908203125f);
0223 }
0224 // degree =  5   => absolute accuracy is  7
0225 template <>
0226 constexpr float approx_atan2s_P<5>(float x) {
0227   auto z = x * x;
0228   return x * ((-10381.9609375f) + z * ((3011.1513671875f) + z * (-827.538330078125f)));
0229 }
0230 // degree =  7   => absolute accuracy is  2
0231 template <>
0232 constexpr float approx_atan2s_P<7>(float x) {
0233   auto z = x * x;
0234   return x * ((-10422.177734375f) + z * (3349.97412109375f + z * ((-1525.589599609375f) + z * 406.64190673828125f)));
0235 }
0236 // degree =  9   => absolute accuracy is 1
0237 template <>
0238 constexpr float approx_atan2s_P<9>(float x) {
0239   auto z = x * x;
0240   return x * ((-10428.984375f) + z * (3445.20654296875f + z * ((-1879.137939453125f) +
0241                                                                z * (888.22314453125f + z * (-217.42669677734375f)))));
0242 }
0243
0244 template <int DEGREE>
0245 constexpr short unsafe_atan2s_impl(float y, float x) {
0246   constexpr int maxshort = (int)(std::numeric_limits<short>::max()) + 1;
0247   constexpr short pi4 = short(maxshort / 4);
0248   constexpr short pi34 = short(3 * maxshort / 4);
0249
0250   auto r = (std::abs(x) - std::abs(y)) / (std::abs(x) + std::abs(y));
0251   if (x < 0)
0252     r = -r;
0253
0254   auto angle = (x >= 0) ? pi4 : pi34;
0255   angle += short(approx_atan2s_P<DEGREE>(r));
0256
0257   return (y < 0) ? -angle : angle;
0258 }
0259
0260 template <int DEGREE>
0261 constexpr short unsafe_atan2s(float y, float x) {
0262   return unsafe_atan2s_impl<DEGREE>(y, x);
0263 }
0264
0265 constexpr int phi2int(float x) {
0266   constexpr float p2i = ((long long)(std::numeric_limits<int>::max()) + 1LL) / M_PI;
0267   return std::round(x * p2i);
0268 }
0269
0270 constexpr float int2phi(int x) {
0271   constexpr float i2p = M_PI / ((long long)(std::numeric_limits<int>::max()) + 1LL);
0272   return float(x) * i2p;
0273 }
0274
0275 constexpr double int2dphi(int x) {
0276   constexpr double i2p = M_PI / ((long long)(std::numeric_limits<int>::max()) + 1LL);
0277   return x * i2p;
0278 }
0279
0280 constexpr short phi2short(float x) {
0281   constexpr float p2i = ((int)(std::numeric_limits<short>::max()) + 1) / M_PI;
0282   return std::round(x * p2i);
0283 }
0284
0285 constexpr float short2phi(short x) {
0286   constexpr float i2p = M_PI / ((int)(std::numeric_limits<short>::max()) + 1);
0287   return float(x) * i2p;
0288 }
0289
0290 #endif