Math/interface/approx_asin.h

0001 #ifndef DataFormatsMathapprox_asin_h
0002 #define DataFormatsMathapprox_asin_h
0003
0004 #include <cmath>
0005 template <int DEGREE>
0006 constexpr float approx_asin_P(float z);
0007
0008 // degree =  3   => absolute accuracy is  8 bits
0009 template <>
0010 constexpr float approx_asin_P<3>(float z) {
0011   return 1.f + z * 0.2133004963397979736328125f;
0012 }
0013 // degree =  5   => absolute accuracy is  11 bits
0014 template <>
0015 constexpr float approx_asin_P<5>(float z) {
0016   return 1.f + z * (0.154711246490478515625f + z * 0.1322681009769439697265625f);
0017 }
0018 // degree =  7   => absolute accuracy is  15 bits
0019 template <>
0020 constexpr float approx_asin_P<7>(float z) {
0021   return 1.f + z * (0.169519245624542236328125f + z * (4.9031913280487060546875e-2f + z * 0.10941398143768310546875));
0022 }
0023 // degree =  9   => absolute accuracy is  18 bits
0024 template <>
0025 constexpr float approx_asin_P<9>(float z) {
0026   return 1.f + z * (0.166020572185516357421875f +
0027                     z * (8.44048559665679931640625e-2f +
0028                          z * (1.11602735705673694610595703125e-3f + z * 0.103476583957672119140625f)));
0029 }
0030 // degree =  11   => absolute accuracy is  21 bits
0031 template <>
0032 constexpr float approx_asin_P<11>(float z) {
0033   return 1.f + z * (0.1668075025081634521484375f +
0034                     z * (7.20207393169403076171875e-2f +
0035                          z * (6.607978045940399169921875e-2f +
0036                               z * ((-3.6048568785190582275390625e-2f) + z * 0.10574872791767120361328125f))));
0037 }
0038
0039 /*
0040    // degree =  3   => absolute accuracy is  8 bits
0041 template<> constexpr float approx_asin_P< 3 >(float z){
0042  return  1.f + z * 0.2114248573780059814453125f;
0043 }
0044    // degree =  5   => absolute accuracy is  12 bits
0045 template<> constexpr float approx_asin_P< 5 >(float z){
0046  return  1.f + z * (0.1556626856327056884765625f + z * 0.1295671761035919189453125f);
0047 }
0048    // degree =  7   => absolute accuracy is  15 bits
0049 template<> constexpr float approx_asin_P< 7 >(float z){
0050  return  1.f + z * (0.1691854894161224365234375f + z * (5.1305986940860748291015625e-2f + z * 0.1058919131755828857421875f));
0051 }
0052    // degree =  9   => absolute accuracy is  18 bits
0053 template<> constexpr float approx_asin_P< 9 >(float z){
0054  return  1.f + z * (0.166119158267974853515625f + z * (8.322779834270477294921875e-2f + z * (5.28236292302608489990234375e-3f + z * 9.89462435245513916015625e-2f)));
0055 }
0056    // degree =  11   => absolute accuracy is  21 bits
0057 template<> constexpr float approx_asin_P< 11 >(float z){
0058  return  1.f + z * (0.1667812168598175048828125f + z * (7.249967753887176513671875e-2f + z * (6.321799755096435546875e-2f + z * ((-2.913488447666168212890625e-2f) + z * 9.9913299083709716796875e-2f))));
0059 }
0060 */
0061
0062 // valid for |x|<0.71
0063 template <int DEGREE>
0064 constexpr float unsafe_asin07(float x) {
0065   auto z = x * x;
0066   return x * approx_asin_P<DEGREE>(z);
0067 }
0068
0069 template <int DEGREE>
0070 constexpr float unsafe_acos07(float x) {
0071   constexpr float pihalf = M_PI / 2;
0072   return pihalf - unsafe_asin07<DEGREE>(x);
0073 }
0074
0075 // for |x|> 0.71 use slower
0076 template <int DEGREE>
0077 constexpr float unsafe_acos71(float x) {
0078   constexpr float pi = M_PI;
0079   auto z = 1.f - x * x;
0080   z = std::sqrt(z) * approx_asin_P<DEGREE>(z);
0081   return x > 0 ? z : pi - z;
0082 }
0083
0084 template <int DEGREE>
0085 constexpr float unsafe_asin71(float x) {
0086   constexpr float pihalf = M_PI / 2;
0087   return pihalf - unsafe_acos71<DEGREE>(x);
0088 }
0089
0090 template <int DEGREE>
0091 constexpr float unsafe_asin(float x) {
0092   return (std::abs(x) < 0.71f) ? unsafe_asin07<DEGREE>(x) : unsafe_asin71<DEGREE>(x);
0093 }
0094
0095 template <int DEGREE>
0096 constexpr float unsafe_acos(float x) {
0097   return (std::abs(x) < 0.71f) ? unsafe_acos07<DEGREE>(x) : unsafe_acos71<DEGREE>(x);
0098 }
0099
0100 #endif