Keyboard shortcuts

Press ← or → to navigate between chapters

Press S or / to search in the book

Press ? to show this help

Press Esc to hide this help

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

このコミットは、Go言語の標準ライブラリである math パッケージ内の src/pkg/math/atan.go ファイルに対する変更です。このファイルは、アークタンジェント（Atan）関数、およびそれに関連するアークサイン（Asin）とアークコサイン（Acos）関数の実装を含んでいます。これらの関数は、与えられた数値の逆三角関数値を計算するために使用されます。

コミット

このコミットは、Go言語の math パッケージにおける Atan、Asin、Acos 関数の精度を向上させることを目的としています。既存のテストでは捕捉できなかった特定の入力値において、旧アルゴリズムが許容範囲を超える誤差を生じていた問題を修正し、より高精度な計算結果を提供します。

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

https://github.com/golang/go/commit/35c3afdb62809e40422b5590d807418697ba8660

元コミット内容

commit 35c3afdb62809e40422b5590d807418697ba8660
Author: Charles L. Dorian <cldorian@gmail.com>
Date:   Sun Jun 24 19:39:07 2012 -0400

    math: improve Atan, Asin and Acos accuracy
    
    pkg/math/all_test.go tests Atan (and therefore Asin and Acos) to a
    relative accuracy of 4e-16, but the test vector misses values where
    the old algorithm was in error by more than that. For example:
    
    x            newError   oldError
    0.414215746  1.41e-16  -4.24e-16
    0.414216076  1.41e-16  -4.24e-16
    0.414217632  1.41e-16  -4.24e-16
    0.414218770  1.41e-16  -4.24e-16
    0.414225466  0         -5.65e-16
    0.414226244  1.41e-16  -4.24e-16
    0.414228756  0         -5.65e-16
    0.414235089  0         -5.65e-16
    0.414237070  0         -5.65e-16
    
    R=rsc, golang-dev
    CC=golang-dev
    https://golang.org/cl/6302093

変更の背景

Go言語の math パッケージに実装されている Atan (アークタンジェント)、Asin (アークサイン)、Acos (アークコサイン) 関数は、数値計算において高い精度が求められます。しかし、このコミット以前のバージョンでは、特定の入力値に対して、既存のテスト (pkg/math/all_test.go) で設定されていた相対精度 4e-16 を超える誤差が発生していました。

コミットメッセージに示されている例では、x の値が 0.4142... 付近で、旧アルゴリズムが −4.24e-16 から −5.65e-16 といった比較的大きな誤差を生じていたのに対し、新しいアルゴリズムでは誤差が 0 から 1.41e-16 に改善されていることが示されています。これは、旧アルゴリズムが特定の領域で数値的に不安定であったか、あるいは近似式の係数が最適ではなかったことを示唆しています。

このような精度問題は、科学技術計算、グラフィックス、物理シミュレーションなど、浮動小数点数の精度が結果に大きく影響するアプリケーションにおいて、誤った結果や不安定な挙動を引き起こす可能性があります。そのため、より堅牢で高精度な数学関数の実装が求められ、このコミットに至りました。

前提知識の解説

1. 逆三角関数 (Atan, Asin, Acos)

  • アークタンジェント (Atan): y = tan(x) の逆関数で、x のタンジェント値が y となる角度 x をラジアンで返します。通常、戻り値の範囲は [-π/2, π/2] です。
  • アークサイン (Asin): y = sin(x) の逆関数で、x のサイン値が y となる角度 x をラジアンで返します。戻り値の範囲は [-π/2, π/2] です。
  • アークコサイン (Acos): y = cos(x) の逆関数で、x のコサイン値が y となる角度 x をラジアンで返します。戻り値の範囲は [0, π] です。

これらの関数は、幾何学、物理学、工学など、様々な分野で角度や方向を計算するために不可欠です。Asin や Acos は、しばしば Atan を用いて実装されます。例えば、asin(x) = atan(x / sqrt(1 - x^2)) のような関係があります。そのため、Atan の精度向上は、Asin や Acos の精度向上にも直接的に寄与します。

2. 浮動小数点数精度と誤差

コンピュータにおける数値計算は、有限のビット数で実数を近似する浮動小数点数で行われます。この近似には常に誤差が伴います。

  • 丸め誤差 (Rounding Error): 浮動小数点数で表現できない実数を最も近い表現可能な値に丸める際に発生する誤差です。
  • 打ち切り誤差 (Truncation Error): 無限級数や反復計算を途中で打ち切る際に発生する誤差です。数学関数の近似計算でよく見られます。
  • 相対誤差 (Relative Error): 誤差の絶対値を真値で割ったもので、真値に対する誤差の割合を示します。コミットメッセージの 4e-16 は相対誤差の許容範囲を示しています。

