fast_sin_cos

sin,cosの高速化

元のURLがみれなくなっていた.webarchiveやstackoverflowに残っていたので写経しておく.

$\sin$関数は $$ \sin0 = 0 \\ \sin\frac{\pi}{2} = 1 \\ \sin\pi = 0 $$ となる.$\sin$関数の$0 \thicksim \pi$について着目する.この範囲を二次関数に近似してみる.
$$ y = a+bx+cx^2 $$ について $$ a+b\times0+c\times0^2=0\\ a+b\times\frac{\pi}{2}+c\times(\frac{\pi}{2})^2=1\\ a+b\times\pi+c\times\pi^2=0 $$ を満たす$a, b, c$は $$ a = 0\\ b = \frac{4}{\pi}\\ c = -\frac{4}{\pi^2} $$ となるので $$ \sin\theta \thickapprox \frac{4}{\pi}\theta-\frac{4}{\pi^2}{\theta}^2 $$ $\sin$は奇関数なので,入力が負の場合は符号を反転すればよい.

if(x > 0)
{
    y = 4/pi x - 4/pi^2 x^2;
}
else
{
    y = 4/pi x + 4/pi^2 x^2;
}

$\cos$は三角関数の周期性を利用する.$\cos\theta= \sin(\theta+\frac{\pi}{2})$として計算すればよい.$\frac{\pi}{2}$を加えて$\pi$以上になるのであれば$2\pi$を減算して範囲を$-\pi \thicksim \pi$に正規化して計算する.

x += pi/2;
 
if(x > pi)   // Original x > pi/2
{
    x -= 2 * pi;   // Wrap: cos(x) = cos(x - 2 pi)
}
 
y = sine(x);

放物線近似による絶対最大誤差は0.056なので5.6%,用途によっては精度が十分かもしれないし足りないかもしれない.
グラフを書いてみると,$0, \frac{pi}{2}, \pi$を除いて近似が実際の$\sin$よりも大きい方をとっている.つまり,スケーリングすれば形状を近づけることができそう.2つの加重平均を利用してさらに近似する. $$ Q(\frac{4}{\pi}\theta-\frac{4}{\pi^2}{\theta}^2)+P(\frac{4}{\pi}\theta-\frac{4}{\pi^2}{\theta}^2) $$ ここで, $$ Q+P=1 $$ として,$Q$(と$P$)の値を設定することで,誤差を最小化できる.
絶対誤差を最小化 $$ Q = 0.775\\ P = 0.225 $$ 相対誤差を最小化 $$ Q = 0.782\\ P = 0.218 $$ ここでは絶対誤差の組み合わせを採用する.加重平均を加えた近似の最大誤差は0.001, 0.1%となり放物線のみの近似より50倍ほど精度が上がった.

最終的な実装

float sine(float x)
{
    const float B = 4/pi;
    const float C = -4/(pi*pi);
 
    float y = B * x + C * x * abs(x);
 
    #ifdef EXTRA_PRECISION
    //  const float Q = 0.775;
        const float P = 0.225;
 
        y = P * (y * abs(y) - y) + y;   // Q * y + P * y * abs(y)
    #endif
}

namespace nick_sincos {
 
template <typename T>
inline T sin(T theta) {
  static_assert(std::is_same<T, float>::value || std::is_same<T, double>::value,
                "T must be float or double");
 
  T sin;
  constexpr T pi = 3.14159265358979323846264338327950288;
  constexpr T two_pi = T(2.0) * pi;
  constexpr T B = T(4.0) / pi;
  constexpr T C = T(4.0) / (pi * pi);
  constexpr T zero = T(0.0);
  constexpr T P = T(0.225);
 
  while (theta < -pi) {
    theta += two_pi;
  }
  while (theta > pi) {
    theta -= two_pi;
  }
 
  // sin
  if (theta < zero) {
    sin = theta * B + C * theta * theta;
    if (sin < zero) {
      sin = P * (sin * -sin - sin) + sin;
    } else {
      sin = P * (sin * sin - sin) + sin;
    }
  } else {
    sin = theta * B - C * theta * theta;
    if (sin < zero) {
      sin = P * (sin * -sin - sin) + sin;
    } else {
      sin = P * (sin * sin - sin) + sin;
    }
  }
  return sin;
}
 
template <typename T>
inline T cos(T theta) {
  static_assert(std::is_same<T, float>::value || std::is_same<T, double>::value,
                "T must be float or double");
 
  T cos;
  constexpr T pi = 3.14159265358979323846264338327950288;
  constexpr T half_pi = pi / T(2.0);
  constexpr T two_pi = T(2.0) * pi;
  constexpr T B = T(4.0) / pi;
  constexpr T C = T(4.0) / (pi * pi);
  constexpr T zero = T(0.0);
  constexpr T P = T(0.225);
 
  while (theta < -pi) {
    theta += two_pi;
  }
  while (theta > pi) {
    theta -= two_pi;
  }
  theta += half_pi;
  if (theta > pi) {
    theta -= two_pi;
  }
  // cos
  if (theta < zero) {
    cos = theta * B + C * theta * theta;
    if (cos < zero) {
      cos = P * (cos * -cos - cos) + cos;
    } else {
      cos = P * (cos * cos - cos) + cos;
    }
  } else {
    cos = theta * B - C * theta * theta;
    if (cos < zero) {
      cos = P * (cos * -cos - cos) + cos;
    } else {
      cos = P * (cos * cos - cos) + cos;
    }
  }
  return cos;
}
 
template <typename T>
inline std::tuple<T, T> sincos(T theta) {
  static_assert(std::is_same<T, float>::value || std::is_same<T, double>::value,
                "T must be float or double");
 
  T sin_out = sin(theta);
  T cos_tou = cos(theta);
  return {sin_out, cos_tou};
}
}  // namespace nick_sincos

参考文献

  • fast_sin_cos.txt
  • 最終更新: 2024/06/24
  • by yuqlid