Go 言語の math.Sin は計算結果がおかしいことがある

# math.Sin

Go で三角関数を利用していたら、 `math.Sin(x)` の結果が 1 を超えることがあった。

調べて見ると、計算結果がおかしくなるのは引数が `1e15` などの巨大な数の場合のみみたいなので困ることは少ないと思う。

# 他言語との比較

各言語での計算結果を下表にまとめた。
`"%.3g"` で文字列化した結果が C++ と異なる場合に太字にしてある。

|θ|c++|ruby|python3|node|groovy|go|
|:--|:--|:--|:--|:--|:--|:--|
|2**1|0.909|0.909|0.909|0.909|0.909|0.909|
|2**40|-0.406|-0.406|-0.406|-0.406|-0.406|-0.406|
|2**41|0.742|0.742|0.742|0.742|0.742|0.742|
|2**42|0.995|0.995|0.995|0.995|0.995|0.995|
|2**43|-0.199|-0.199|-0.199|-0.199|-0.199|-0.199|
|2**44|0.39|0.39|0.39|0.39|0.39|0.39|
|2**45|0.718|0.718|0.718|0.718|0.718|**0.719**|
|2**46|0.999|0.999|0.999|0.999|0.999|0.999|
|2**47|-0.0649|-0.0649|-0.0649|-0.0649|-0.0649|**-0.0673**|
|2**48|0.129|0.129|0.129|0.129|0.129|**0.134**|
|2**49|0.257|0.257|0.257|0.257|0.257|**0.266**|
|2**50|0.496|0.496|0.496|0.496|0.496|**0.513**|
|2**51|0.862|0.862|0.862|0.862|0.862|**0.912**|
|2**52|0.874|0.874|0.874|0.874|0.874|**0.792**|
|2**53|-0.849|-0.849|-0.849|-0.849|-0.849|**-0.967**|
|2**54|0.897|0.897|0.897|0.897|0.897|**0.492**|
|2**55|-0.792|-0.792|-0.792|-0.792|-0.792|**-0.816**|
|2**56|0.967|0.967|0.967|0.967|0.967|**1.09e+03**|
|2**57|-0.493|-0.493|-0.493|-0.493|-0.493|**1.79e+07**|
|2**58|0.858|0.858|0.858|0.858|0.858|**1.86e+11**|
|2**59|0.882|0.882|0.882|0.882|0.882|**1.63e+15**|
|2**60|-0.831|-0.831|-0.831|-0.831|-0.831|**1.36e+19**|
|2**61|0.925|0.925|0.925|0.925|0.925|**1.12e+23**|
|2**62|-0.703|-0.703|-0.703|-0.703|-0.703|**9.16e+26**|
|2**63|1|1|1|1|1|**1.04e+240**|
|2**64|0.0236|0.0236|0.0236|0.0236|0.0236|**3.38e+242**|

`2**45` でずれ始め、`2**56` で完全におかしくなっている。

# src/math/sin.go