数値計算のアルゴリズム設計では、これらの誤差が累積しないように、また特定の入力値で誤差が急増しないように、細心の注意を払う必要があります。特に、非常に小さい数と大きい数の演算、あるいはほぼ等しい数の引き算（桁落ち）は、精度問題を引き起こしやすいです。

3. Cephes Math Library

Cephes Math Libraryは、Stephen L. Moshierによって開発された、高品質な数学関数のC言語実装のコレクションです。これは、数値計算の分野で広く利用されており、その精度と堅牢性で知られています。多くの科学技術計算ライブラリやシステムで、そのアルゴリズムや実装が参考にされています。

Cephesライブラリの関数は、通常、多項式近似や有理関数近似（多項式の比）を用いて実装されており、特定の範囲で高い精度を達成するように設計されています。このコミットでは、Atan 関数の実装をCephesライブラリの atan.c に基づくものに置き換えることで、精度向上を図っています。

4. 有理関数近似 (Rational Function Approximation)

数学関数をコンピュータで計算する際、多くの場合、無限級数（テイラー級数など）や多項式、あるいは有理関数（多項式の比）を用いて近似します。

  • 多項式近似: P(x) = a_n x^n + ... + a_1 x + a_0 の形で関数を近似します。
  • 有理関数近似: R(x) = P(x) / Q(x) の形で関数を近似します。ここで P(x) と Q(x) は多項式です。

有理関数近似は、多項式近似よりも少ない項数で同等以上の精度を達成できる場合があり、特に広範囲で関数を近似する際に有効です。コミットメッセージのコメントにも「rational function of degree 4/5」と記載されており、これは分子が4次、分母が5次の多項式による有理関数近似を意味します。

技術的詳細

このコミットにおける Atan 関数の精度向上は、主に以下の2つの関数 xatan と satan の実装変更によって実現されています。

1. xatan 関数の変更

xatan 関数は、Atan 計算の中核をなす部分で、特定の狭い範囲（旧版では [-0.414..., +0.414...]、新版では [0, 0.66]）におけるアークタンジェント値を計算します。

  • 旧実装: Hart & Cheneyの係数 #5077 (19.56D) に基づく多項式近似を使用していました。

    sq := arg * arg
value := ((((P4*sq+P3)*sq+P2)*sq+P1)*sq + P0)
value = value / (((((sq+Q4)*sq+Q3)*sq+Q2)*sq+Q1)*sq + Q0)
return value * arg

    このコードは、arg の2乗 sq を用いた多項式 P(sq) と Q(sq) を計算し、P(sq)/Q(sq) * arg の形で近似を行っていました。

  • 新実装: Cephes Math Libraryの atan.c に基づく有理関数近似に完全に置き換えられました。新しい係数 P0 から P4、Q0 から Q4 が導入されています。

    z := x * x
z = z * ((((P0*z+P1)*z+P2)*z+P3)*z + P4) / (((((z+Q0)*z+Q1)*z+Q2)*z+Q3)*z + Q4)
z = x*z + x
return z

    この新しい有理関数は、x の2乗 z を用いて計算され、x + x^3 P(x^2)/Q(x^2) の形式（コミットコメントでは x + x**3 P(x)/Q(x) と記載）で近似されます。Cephesライブラリは、その精度と堅牢性で知られており、この変更により、より広範囲の入力値に対して安定した高精度な結果が得られるようになりました。特に、旧アルゴリズムで誤差が大きかった x の値に対して、大幅な改善が見られます。

2. satan 関数の変更

satan 関数は、入力された引数（正の値）を xatan 関数が処理できる範囲に「範囲縮小 (range reduction)」する役割を担っています。

  • 旧実装: 引数を [0, Sqrt2-1] (約 [0, 0.414...]) の範囲に縮小していました。

    • arg < Sqrt2-1 の場合: そのまま xatan(arg) を呼び出し。
    • arg > Sqrt2+1 の場合: Pi/2 - xatan(1/arg) を計算。これは atan(x) + atan(1/x) = Pi/2 の関係を利用しています。
    • それ以外の場合: Pi/4 + xatan((arg-1)/(arg+1)) を計算。これは atan(x) = Pi/4 + atan((x-1)/(x+1)) の関係を利用しています。

  • 新実装: 範囲縮小のターゲットが [0, 0.66] に変更されました。また、Morebits と Tan3pio8 という新しい定数が導入されています。

    const (
    Morebits = 6.123233995736765886130e-17 // pi/2 = PIO2 + Morebits
    Tan3pio8 = 2.41421356237309504880      // tan(3*pi/8)
)
if x <= 0.66 {
    return xatan(x)
}
if x > Tan3pio8 {
    return Pi/2 - xatan(1/x) + Morebits
}
return Pi/4 + xatan((x-1)/(x+1)) + 0.5*Morebits

    Morebits は Pi/2 のより正確な表現を補正するための微小な値であり、Tan3pio8 は tan(3π/8) の値です。これらの定数と新しい範囲縮小ロジックは、Cephesライブラリの設計思想に従っており、xatan の新しい入力範囲 [0, 0.66] に合わせて調整されています。これにより、Atan 関数全体としての精度と安定性が向上します。

