[インデックス 10156] ファイルの概要

このコミットは、Go言語のmathパッケージにおけるSin（サイン）およびCos（コサイン）関数の精度向上を目的としています。特に、既存の実装が抱えていた精度に関する問題を解決し、より広範な入力値に対して正確な結果を返すように改善されています。

コミット

commit 739c442e42312ee843f4736d10f73f9c7d292226
Author: Charles L. Dorian <cldorian@gmail.com>
Date:   Mon Oct 31 14:26:05 2011 -0400

    math: Improved accuracy for Sin and Cos.
    
    Fixes #1564.
    
    R=rsc, dchest
    CC=golang-dev
    https://golang.org/cl/5320056

GitHub上でのコミットページへのリンク

https://github.com/golang/go/commit/739c442e42312ee843f4736d10f73f9c7d292226

元コミット内容

このコミットの元のメッセージは以下の通りです。

math: Improved accuracy for Sin and Cos.

Fixes #1564.

R=rsc, dchest
CC=golang-dev
https://golang.org/cl/5320056

これは、mathパッケージのSinおよびCos関数の精度が改善されたことを示しています。Fixes #1564は、GoのIssueトラッカーで報告されていたIssue 1564を修正したことを意味します。

変更の背景

この変更の背景には、Go言語のmathパッケージにおける三角関数（特にサインとコサイン）の計算精度に関する問題がありました。Issue #1564で報告されたように、既存の実装では特定の入力値に対して期待される精度が得られないケースが存在していました。特に、大きな数値や特定の境界値において、浮動小数点演算の性質上、誤差が蓄積しやすくなることが課題でした。

このような数値計算の精度問題は、科学技術計算、グラフィックス、物理シミュレーションなど、正確な浮動小数点演算が求められる多くのアプリケーションにおいて致命的な影響を与える可能性があります。そのため、標準ライブラリとして提供されるmathパッケージの関数は、可能な限り高い精度とロバスト性を持つことが求められます。

このコミットは、より信頼性の高い数値計算を提供するために、既存のアルゴリズムをより高精度なものに置き換えることで、この課題を解決しようとしたものです。

前提知識の解説

浮動小数点数と精度

コンピュータにおける浮動小数点数は、実数を近似的に表現するための形式です。IEEE 754規格が広く用いられており、Go言語のfloat64型もこれに準拠しています。浮動小数点数は、有限のビット数で表現されるため、全ての実数を正確に表現することはできません。このため、計算には常に「丸め誤差」が伴います。

三角関数のような超越関数を計算する場合、テイラー展開やチェビシェフ近似などの多項式近似が用いられます。これらの近似は、ある範囲内で関数を多項式で近似することで計算を効率化しますが、近似の次数や係数の選択、そして入力値の範囲（定義域）によって精度が大きく左右されます。特に、入力値が大きくなると、周期関数の性質上、引数を適切な範囲に「還元」する処理が必要になりますが、この還元処理自体にも誤差が伴い、最終的な結果の精度に影響を与えます。

範囲還元 (Range Reduction)

三角関数は周期関数であるため、任意の入力値xに対して、xを2πの倍数で割った余り（主値）を計算し、その主値に対して関数を評価することで、計算の範囲を限定することができます。このプロセスを「範囲還元」と呼びます。

例えば、sin(x)を計算する際に、xが非常に大きな値であっても、x = 2nπ + θ（ここで0 <= θ < 2π）と表現できるため、sin(x) = sin(θ)として計算できます。しかし、xが非常に大きい場合、xから2nπを引く際に、xと2nπの有効桁数が大きく異なるため、減算によって有効桁数が失われ、結果としてθの精度が著しく低下する可能性があります。これを「桁落ち」と呼びます。

高精度な範囲還元を実現するためには、πの値を非常に高い精度で保持し、多倍長演算のような手法を用いて桁落ちを防ぐ必要があります。

Cephes Math Library

Cephes Math Libraryは、Stephen L. Moshierによって開発された、C言語で書かれた数値計算ライブラリです。このライブラリは、様々な特殊関数や超越関数（三角関数、指数関数、対数関数、ガンマ関数など）の高精度な実装を提供しており、その堅牢性と精度で知られています。多くの科学技術計算ソフトウェアやライブラリで、そのアルゴリズムや実装が参考にされています。

このコミットでは、SinおよびCos関数の実装において、Cephes Math Libraryのアルゴリズムと係数が採用されています。これは、Cephesライブラリが提供する高精度な計算手法を取り入れることで、Goのmathパッケージの精度を向上させることを意図しています。

技術的詳細

