ラジアンの引数

# 目立つ誤差

整数 $n$ に対する $sin\left(n \pi \div 2 \right)$ を表示するプログラムと実行結果

```c:sample.c
#include <stdio.h>
#include <math.h>

void print_fp64(double x)
{
  union {
    double f;
    unsigned long long i;
  } t;

  unsigned long long u;
  unsigned long long m;
  int s, e;

  t.f = x;
  u = t.i;

  s = (u >> 63);
  e = (u >> 52) & 0x7ff;
  m = (u & ((1ULL << 52) - 1));
  printf("(%c0x%d.%015llx * 2^%-5d)",
         (s ? '-' : '+'), (e ? 1 : 0), m, (e - 1023));
}

int main(int argc, char **argv)
{
  int n;
  for (n = 0; n < 10; n++)
    {
      double r = M_PI/2 * n;
      double s = sin(r);

      printf("%2d", n);
      printf(":  r = %+e  ", r); print_fp64(r);
      printf(":  s = %+e  ", s); print_fp64(s);
      printf("\n");
    }
  return 0;
}
```

環境： Apple clang version 11.0.0 (clang-1100.0.33.17), Target: x86_64-apple-darwin19.3.0

```console:
$ ./sample
 0:  r = +0.000000e+00  (+0x0.000000000000000 * 2^-1023):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
 1:  r = +1.570796e+00  (+0x1.00921fb54442d18 * 2^0    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
 2:  r = +3.141593e+00  (+0x1.00921fb54442d18 * 2^1    ):  s = +1.224647e-16  (+0x1.001a62633145c07 * 2^-53  )
 3:  r = +4.712389e+00  (+0x1.002d97c7f3321d2 * 2^2    ):  s = -1.000000e+00  (-0x1.000000000000000 * 2^0    )
 4:  r = +6.283185e+00  (+0x1.00921fb54442d18 * 2^2    ):  s = -2.449294e-16  (-0x1.001a62633145c07 * 2^-52  )
 5:  r = +7.853982e+00  (+0x1.00f6a7a2955385e * 2^2    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
 6:  r = +9.424778e+00  (+0x1.002d97c7f3321d2 * 2^3    ):  s = +3.673940e-16  (+0x1.00a79394c9e8a0a * 2^-52  )
 7:  r = +1.099557e+01  (+0x1.005fdbbe9bba775 * 2^3    ):  s = -1.000000e+00  (-0x1.000000000000000 * 2^0    )
 8:  r = +1.256637e+01  (+0x1.00921fb54442d18 * 2^3    ):  s = -4.898587e-16  (-0x1.001a62633145c07 * 2^-51  )
 9:  r = +1.413717e+01  (+0x1.00c463abeccb2bb * 2^3    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
```

