SinEstimate.h 4.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136
  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 sin(x). The polynomial p(x) of
  10. // degree D has only odd-power terms, is required to have linear term x,
  11. // and p(pi/2) = sin(pi/2) = 1. It minimizes the quantity
  12. // maximum{|sin(x) - p(x)| : x in [-pi/2,pi/2]} over all polynomials of
  13. // degree D subject to the constraints mentioned.
  14. namespace WwiseGTE
  15. {
  16. template <typename Real>
  17. class SinEstimate
  18. {
  19. public:
  20. // The input constraint is x in [-pi/2,pi/2]. For example,
  21. // float x; // in [-pi/2,pi/2]
  22. // float result = SinEstimate<float>::Degree<3>(x);
  23. template <int D>
  24. inline static Real Degree(Real x)
  25. {
  26. return Evaluate(degree<D>(), x);
  27. }
  28. // The input x can be any real number. Range reduction is used to
  29. // generate a value y in [-pi/2,pi/2] for which sin(y) = sin(x).
  30. // For example,
  31. // float x; // x any real number
  32. // float result = SinEstimate<float>::DegreeRR<3>(x);
  33. template <int D>
  34. inline static Real DegreeRR(Real x)
  35. {
  36. return Degree<D>(Reduce(x));
  37. }
  38. private:
  39. // Metaprogramming and private implementation to allow specialization
  40. // of a template member function.
  41. template <int D> struct degree {};
  42. inline static Real Evaluate(degree<3>, Real x)
  43. {
  44. Real xsqr = x * x;
  45. Real poly;
  46. poly = (Real)GTE_C_SIN_DEG3_C1;
  47. poly = (Real)GTE_C_SIN_DEG3_C0 + poly * xsqr;
  48. poly = poly * x;
  49. return poly;
  50. }
  51. inline static Real Evaluate(degree<5>, Real x)
  52. {
  53. Real xsqr = x * x;
  54. Real poly;
  55. poly = (Real)GTE_C_SIN_DEG5_C2;
  56. poly = (Real)GTE_C_SIN_DEG5_C1 + poly * xsqr;
  57. poly = (Real)GTE_C_SIN_DEG5_C0 + poly * xsqr;
  58. poly = poly * x;
  59. return poly;
  60. }
  61. inline static Real Evaluate(degree<7>, Real x)
  62. {
  63. Real xsqr = x * x;
  64. Real poly;
  65. poly = (Real)GTE_C_SIN_DEG7_C3;
  66. poly = (Real)GTE_C_SIN_DEG7_C2 + poly * xsqr;
  67. poly = (Real)GTE_C_SIN_DEG7_C1 + poly * xsqr;
  68. poly = (Real)GTE_C_SIN_DEG7_C0 + poly * xsqr;
  69. poly = poly * x;
  70. return poly;
  71. }
  72. inline static Real Evaluate(degree<9>, Real x)
  73. {
  74. Real xsqr = x * x;
  75. Real poly;
  76. poly = (Real)GTE_C_SIN_DEG9_C4;
  77. poly = (Real)GTE_C_SIN_DEG9_C3 + poly * xsqr;
  78. poly = (Real)GTE_C_SIN_DEG9_C2 + poly * xsqr;
  79. poly = (Real)GTE_C_SIN_DEG9_C1 + poly * xsqr;
  80. poly = (Real)GTE_C_SIN_DEG9_C0 + poly * xsqr;
  81. poly = poly * x;
  82. return poly;
  83. }
  84. inline static Real Evaluate(degree<11>, Real x)
  85. {
  86. Real xsqr = x * x;
  87. Real poly;
  88. poly = (Real)GTE_C_SIN_DEG11_C5;
  89. poly = (Real)GTE_C_SIN_DEG11_C4 + poly * xsqr;
  90. poly = (Real)GTE_C_SIN_DEG11_C3 + poly * xsqr;
  91. poly = (Real)GTE_C_SIN_DEG11_C2 + poly * xsqr;
  92. poly = (Real)GTE_C_SIN_DEG11_C1 + poly * xsqr;
  93. poly = (Real)GTE_C_SIN_DEG11_C0 + poly * xsqr;
  94. poly = poly * x;
  95. return poly;
  96. }
  97. // Support for range reduction.
  98. inline static Real Reduce(Real x)
  99. {
  100. // Map x to y in [-pi,pi], x = 2*pi*quotient + remainder.
  101. Real quotient = (Real)GTE_C_INV_TWO_PI * x;
  102. if (x >= (Real)0)
  103. {
  104. quotient = (Real)((int)(quotient + (Real)0.5));
  105. }
  106. else
  107. {
  108. quotient = (Real)((int)(quotient - (Real)0.5));
  109. }
  110. Real y = x - (Real)GTE_C_TWO_PI * quotient;
  111. // Map y to [-pi/2,pi/2] with sin(y) = sin(x).
  112. if (y > (Real)GTE_C_HALF_PI)
  113. {
  114. y = (Real)GTE_C_PI - y;
  115. }
  116. else if (y < (Real)-GTE_C_HALF_PI)
  117. {
  118. y = (Real)-GTE_C_PI - y;
  119. }
  120. return y;
  121. }
  122. };
  123. }