このコミットの主要な変更点は、src/pkg/math/sin.goファイルにおけるSinおよびCos関数の実装が、従来のHart & Cheneyの係数を用いたものから、Cephes Math Libraryのアルゴリズムと係数を用いたものへと全面的に刷新されたことです。

従来のsinus関数（削除されたコード）

変更前のsinus関数は、xを2/Piで乗算して正規化し、Modf関数を使って整数部と小数部に分け、象限（quad）に基づいて符号と値を調整していました。そして、PとQという係数を用いた有理関数近似（多項式の比）によって最終的な値を計算していました。

func sinus(x float64, quad int) float64 {
	const (
		P0 = .1357884097877375669092680e8
		// ...
		Q0 = .8644558652922534429915149e7
		// ...
	)
	// ... 範囲還元と象限の調整 ...
	yy := y * y
	temp1 := ((((P4*yy+P3)*yy+P2)*yy+P1)*yy + P0) * y
	temp2 := ((((yy+Q3)*yy+Q2)*yy+Q1)*yy + Q0)
	return temp1 / temp2
}

この実装は、特定の範囲では機能しますが、特に大きな入力値に対する範囲還元において精度が問題となる可能性がありました。

新しいSinおよびCos関数の実装（追加されたコード）

新しい実装は、Cephes Math Libraryの設計思想に基づいています。

1. 特殊ケースのハンドリング: NaNやInfといった特殊な浮動小数点値に対する挙動が明示的に定義されています。

   switch {
   case x != x || x < -MaxFloat64 || x > MaxFloat64: // IsNaN(x) || IsInf(x, 0):
       return NaN()
   }
   

   Sin関数ではx == 0の場合の±0のハンドリングも追加されています。

2. 引数の正数化と符号の保存: 入力xが負の場合、絶対値を取り、最終的な結果の符号を決定するためにsignフラグを立てます。

   sign := false
   if x < 0 {
       x = -x
       sign = true // Sinの場合
       // Cosの場合、xが負でもCos(x) = Cos(-x)なので、符号は変わらない
   }
   

3. 高精度な範囲還元: ここが最も重要な変更点です。従来の2/Piによる単純な正規化ではなく、PI4A, PI4B, PI4Cという3つの定数に分割されたπ/4と、M4PI（4/π）を用いて、拡張精度（extended precision）でのモジュラ演算を行います。

   const (
       PI4A = 7.85398125648498535156E-1 // 0x3fe921fb40000000, Pi/4 split into three parts
       PI4B = 3.77489470793079817668E-8 // 0x3e64442d00000000,
       PI4C = 2.69515142907905952645E-15 // 0x3ce8469898cc5170,
       M4PI = 1.273239544735162542821171882678754627704620361328125 // 4/pi
   )
   
   j := int64(x * M4PI) // integer part of x/(Pi/4)
   y := float64(j)      // integer part as float
   
   // map zeros to origin
   if j&1 == 1 {
       j += 1
       y += 1
   }
   
   j &= 7 // octant modulo 2Pi radians (360 degrees)
   // reflect in x axis (Sinの場合)
   if j > 3 {
       sign = !sign
       j -= 4
   }
   // Cosの場合の象限調整
   if j > 3 {
       j -= 4
       sign = !sign
   }
   if j > 1 {
       sign = !sign
   }
   
   z := ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
   zz := z * z
   

   このz := ((x - y*PI4A) - y*PI4B) - y*PI4Cという計算は、xからy * (PI4A + PI4B + PI4C)を引く操作を、複数の浮動小数点数で段階的に行うことで、桁落ちによる精度損失を最小限に抑えるためのテクニックです。PI4A, PI4B, PI4Cはπ/4を多倍長精度で表現したものであり、これによりxをπ/4の倍数で割った余りzを非常に高い精度で計算できます。

4. 多項式近似: 還元されたzとzz（z*z）を用いて、定義域[0, π/4]におけるサインまたはコサインの多項式近似を計算します。

   // sin coefficients
   var _sin = [...]float64{
       1.58962301576546568060E-10, // ...
       // ...
   }
   // cos coefficients
   var _cos = [...]float64{
       -1.13585365213876817300E-11, // ...
       // ...
   }
   
   if j == 1 || j == 2 { // Cos(z)を計算する象限
       y = 1.0 - 0.5*zz + zz*zz*((((((_cos[0]*zz)+_cos[1])*zz+_cos[2])*zz+_cos[3])*zz+_cos[4])*zz+_cos[5])
   } else { // Sin(z)を計算する象限
       y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
   }
   

   _sinと_cosは、Cephes Math Libraryから導出された多項式近似の係数です。jの値（象限）に応じて、sin(z)またはcos(z)のどちらかの多項式を評価します。これは、sin(x)とcos(x)がπ/2の位相差を持つことを利用し、常に[0, π/4]の範囲で最も収束の良い多項式近似を使用するためです。例えば、sin(x)を計算する際にxがπ/4からπ/2の範囲にある場合、cos(π/2 - x)としてコサインの多項式で計算する方が精度が良い、といった最適化が行われています。