となりました。括弧の中は、[浮動小数点データ形式](https://qiita.com/ikiuo/items/6c02524a35f676ec753b#%E6%B5%AE%E5%8B%95%E5%B0%8F%E6%95%B0%E7%82%B9%E3%83%87%E3%83%BC%E3%82%BF%E5%BD%A2%E5%BC%8F)の数値です。

s は 0,+1,0,-1,0,... の繰り返しのはずですが演算誤差がでています。ここでの誤差 $\left(10^{-16} 程度\right)$ は十分小さいと思いますが、もともと切りの良い数字は **誤差が目立ちます** 。

# 近似値になる引数

$\pi$ は一般的な CPU の演算型では正確に表現できません。(参考：[π÷2 の2進数](https://qiita.com/ikiuo/items/6c02524a35f676ec753b#%CF%802-%E3%81%AE2%E9%80%B2%E6%95%B0))

使用したい $\pi \div 2$ は、double 型の
　　1.5707963267948965579989817342720925807952880859375
と
　　1.5707963267948967800435866593034006655216217041015625
の間にあります。

引数は前者の近似値を整数倍しているので、意図したものと誤差があります。

引数の時点で誤差があるので、結果に誤差があっても仕方がない気がしますが、引数 $r$ の仮数部を見ると 1,2,4,8 倍のときは同じだから、引数 $r$ を受け取った側での改善余地がありそうです。単純に

```c:
#define sin(r) sin(fmod((r), M_PI*2))
```

で正規化してみると

```console:
 0:  r = +0.000000e+00  (+0x0.000000000000000 * 2^-1023):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
 1:  r = +1.570796e+00  (+0x1.00921fb54442d18 * 2^0    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
 2:  r = +3.141593e+00  (+0x1.00921fb54442d18 * 2^1    ):  s = +1.224647e-16  (+0x1.001a62633145c07 * 2^-53  )
 3:  r = +4.712389e+00  (+0x1.002d97c7f3321d2 * 2^2    ):  s = -1.000000e+00  (-0x1.000000000000000 * 2^0    )
 4:  r = +6.283185e+00  (+0x1.00921fb54442d18 * 2^2    ):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
 5:  r = +7.853982e+00  (+0x1.00f6a7a2955385e * 2^2    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
 6:  r = +9.424778e+00  (+0x1.002d97c7f3321d2 * 2^3    ):  s = +1.224647e-16  (+0x1.001a62633145c07 * 2^-53  )
 7:  r = +1.099557e+01  (+0x1.005fdbbe9bba775 * 2^3    ):  s = -1.000000e+00  (-0x1.000000000000000 * 2^0    )
 8:  r = +1.256637e+01  (+0x1.00921fb54442d18 * 2^3    ):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
 9:  r = +1.413717e+01  (+0x1.00c463abeccb2bb * 2^3    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
...
22:  r = +3.455752e+01  (+0x1.001475cc9eedf00 * 2^5    ):  s = +3.675178e-15  (+0x1.0008d313198a2e0 * 2^-48  )
...
44:  r = +6.911504e+01  (+0x1.001475cc9eedf00 * 2^6    ):  s = -7.350357e-15  (-0x1.0008d313198a2e0 * 2^-47  )
...
```

切りの良い数字でないのは $\left(4k+2\right)$ 倍のときになり、22倍未満までは 0,1,2,3 倍の繰り返しに改善しました。さらに、象限まで分解したものを挟んで

```c:
double xsin(double r)
{
  double pi = M_PI;
  double pi2 = pi * 2;
  double hpi = pi / 2;
  double s = 1;
  int n;

  if (r < 0)
    {
      s = -s;
      r = -r;
    }
  n = r / pi2;
  r -= pi2 * n;
  if (r >= pi2)
    r -= pi2;
  if (r >= pi)
    {
      s = -s;
      r -= pi;
    }
  if (r >= hpi)
    r = pi - r;
  return (s * sin(r));
}
#define sin(r) xsin(r)
```

とすると

```console
 0:  r = +0.000000e+00  (+0x0.000000000000000 * 2^-1023):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
 1:  r = +1.570796e+00  (+0x1.00921fb54442d18 * 2^0    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
 2:  r = +3.141593e+00  (+0x1.00921fb54442d18 * 2^1    ):  s = -0.000000e+00  (-0x0.000000000000000 * 2^-1023)
 3:  r = +4.712389e+00  (+0x1.002d97c7f3321d2 * 2^2    ):  s = -1.000000e+00  (-0x1.000000000000000 * 2^0    )
 4:  r = +6.283185e+00  (+0x1.00921fb54442d18 * 2^2    ):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
 5:  r = +7.853982e+00  (+0x1.00f6a7a2955385e * 2^2    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
 6:  r = +9.424778e+00  (+0x1.002d97c7f3321d2 * 2^3    ):  s = -0.000000e+00  (-0x0.000000000000000 * 2^-1023)
 7:  r = +1.099557e+01  (+0x1.005fdbbe9bba775 * 2^3    ):  s = -1.000000e+00  (-0x1.000000000000000 * 2^0    )
 8:  r = +1.256637e+01  (+0x1.00921fb54442d18 * 2^3    ):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
 9:  r = +1.413717e+01  (+0x1.00c463abeccb2bb * 2^3    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
10:  r = +1.570796e+01  (+0x1.00f6a7a2955385e * 2^3    ):  s = -0.000000e+00  (-0x0.000000000000000 * 2^-1023)
11:  r = +1.727876e+01  (+0x1.001475cc9eedf00 * 2^4    ):  s = -1.000000e+00  (-0x1.000000000000000 * 2^0    )
12:  r = +1.884956e+01  (+0x1.002d97c7f3321d2 * 2^4    ):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
13:  r = +2.042035e+01  (+0x1.0046b9c347764a4 * 2^4    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
14:  r = +2.199115e+01  (+0x1.005fdbbe9bba775 * 2^4    ):  s = -0.000000e+00  (-0x0.000000000000000 * 2^-1023)
15:  r = +2.356194e+01  (+0x1.0078fdb9effea46 * 2^4    ):  s = -1.000000e+00  (-0x1.000000000000000 * 2^0    )
16:  r = +2.513274e+01  (+0x1.00921fb54442d18 * 2^4    ):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
17:  r = +2.670354e+01  (+0x1.00ab41b09886fea * 2^4    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
18:  r = +2.827433e+01  (+0x1.00c463abeccb2bb * 2^4    ):  s = -0.000000e+00  (-0x0.000000000000000 * 2^-1023)
19:  r = +2.984513e+01  (+0x1.00dd85a7410f58c * 2^4    ):  s = -1.000000e+00  (-0x1.000000000000000 * 2^0    )
20:  r = +3.141593e+01  (+0x1.00f6a7a2955385e * 2^4    ):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
21:  r = +3.298672e+01  (+0x1.0007e4cef4cbd98 * 2^5    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
22:  r = +3.455752e+01  (+0x1.001475cc9eedf00 * 2^5    ):  s = +3.552714e-15  (+0x1.000000000000000 * 2^-48  )
23:  r = +3.612832e+01  (+0x1.002106ca4910069 * 2^5    ):  s = -1.000000e+00  (-0x1.000000000000000 * 2^0    )
24:  r = +3.769911e+01  (+0x1.002d97c7f3321d2 * 2^5    ):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
25:  r = +3.926991e+01  (+0x1.003a28c59d5433b * 2^5    ):  s = +1.000000e+00  (+0x1.000000000000000 * 2^0    )
26:  r = +4.084070e+01  (+0x1.0046b9c347764a4 * 2^5    ):  s = -3.552714e-15  (-0x1.000000000000000 * 2^-48  )
27:  r = +4.241150e+01  (+0x1.00534ac0f19860c * 2^5    ):  s = -1.000000e+00  (-0x1.000000000000000 * 2^0    )
28:  r = +4.398230e+01  (+0x1.005fdbbe9bba775 * 2^5    ):  s = +0.000000e+00  (+0x0.000000000000000 * 2^-1023)
...
```

22倍未満までは、0,+1,0,-1,0,... になりました。

引数を受け取る側で、ある程度の改善が可能ですが、引数が大きくなるに従って、整数部の情報が増えて小数部の情報が減るので、正規化の誤差も次第に大きくなっていく状況は避けられません。


# 引数はラジアンではないものを使う

浮動小数点演算が使えなかった 16bit 機では、整数の
　　$0 \le x < 4096$
範囲を、一回転分として
　　$0 \le x < 2 \pi$
の代わりにしていました。範囲の正規化は
　　(x & 4095)
の１ステップで完了です。sin,cos 関数などはテーブル参照です。

角分解能は $0.1^\circ$ 弱ですが、これでも作っていたソフトでは十分でした。

浮動小数点形式でも同様に、引数を $\left(2\pi→1.0\right)$ として正規化する方法が考えられますが、「度」の単位でも整数部側に $\left(2\pi→360\right)$ がくるので、馴染みのある「度」で扱うのがいいのかもしれません。

```c:
double deg_sin(double x)
{
  double s = 1;

  if (x < 0)
    {
      s = -s;
      x = -x;
    }
  x = fmod(x, 360.0);
  if (x >= 180.0)
    {
      s = -s;
      x -= 180.0;
    }
  if (x >= 90.0)
    x = 180.0 - x;
  return (s * sin(x * M_PI / 180));
}
```

---