SqrtEstimate.h 5.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161
  1. // David Eberly, Geometric Tools, Redmond WA 98052
  2. // Copyright (c) 1998-2020
  3. // Distributed under the Boost Software License, Version 1.0.
  4. // https://www.boost.org/LICENSE_1_0.txt
  5. // https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
  6. // Version: 4.0.2019.08.13
  7. #pragma once
  8. #include <Mathematics/Math.h>
  9. // Minimax polynomial approximations to sqrt(x). The polynomial p(x) of
  10. // degree D minimizes the quantity maximum{|sqrt(x) - p(x)| : x in [1,2]}
  11. // over all polynomials of degree D.
  12. namespace WwiseGTE
  13. {
  14. template <typename Real>
  15. class SqrtEstimate
  16. {
  17. public:
  18. // The input constraint is x in [1,2]. For example,
  19. // float x; // in [1,2]
  20. // float result = SqrtEstimate<float>::Degree<3>(x);
  21. template <int D>
  22. inline static Real Degree(Real x)
  23. {
  24. Real t = x - (Real)1; // t in [0,1]
  25. return Evaluate(degree<D>(), t);
  26. }
  27. // The input constraint is x >= 0. Range reduction is used to
  28. // generate a value y in [0,1], call Degree(y), and combine the
  29. // output with the proper exponent to obtain the approximation.
  30. // For example,
  31. // float x; // x >= 0
  32. // float result = SqrtEstimate<float>::DegreeRR<3>(x);
  33. template <int D>
  34. inline static Real DegreeRR(Real x)
  35. {
  36. Real adj, y;
  37. int p;
  38. Reduce(x, adj, y, p);
  39. Real poly = Degree<D>(y);
  40. Real result = Combine(adj, poly, p);
  41. return result;
  42. }
  43. private:
  44. // Metaprogramming and private implementation to allow specialization
  45. // of a template member function.
  46. template <int D> struct degree {};
  47. inline static Real Evaluate(degree<1>, Real t)
  48. {
  49. Real poly;
  50. poly = (Real)GTE_C_SQRT_DEG1_C1;
  51. poly = (Real)GTE_C_SQRT_DEG1_C0 + poly * t;
  52. return poly;
  53. }
  54. inline static Real Evaluate(degree<2>, Real t)
  55. {
  56. Real poly;
  57. poly = (Real)GTE_C_SQRT_DEG2_C2;
  58. poly = (Real)GTE_C_SQRT_DEG2_C1 + poly * t;
  59. poly = (Real)GTE_C_SQRT_DEG2_C0 + poly * t;
  60. return poly;
  61. }
  62. inline static Real Evaluate(degree<3>, Real t)
  63. {
  64. Real poly;
  65. poly = (Real)GTE_C_SQRT_DEG3_C3;
  66. poly = (Real)GTE_C_SQRT_DEG3_C2 + poly * t;
  67. poly = (Real)GTE_C_SQRT_DEG3_C1 + poly * t;
  68. poly = (Real)GTE_C_SQRT_DEG3_C0 + poly * t;
  69. return poly;
  70. }
  71. inline static Real Evaluate(degree<4>, Real t)
  72. {
  73. Real poly;
  74. poly = (Real)GTE_C_SQRT_DEG4_C4;
  75. poly = (Real)GTE_C_SQRT_DEG4_C3 + poly * t;
  76. poly = (Real)GTE_C_SQRT_DEG4_C2 + poly * t;
  77. poly = (Real)GTE_C_SQRT_DEG4_C1 + poly * t;
  78. poly = (Real)GTE_C_SQRT_DEG4_C0 + poly * t;
  79. return poly;
  80. }
  81. inline static Real Evaluate(degree<5>, Real t)
  82. {
  83. Real poly;
  84. poly = (Real)GTE_C_SQRT_DEG5_C5;
  85. poly = (Real)GTE_C_SQRT_DEG5_C4 + poly * t;
  86. poly = (Real)GTE_C_SQRT_DEG5_C3 + poly * t;
  87. poly = (Real)GTE_C_SQRT_DEG5_C2 + poly * t;
  88. poly = (Real)GTE_C_SQRT_DEG5_C1 + poly * t;
  89. poly = (Real)GTE_C_SQRT_DEG5_C0 + poly * t;
  90. return poly;
  91. }
  92. inline static Real Evaluate(degree<6>, Real t)
  93. {
  94. Real poly;
  95. poly = (Real)GTE_C_SQRT_DEG6_C6;
  96. poly = (Real)GTE_C_SQRT_DEG6_C5 + poly * t;
  97. poly = (Real)GTE_C_SQRT_DEG6_C4 + poly * t;
  98. poly = (Real)GTE_C_SQRT_DEG6_C3 + poly * t;
  99. poly = (Real)GTE_C_SQRT_DEG6_C2 + poly * t;
  100. poly = (Real)GTE_C_SQRT_DEG6_C1 + poly * t;
  101. poly = (Real)GTE_C_SQRT_DEG6_C0 + poly * t;
  102. return poly;
  103. }
  104. inline static Real Evaluate(degree<7>, Real t)
  105. {
  106. Real poly;
  107. poly = (Real)GTE_C_SQRT_DEG7_C7;
  108. poly = (Real)GTE_C_SQRT_DEG7_C6 + poly * t;
  109. poly = (Real)GTE_C_SQRT_DEG7_C5 + poly * t;
  110. poly = (Real)GTE_C_SQRT_DEG7_C4 + poly * t;
  111. poly = (Real)GTE_C_SQRT_DEG7_C3 + poly * t;
  112. poly = (Real)GTE_C_SQRT_DEG7_C2 + poly * t;
  113. poly = (Real)GTE_C_SQRT_DEG7_C1 + poly * t;
  114. poly = (Real)GTE_C_SQRT_DEG7_C0 + poly * t;
  115. return poly;
  116. }
  117. inline static Real Evaluate(degree<8>, Real t)
  118. {
  119. Real poly;
  120. poly = (Real)GTE_C_SQRT_DEG8_C8;
  121. poly = (Real)GTE_C_SQRT_DEG8_C7 + poly * t;
  122. poly = (Real)GTE_C_SQRT_DEG8_C6 + poly * t;
  123. poly = (Real)GTE_C_SQRT_DEG8_C5 + poly * t;
  124. poly = (Real)GTE_C_SQRT_DEG8_C4 + poly * t;
  125. poly = (Real)GTE_C_SQRT_DEG8_C3 + poly * t;
  126. poly = (Real)GTE_C_SQRT_DEG8_C2 + poly * t;
  127. poly = (Real)GTE_C_SQRT_DEG8_C1 + poly * t;
  128. poly = (Real)GTE_C_SQRT_DEG8_C0 + poly * t;
  129. return poly;
  130. }
  131. // Support for range reduction.
  132. inline static void Reduce(Real x, Real& adj, Real& y, int& p)
  133. {
  134. y = std::frexp(x, &p); // y in [1/2,1)
  135. y = (Real)2 * y; // y in [1,2)
  136. --p;
  137. adj = (1 & p) * (Real)GTE_C_SQRT_2 + (1 & ~p) * (Real)1;
  138. p >>= 1;
  139. }
  140. inline static Real Combine(Real adj, Real y, int p)
  141. {
  142. return adj * std::ldexp(y, p);
  143. }
  144. };
  145. }