5. 最終的な符号の適用: 最初に保存したsignフラグに基づいて、最終結果の符号を調整します。

テストコードの変更

src/pkg/math/all_test.goでは、TestCos, TestSin, TestSincos関数において、結果の比較に使用されるアサーション関数がcloseからverycloseに変更されています。

--- a/src/pkg/math/all_test.go
+++ b/src/pkg/math/all_test.go
@@ -1709,7 +1709,7 @@ func TestCopysign(t *testing.T) {
 
 func TestCos(t *testing.T) {
 	for i := 0; i < len(vf); i++ {
-		if f := Cos(vf[i]); !close(cos[i], f) {
+		if f := Cos(vf[i]); !veryclose(cos[i], f) {
 			t.Errorf("Cos(%g) = %g, want %g", vf[i], f, cos[i])
 		}
 	}
@@ -2192,7 +2192,7 @@ func TestSignbit(t *testing.T) {
 }
 func TestSin(t *testing.T) {
 	for i := 0; i < len(vf); i++ {
-		if f := Sin(vf[i]); !close(sin[i], f) {
+		if f := Sin(vf[i]); !veryclose(sin[i], f) {
 			t.Errorf("Sin(%g) = %g, want %g", vf[i], f, sin[i])
 		}
 	}
@@ -2205,7 +2205,7 @@ func TestSincos(t *testing.T) {
 
 func TestSincos(t *testing.T) {
 	for i := 0; i < len(vf); i++ {
-		if s, c := Sincos(vf[i]); !close(sin[i], s) || !close(cos[i], c) {
+		if s, c := Sincos(vf[i]); !veryclose(sin[i], s) || !veryclose(cos[i], c) {
 			t.Errorf("Sincos(%g) = %g, %g want %g, %g", vf[i], s, c, sin[i], cos[i])
 	}
 }

これは、新しい実装がより高精度になったため、テストにおける許容誤差の基準をより厳しくしたことを示唆しています。veryclose関数は、close関数よりもさらに厳密な浮動小数点数の比較を行うものと推測されます。

コアとなるコードの変更箇所

このコミットのコアとなるコードの変更は、主にsrc/pkg/math/sin.goファイルに集中しています。

1. 著作権表示の更新: // Copyright 2009 The Go Authors. All rights reserved. から // Copyright 2011 The Go Authors. All rights reserved. へ変更。

2. sinus関数の削除: 従来のsinus関数（SinとCosの両方から呼び出されていた内部ヘルパー関数）が完全に削除されました。

3. Cephes Math Libraryに関するコメントの追加: 新しいSinとCos関数の実装がCephes Math Libraryに基づいていることを説明する詳細なコメントブロックが追加されました。これには、Cephesライブラリの概要、精度情報、著作権情報などが含まれています。

4. サイン・コサイン係数の定義: Cephesライブラリから導出された多項式近似のための係数配列_sinと_cosがグローバル変数として定義されました。

   var _sin = [...]float64{
       1.58962301576546568060E-10, // 0x3de5d8fd1fd19ccd
       -2.50507477628578072866E-8, // 0xbe5ae5e5a9291f5d
       // ... (合計6つの係数)
   }
   var _cos = [...]float64{
       -1.13585365213876817300E-11, // 0xbda8fa49a0861a9b
       2.08757008419747316778E-9,   // 0x3e21ee9d7b4e3f05
       // ... (合計6つの係数)
   }
   

5. Cos関数の再実装: Cos関数が、Cephesのアルゴリズムに基づいて全面的に書き直されました。これには、特殊ケースのハンドリング、高精度な範囲還元、そして_cosおよび_sin係数を用いた多項式近似の適用が含まれます。

6. Sin関数の再実装: Sin関数も同様に、Cephesのアルゴリズムに基づいて全面的に書き直されました。Cos関数と非常に似た構造を持ちますが、符号の扱いと多項式近似の選択（jの値に基づく）が異なります。

コアとなるコードの解説

src/pkg/math/sin.goにおけるCos関数とSin関数の新しい実装

新しいCos関数とSin関数は、以下の主要なステップで構成されています。

1. 定数の定義: PI4A, PI4B, PI4Cは、π/4を3つのfloat64値に分割して表現したものです。これにより、π/4をより高い精度で扱うことが可能になります。M4PIは4/πの近似値です。これらの定数は、高精度な範囲還元のために使用されます。

2. 特殊値の処理: 入力xがNaN（非数）やInf（無限大）の場合、Goのmath.NaN()を返します。Sin関数では、xが±0の場合に±0を返すというIEEE 754の規定も満たしています。

3. 引数の正規化と符号の決定: 入力xの絶対値を取り、元のxが負であった場合はsignフラグをtrueに設定します。Cos関数は偶関数なので、Cos(-x) = Cos(x)であるため、Cos関数ではsignフラグは最終結果の符号には影響しません。しかし、Sin関数は奇関数なので、Sin(-x) = -Sin(x)であるため、signフラグが最終結果の符号を決定します。

4. 範囲還元（Range Reduction）: j := int64(x * M4PI): xをπ/4で割ったおおよその整数部分jを計算します。 y := float64(j): jをfloat64に変換します。 if j&1 == 1 { j += 1; y += 1 }: これは、jが奇数の場合（つまりxがπ/4の奇数倍に近い場合）、jを偶数に調整することで、後続の計算でzが0に近くなるようにするための最適化です。これにより、多項式近似の精度が向上します。 j &= 7: jを8で割った余りを取ります。これは、xが2π（8 * π/4）のどの「八分円」（octant）に属するかを特定するためです。これにより、三角関数の周期性を利用して、計算を[0, 2π]の範囲にマッピングします。 if j > 3 { j -= 4; sign = !sign }: jが4以上の場合、xがπより大きいことを意味します。この場合、xをπだけずらしてjを調整し、結果の符号を反転させます。これは、sin(x + π) = -sin(x)やcos(x + π) = -cos(x)といった関係を利用したものです。 if j > 1 { sign = !sign }: Cos関数にのみ存在するこの行は、jが2または3の場合（つまりxがπ/2からπの範囲にある場合）に符号を反転させます。これは、cos(x + π/2) = -sin(x)やcos(x + π)の関係を利用したものです。

   z := ((x - y*PI4A) - y*PI4B) - y*PI4C: これが「拡張精度モジュラ演算」の核心です。xからy * (PI4A + PI4B + PI4C)を引くことで、xをπ/4の倍数で割った余りzを非常に高い精度で計算します。複数のfloat64値でπ/4を表現し、段階的に減算を行うことで、通常の浮動小数点減算で発生しやすい桁落ちを防ぎ、zの精度を保ちます。

5. 多項式近似の選択と評価: zz := z * z: zの二乗を計算します。多項式近似は通常、x^2の項で構成されるため、事前に計算しておきます。 if j == 1 || j == 2: この条件は、還元されたzが[0, π/4]の範囲にあり、かつ元のxがπ/4からπ/2の範囲（またはそれに対応する象限）にある場合に真となります。この場合、sin(x)はcos(π/2 - x)として、cos(x)はsin(π/2 - x)として計算するのが精度上有利です。 y = 1.0 - 0.5*zz + zz*zz*((((((_cos[0]*zz)+_cos[1])*zz+_cos[2])*zz+_cos[3])*zz+_cos[4])*zz+_cos[5]): この式は、cos(z)のテイラー展開（またはチェビシェフ近似）の多項式を評価しています。係数_cosが使用されます。 else: 上記の条件が偽の場合、つまり還元されたzが[0, π/4]の範囲にあり、かつ元のxが[0, π/4]の範囲（またはそれに対応する象限）にある場合です。 y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5]): この式は、sin(z)のテイラー展開（またはチェビシェフ近似）の多項式を評価しています。係数_sinが使用されます。

6. 最終結果の符号適用: if sign { y = -y }: 最初に保存したsignフラグがtrueであれば、最終結果の符号を反転させます。

この新しい実装は、Cephes Math Libraryの知見を取り入れることで、特に大きな入力値に対する三角関数の計算精度を大幅に向上させています。高精度な範囲還元と、定義域の異なる多項式近似の使い分けが、その精度向上の鍵となっています。

関連リンク

• Go Issue #1564: https://code.google.com/p/go/issues/detail?id=1564 (古いGoのIssueトラッカーのリンクですが、現在はGitHubに移行しています)
• Go CL 5320056: https://golang.org/cl/5320056 (Goのコードレビューシステムにおける変更リスト)

参考にした情報源リンク

• Cephes Math Library: http://netlib.sandia.gov/cephes/
• sin.c from Cephes: http://netlib.sandia.gov/cephes/cmath/sin.c
• cos.c from Cephes: http://netlib.sandia.gov/cephes/cmath/cos.c
• IEEE 754 浮動小数点数: 浮動小数点数の表現と演算に関する国際標準。
• 数値計算における桁落ち: 浮動小数点演算で発生する誤差の一種。
• 多項式近似: 関数を多項式で近似する数学的手法。