3. Atan, Asin, Acos への影響

このコミットは Atan 関数の内部実装を直接変更していますが、Asin と Acos 関数も Atan を利用して実装されていることが多いため、間接的にこれらの関数の精度も向上します。例えば、Go言語の math パッケージでは、Asin は atan2(x, Sqrt(1-x*x)) を、Acos は atan2(Sqrt(1-x*x), x) を利用して実装されています（atan2 は Atan を基盤とする）。したがって、Atan の精度が向上すれば、これらの関数もより正確な結果を返すようになります。

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

src/pkg/math/atan.go ファイルにおいて、主に以下の部分が変更されました。

  1. コメントの追加と変更:

    • ファイルの冒頭に、Cephes Math Libraryの atan.c に基づく実装であること、およびそのライセンス情報に関する詳細なコメントが追加されました。
    • xatan 関数のコメントが「evaluates a series valid in the range [0, 0.66]」に変更されました。

  2. xatan 関数の係数と計算ロジックの変更:

    • 旧係数 P4 から P0、Q4 から Q0 が、新しい値に置き換えられました。
    • value の計算ロジックが、新しい有理関数近似の形式に完全に変更されました。
    --- a/src/pkg/math/atan.go
+++ b/src/pkg/math/atan.go
@@ -6,51 +6,92 @@ package math

 /*
 	Floating-point arctangent.
-
-	Atan returns the value of the arctangent of its
-	argument in the range [-pi/2,pi/2].
-	There are no error returns.
-	Coefficients are #5077 from Hart & Cheney. (19.56D)
 */

-// xatan evaluates a series valid in the
-// range [-0.414...,+0.414...]. (tan(pi/8))
-func xatan(arg float64) float64 {
+// The original C code, the long comment, and the constants below were
+// from http://netlib.sandia.gov/cephes/cmath/atan.c, available from
+// http://www.netlib.org/cephes/cmath.tgz.
+// The go code is a version of the original C.
+//
+// atan.c
+// Inverse circular tangent (arctangent)
+//
+// SYNOPSIS:
+// double x, y, atan();
+// y = atan( x );
+//
+// DESCRIPTION:
+// Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
+//
+// Range reduction is from three intervals into the interval from zero to 0.66.
+// The approximant uses a rational function of degree 4/5 of the form
+// x + x**3 P(x)/Q(x).
+//
+// ACCURACY:
+//                      Relative error:
+// arithmetic   domain    # trials  peak     rms
+//    DEC       -10, 10   50000     2.4e-17  8.3e-18
+//    IEEE      -10, 10   10^6      1.8e-16  5.0e-17
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+//    Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+//   The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+//   Stephen L. Moshier
+//   moshier@na-net.ornl.gov
+
+// xatan evaluates a series valid in the range [0, 0.66].
+func xatan(x float64) float64 {
 	const (
-		P4 = .161536412982230228262e2
-		P3 = .26842548195503973794141e3
-		P2 = .11530293515404850115428136e4
-		P1 = .178040631643319697105464587e4
-		P0 = .89678597403663861959987488e3
-		Q4 = .5895697050844462222791e2
-		Q3 = .536265374031215315104235e3
-		Q2 = .16667838148816337184521798e4
-		Q1 = .207933497444540981287275926e4
-		Q0 = .89678597403663861962481162e3
+		P0 = -8.750608600031904122785e-01
+		P1 = -1.615753718733365076637e+01
+		P2 = -7.500855792314704667340e+01
+		P3 = -1.228866684490136173410e+02
+		P4 = -6.485021904942025371773e+01
+		Q0 = +2.485846490142306297962e+01
+		Q1 = +1.650270098316988542046e+02
+		Q2 = +4.328810604912902668951e+02
+		Q3 = +4.853903996359136964868e+02
+		Q4 = +1.945506571482613964425e+02
 	)
-	sq := arg * arg
-	value := ((((P4*sq+P3)*sq+P2)*sq+P1)*sq + P0)
-	value = value / (((((sq+Q4)*sq+Q3)*sq+Q2)*sq+Q1)*sq + Q0)
-	return value * arg
+	z := x * x
+	z = z * ((((P0*z+P1)*z+P2)*z+P3)*z + P4) / (((((z+Q0)*z+Q1)*z+Q2)*z+Q3)*z + Q4)
+	z = x*z + x
+	return z
 }

  3. satan 関数の範囲縮小ロジックと定数の変更:

    • Morebits と Tan3pio8 という新しい定数が追加されました。
    • 範囲縮小の条件と計算式が変更され、xatan の新しい入力範囲 [0, 0.66] に対応するように調整されました。
    --- a/src/pkg/math/atan.go
+++ b/src/pkg/math/atan.go
@@ -59,18 +100,24 @@
 }

 // satan reduces its argument (known to be positive)