[src/math/sin.go](https://golang.org/src/math/sin.go) を見ると

> Results may be meaningless for x > 2**49 = 5.6e14.

と書いてある。
meaningless な値を返すぐらいなら、仕様に有効な値の範囲を明記して、それ以外の入力に対しては非数を返してほしいと思う。

# sin 以外の三角関数について

調べてみたら、sin, cos, tan は go は他の言語と違う結果になる感じ。
atan, atan2 は大丈夫そう。
asin, acos は、大きな入力を受け付けないので調べていない。

結果は以下の通り：


## t = 1e+12

|処理系|sin(t)|cos(t)|tan(t)|atan(t)|atan2(t,t/2)|
|:--|:--:|:--:|:--:|:--:|:--:|
|c++|-0.611|0.791|-0.772|1.57|1.11|
|node|-0.611|0.791|-0.772|1.57|1.11|
|python3|-0.611|0.791|-0.772|1.57|1.11|
|groovy|-0.611|0.791|-0.772|1.57|1.11|
|go|-0.611|0.791|-0.772|1.57|1.11|

## t = 1e+14

|処理系|sin(t)|cos(t)|tan(t)|atan(t)|atan2(t,t/2)|
|:--|:--:|:--:|:--:|:--:|:--:|
|c++|-0.209|-0.978|0.214|1.57|1.11|
|node|-0.209|-0.978|0.214|1.57|1.11|
|python3|-0.209|-0.978|0.214|1.57|1.11|
|groovy|-0.209|-0.978|0.214|1.57|1.11|
|go|**-0.217**|**-0.976**|**0.222**|1.57|1.11|

## t = 1e+16

|処理系|sin(t)|cos(t)|tan(t)|atan(t)|atan2(t,t/2)|
|:--|:--:|:--:|:--:|:--:|:--:|
|c++|0.78|-0.626|-1.25|1.57|1.11|
|node|0.78|-0.626|-1.25|1.57|1.11|
|python3|0.78|-0.626|-1.25|1.57|1.11|
|groovy|0.78|-0.626|-1.25|1.57|1.11|
|go|**0.969**|**0.248**|**3.91**|1.57|1.11|

## t = 1e+308

|処理系|sin(t)|cos(t)|tan(t)|atan(t)|atan2(t,t/2)|
|:--|:--:|:--:|:--:|:--:|:--:|
|c++|0.453|-0.891|-0.509|1.57|1.11|
|node|0.453|-0.891|-0.509|1.57|1.11|
|python3|0.453|-0.891|-0.509|1.57|1.11|
|groovy|0.453|-0.891|-0.509|1.57|1.11|
|go|**Inf**|**-Inf**|**NaN**|1.57|1.11|

まあ、1e16 を超える数を sin, cos, tan に入れても意味が無いのはわかるけど、sin, cos なら絶対値が 1 以下の値がかえってきてほしいと思う。あと、$(sin(t))^{2} + (cos(t))^{2}$ は (多少の誤差は許容するけど) 1 になってほしい。

# 実行速度

go の sin, cos は独自実装のようなので、実行速度がちょっと気になった。

で。
測ってみた。

## C++ のソースコード

```c++:c++
#include <iostream>
#include <iomanip>
#include <cmath>

double run(int size){
  double sum = 0;
  for( int y=0 ; y<size ; ++y ){
    for( int x=0 ; x<size ; ++x ){
      sum += std::sin(x+y+sum);
      sum += std::cos(x-y+sum);
    }
  }
  return sum;
}

int main(){
  int size;
  std::cin >> size;
  std::cout << std::setprecision(20) << run(size) << std::endl;
  return 0;
}
```

## go のソースコード

```go:go
package main

import (
	"bufio"
	"fmt"
	"math"
	"os"
	"strconv"
)

func run(size int) float64 {
	sum := 0.0
	for y := 0; y < size; y++ {
		for x := 0; x < size; x++ {
			sum += math.Sin(float64(x+y) + sum)
			sum += math.Cos(float64(x-y) + sum)
		}
	}
	return sum
}

func main() {
	stdin := bufio.NewScanner(os.Stdin)
	stdin.Scan()
	text := stdin.Text()
	size, _ := strconv.ParseInt(text, 10, 32)
	fmt.Printf("%.20f\n", run(int(size)))
}
```

## 測定

```
echo 5000|time a.out
```

みたいな感じ。user と real がほぼ一致し、sys がゼロだったので、user の値のみを表にまとめた。

## 測定結果


|コンパイラ|最適化オプション|時間|
|:--|:--|:--|
|clang++|-O0|2.24|
|clang++|-O1|2.13|
|clang++|-O2|2.14|
|clang++|-O3|2.14|
|clang++|-Ofast|2.14|
|g++-8|-O0|2.21|
|g++-8|-O1|2.19|
|g++-8|-O2|2.15|
|g++-8|-O3|2.14|
|g++-8|-Ofast|2.21|
|go|-gcflags '-N -l'|1.70|
|go||1.60|

※ C++ のコンパイルオプションは他に「`-mtune=native -march=native`」を指定している

go は独自実装によって高速化を実現したのかな。
この実験だけではっきりしたことは言えないけれど。