-// to the range [0,0.414...] and calls xatan.
-func satan(arg float64) float64 {
-	if arg < Sqrt2-1 {
-		return xatan(arg)
+// to the range [0, 0.66] and calls xatan.
+func satan(x float64) float64 {
+	const (
+		Morebits = 6.123233995736765886130e-17 // pi/2 = PIO2 + Morebits
+		Tan3pio8 = 2.41421356237309504880      // tan(3*pi/8)
+	)
+	if x <= 0.66 {
+		return xatan(x)
 	}
-	if arg > Sqrt2+1 {
-		return Pi/2 - xatan(1/arg)
+	if x > Tan3pio8 {
+		return Pi/2 - xatan(1/x) + Morebits
 	}
-	return Pi/4 + xatan((arg-1)/(arg+1))
+	return Pi/4 + xatan((x-1)/(x+1)) + 0.5*Morebits
 }

 // Atan returns the arctangent of x.
 //
 // Special cases are:
-//	Atan(±0) = ±0
-//	Atan(±Inf) = ±Pi/2
+//      Atan(±0) = ±0
+//      Atan(±Inf) = ±Pi/2
 func Atan(x float64) float64

コアとなるコードの解説

xatan 関数の解説

xatan 関数は、アークタンジェントの計算において、最も重要な近似計算を行う部分です。旧バージョンではHart & Cheneyの係数を用いた多項式近似が使われていましたが、新バージョンではCephes Math Libraryの atan.c に基づく、より洗練された有理関数近似が導入されました。

新しい xatan の計算式は以下の形式です。 z = x^2 z = z * P(z) / Q(z) result = x * z + x

ここで、P(z) と Q(z) はそれぞれ5次の多項式です（P0からP4、Q0からQ4がその係数）。この有理関数近似は、特定の範囲 [0, 0.66] において、非常に高い精度でアークタンジェント関数を近似するように設計されています。特に、x が小さい値のときに x に近い値になるように x*z + x の形をとることで、原点付近での精度を保ちつつ、より広い範囲で効率的な近似を実現しています。

satan 関数の解説

satan 関数は、Atan 関数に渡された引数 x を、xatan 関数が最も効率的かつ高精度に処理できる狭い範囲 [0, 0.66] に変換（範囲縮小）する役割を担っています。

この関数は、入力 x の値に応じて3つのケースに分岐します。

  1. x <= 0.66 の場合:

    • 引数がすでに xatan の処理範囲内にあるため、そのまま xatan(x) を呼び出します。

  2. x > Tan3pio8 (約 2.414) の場合:

    • このケースは、x が比較的大きい場合に対応します。atan(x) + atan(1/x) = Pi/2 という三角関数の恒等式を利用して、Pi/2 - xatan(1/x) を計算します。1/x は x が大きいほど小さくなるため、xatan の処理範囲に収まります。
    • Morebits (6.123233995736765886130e-17) は、Pi/2 の浮動小数点表現の精度を補正するための微小な値です。これにより、計算結果の精度がさらに向上します。

  3. それ以外の場合 ( 0.66 < x <= Tan3pio8 ):

    • この中間的な範囲では、atan(x) = Pi/4 + atan((x-1)/(x+1)) という恒等式を利用します。 (x-1)/(x+1) は x が 0.66 から Tan3pio8 の範囲にあるとき、xatan の処理範囲 [0, 0.66] に収まるように設計されています。
    • ここでも 0.5*Morebits が追加されており、Pi/4 の精度補正に寄与しています。

これらの範囲縮小戦略と、xatan の高精度な有理関数近似を組み合わせることで、Atan 関数は広範囲の入力値に対して、非常に高い精度と安定性を提供できるようになりました。

関連リンク

参考にした情報源リンク