57

第 6 章 初等関数の計算
しかし著者がとくに強調したいのは，初等関数の数値計算というきわめて「初等的」な，狭い，そしてすでに研究しつくされていると思われる分野でさえも，対象を見る目をかえていじくれば，単なる落ち穂拾いというのではすまされないようなおもしろい話題が，たくさんあるということである。
一松信「初等関数の数値計算」(教育出版)

本章では，多項式関数 p(x) =Pni
i=0 ai x ，平方根√x，三角関数 sin x, cos x，指数関数 exp(x)，(自然) 対数関数 log x の計算方法について解説する。三角関数，指数関数，対数関数については Taylor
展開に基づいた計算方法を解説するが，現実に使用されているアルゴリズムに比べて効率が悪いことが知られている。従って，これらの方法はあくまで数値計算のエッセンスを理解するための一例題としてとらえて欲しい。

6.1 Horner 法による多項式関数の計算
定義 6.1.1 (多項式関数)
以下のような関数を，複素係数の多項式関数 (polynomial function) と呼ぶ。

p(x) =nX

ai xi (ai ∈ C, x ∈ C)

(6.1)

i=0

多項式関数をプログラム言語で扱うには，係数のみを一次元配列に格納しておくだけで済む。特定の実数値 x = α における多項式関数の値 p(α) が必要になれば，以下の Horner 法を使って計算すればよい。このアルゴリズムは四則演算の範囲内で実行可能である。
アルゴリズム 1 (Horner 法)

1. val := an とする。
2. 以下 i = n − 1, ..., 1, 0 に対し，以下を繰り返す。
(a) val := val × x + ai
3. p(x) = val となる。
前述したように，数値計算においては一般に，加減算より乗除算の方が多くの計算時間を必要とする。Horner 法は乗算の回数が最小で済むため，多項式関数の値を評価する標準的な方法となっている。

第6章

58

初等関数の計算問題 6.1.1n 次多項式関数 p(x) の値 p(α) を求めるのに必要な計算量を求めよ。またこの多項式関数の微係数p′ (α) を計算するためには，アルゴリズム 1 をどのように変更すればよいか？

6.2 Newton 法に基づく平方根の計算
ここでは正の引数 a ∈ R に対する平方根√a の計算方法を考える。平方根の計算方法は様々なものが提案されているが，最も頻繁に取上げられるのは Newton 法に基づく方法である。
√2 次方程式 x2 − a = 0 の正の解は a である。従って，この方程式の解を求めるアルゴリズムが√あれば，それが a を求めるアルゴリズムとなる。

Newton 法について第 12 章で詳しく述べるが，この場合は次のような数列 {xi }∞i=0 を生成するアルゴリズムになる。このように同じ計算を繰り返し実行して近似値の精度を高めていく手法を総称して反復法 (iteration method, iterative method) と呼ぶ。これを使用する際には停止則 (第 5 章参照)に注意を払う必要がある。
アルゴリズム 2 (Newton 法による平方根の計算)

1. 初期値 x0 を決める。例えば x0 := a とする。
2. i = 0, 1, ... に対して次の計算を行う。
1axi+1 :=xi +2xi
このアルゴリズムに基づいた√

!
(6.2)

2 を計算した例を以下に示す。

√
例題 6.2.1 ( 2 = 1.41421356237309504 · · · の計算)
アルゴリズム 2 を使って，初期値を x0 = 2 としたもの (2 列目) と，x0 = (1 + 2)/2 としたもの (3 列目) を，IEEE754 倍精度浮動小数点計算を用いて計算した結果が以下の表である。また相対誤差をプロットしたものを図 6.1 に示す。

xix0 = 2 の場合

x0 = (1 + 2)/2 の場合 (相対誤差)(相対誤差)x0−12.00000000000000000(4.14 × 10 )1.50000000000000000(6.07 × 10−2 )x11.50000000000000000(6.07 × 10−2 )1.41666666666666674(1.73 × 10−3 )x2

1.41666666666666674(1.73 × 10−3 )1.41421568627450989(1.50 × 10−6 )x3

1.41421568627450989(1.50 × 10−6 )1.41421356237468987(1.13 × 10−12 )x4−121.41421356237468987(1.13 × 10

x0 = (1 + 2)/2 とした方がより早く真値√)

2 に接近していることが分かる。

問題 6.2.1

a = 10, 100, 1000 に対してアルゴリズム 2 を適用し，10 進 5 桁程度の近似解を求めよ。

6.3. Taylor 展開に基づく初等関数の計算

59

11e-05

初期値 = 2相対誤差 1e-10 初期値 = (1 + 2) /2
1e-15

1e-200123

反復回数45

図 6.1: Newton 法による近似値の相対誤差6.3 Taylor 展開に基づく初等関数の計算
まず微分積分の復習として Talyor の定理を示す。
定理 6.3.1 (Taylor の定理, Taylor 展開, Maclaurin 展開)
実関数 f (x) が閉区間 [a, b] で m 回連続微分可能かつ開区間 (a, b) で m + 1 回微分可能である時，∀x ∈ (a, b) に対して
f (x) = f (a) +f ′ (a)f ′′ (a)f (m) (a)f (m+1) (a + θ(x − a))(x − a) +(x − a)2 + · · · +(x − a)m +(x − a)m+1 (6.3)
1!
2!
m!
(m + 1)!

を満足する 0 < θ < 1 が存在する。特に無限回微分可能であれば

f (x) = f (a) +f ′ (a)f ′′ (a)f (m) (a)(x − a) +(x − a)2 + · · · +(x − a)m + · · ·
1!
2!
m!

(6.4)と x に関する無限級数で表わすことができ，この式 (6.4) を関数 f (x) の Taylor 展開 (Taylor expansion)と呼ぶ。また特に a = 0 の時を Maclaurin 展開と呼ぶ。
この定理は平均値の定理を繰り返し適用することで証明出来る。極めて綺麗な定理であり，応用も範囲も広い。最初の節で述べたように，多項式関数は Horner 法によって四則演算だけでその値を計算することが出来る。即ち，この定理の条件に当てはまる関数で，微係数 f (i) (a) が判明しているものであれば，その関数の近似値を (6.3) 式の多項式関数によって，四則演算の範囲で得ることが可能となる。そして，一般に初等関数と呼ばれる三角関数，指数関数，対数関数は全てこれらの条件に当てはまる上，微係数も容易に求めることが出来るのである。
但し，章の最初にも述べたように，Taylor 展開 (Maclaurin) 展開をそのまま適用して初等関数の値を計算するという手法は計算量の観点からあまり得策ではないことが知られている。実用に供

第6章

60

初等関数の計算されている初等関数のアルゴリズムはこれとは別のもの (最良近似多項式，有理近似式，CORDIC等) であることを付け加えておく。
例題 6.3.2 (代表的な初等関数の Maclaurin 展開)
一般に初等関数 (elementary functions) と呼ばれる，三角関数，指数関数，対数関数は R 全体もしくは特定の区間で無限回連続微分可能である。従って，Taylor 展開 (Maclaurin 展開) が存在する。

exp(x) =
sin x

=

cos x

=log(1 + x) =x
x2
xn++ ··· ++ ···1! 2!
n!
x(2n−1)x3+ · · · + (−1)n−1+ ···x−
3!
(2n − 1)!
x2
x2n1−+ · · · + (−1)n+ ···
2!
(2n)!
x2
xn+1
x−+ · · · + (−1)n+ · · · (ここで (1 + x) > 0)
2n+1
1+

(6.5)
(6.6)
(6.7)
(6.8)

但し，これらの初等関数の Maclaurin 展開式を多項式関数として計算するには，引数 x に応じた配慮をする必要がある。以下，exp(x), sin(x), log(x) について具体的に計算方法を詰めていくことにする。

6.3.1 e = exp(1) の計算と誤差解析
丸め誤差が実数を有限桁の浮動小数点数で近似した結果生じた誤差であるのに対し，打ち切り誤差は無限級数や極限値のような無限回の演算を必要とする解析表現を，有限回の演算で打ち切る(truncate) ことによって生じる誤差である。丸め誤差は数値によって変動し精密な予測が難しいのに対し，打ち切り誤差は解析表現が明らかであれば，それに基づいて予測することが可能である。
故に，打ち切り誤差は理論誤差とも呼ばれる。e = exp(1) の計算を例に，この打ち切り誤差について見ていくことにしよう。

e は指数関数 exp(x)(= e x ) の Maclaurin 展開式 (6.5) 式によって
e=1+

1
1
1++ ··· ++ ···1! 2!
n!という無限級数の形で表現される。しかし，いかに高速なコンピュータといえども無限級数を計算することはできないため，どこかの項 1/m! で計算を打ち切る必要がある。この項までの有限和をêm と書くことにする。すなわち，êm = 1 +111++ ··· +1! 2!
m!

である。この時，打ち切り誤差はe − êm =11++ ···(m + 1)! (m + 2)!となる。
この形では評価が難しいので有限和の Maclaurin 展開公式

e=1+1
1
1exp(θ)++ ··· ++1! 2!
m! (m + 1)!

6.3. Taylor 展開に基づく初等関数の計算

61を用いることにする。ここで θ は 0 < θ < 1 となる定数である。
これを用いると打ち切り誤差はe − êm =exp(θ)(m + 1)!となる。右辺の絶対値を取ればexp(θ)e≤(m + 1)!
(m + 1)!となるので，相対打ち切り誤差を取るとe − êm1≤e(m + 1)!

(6.9)となり，m が決まれば打ち切り誤差の上限を評価することが可能となる。
一般に，打ち切り誤差は計算回数さえ増やせば減らすことができるが，使用する浮動小数点数の丸め誤差の最小単位より過度に小さくしても，コンピュータ資源の無駄遣いにしかならない。実際，IEEE754 倍精度計算を行い，項数を増やしつつ計算してその相対誤差をプロットしてみると図6.2 のようになり，20 項以上取ってもマシンイプシロン ε M 以上の精度を得ることは出来ないことが分かる。

相1対誤
差
1e-05

1e-10

1e-15010

項数2030

図 6.2: Maclaurin 展開に基づく exp(1) の近似値の相対誤差
例えば，10 進 7 桁の浮動小数点数を用いて e を計算するのであれば，先の評価式 (6.9) を用いて

e − êm11≤≈ · 10−6e(m + 1)! 2程度になる m まで計算するのが最適と言える。この場合，

11≈ 2.8 × 10−6 ,≈ 2.8 × 10−7(8 + 1)!
(9 + 1)!

第6章

62

初等関数の計算であるから，m = 9 ，すなわち111++ ··· +1! 2!
9!
程度まで計算しておけば十分である。実際に計算してみると
ê9 = 1 +ê9 = 2.71828152557 · · ·であり，下線部の 7 桁分が真値と一致していることが分かる。
問題 6.3.1e を 10 進 15 桁の精度を得るために必要な項数 m を求め，実際に êm を計算せよ。

6.3.2 exp(x) の計算exp(x) を Maclaurin 展開級数 (6.5) に基づいて計算するには次の 2 点を勘案しなくてはならない。
x < 0 の場合 引数 x が負の場合，exp(x) の絶対値は小さくなる。特に |x| が大きくなると，Maclaurin
展開式の項の絶対値 |xi /i!| も大きくなる。しかしそれらの和を取った結果の絶対値は小さくなるのだから，これは和の計算において maxi |xi /i!| と | exp(x)| の order の差だけ「桁落ち」が起きることを示している。従って，引数が負の時は，exp(x) = 1/ exp(−x) という関係を使い，正の引数 exp(−x) を計算し，しかる後にその逆数を取る，という手順を取る必要がある。

x >> 1 の場合 理論上，この無限級数は任意の x ∈ R について収束することになっているが，数値計算上は x >> 1 の場合収束が遅くなり，必要な項数が増えてしまう。即ち xn /n! << 1 となる n ∈ N が大きくなってしまうことになる。従って，無限級数の計算を行う x の範囲を，例えば 0 < x − ⌊x⌋ < 1 に限定し，それを越える分については別途 exp(⌊x⌋) を計算して掛け合わせるようにすればよい。
以上をまとめると，次のようなアルゴリズムとなる。
アルゴリズム 3 (exp(x) の計算)

(a) 0 ≤ x ≤ 1 の時には Maclaurin 展開級数 (6.5) を使用する。但し x = 0 の時は 1 を， x = 1 の時は 2.718281 · · · を返すようあらかじめ定数を設定しておく。

(b) x > 1 の時は， x′ := x − ⌊x⌋ として，exp(x′ ) exp(⌊x⌋) を計算する。当然 exp(⌊x⌋) の部分は定数2.71828 · · · の ⌊x⌋ 乗として計算する。
(c) x < 0 の時は exp(|x|) を (a) 及び (b) を用いて計算し，その逆数 1/ exp(|x|) を取る。
結局，実際に計算するのは 0 < x < 1 の範囲における exp(x) の Maclaurin 展開級数である。では，どの程度の項数を取れば「収束」し必要な精度を得られるのか。(6.5) の右辺を m 項で打ち切った
有限和を

dm (x) =expmX
1i!
i=0xi6.3. Taylor 展開に基づく初等関数の計算

63と書くことにすれば，e の計算同様，打ち切り誤差はdm (x) =exp(x) − expとなるから，右辺はexp(θx) m+1x
(m + 1)!

exp(θx) m+1 e · exp(x) m+1x≤x(m + 1)!
(m + 1)!と抑えられるので，相対打ち切り誤差を取り

dm (x)
exp(x) − exp
e · xm+1≤exp(x)(m + 1)!

(6.10)に基づき，0 < x < 1 が決まれば打ち切り誤差の上限を評価することが可能となる。
問題 6.3.2

10 進 15 桁の有効桁を持つよう exp(3.2) を計算せよ。(6.10) 式を用いて第何項まで計算すればよいかも予め評価し，実際に計算してその相対誤差を求めよ。

6.3.3

sin x の計算

正弦関数 sin x の計算も，なるべく小さい x を使って Maclaurin 展開式 (6.6) を計算できるように配慮する必要がある。そのため，sin x の性質を利用して次のように計算すると良い。
アルゴリズム 4 (sin(x) の計算)

(a) 0 ≤ x ≤ π/2 の時には Maclaurin 展開級数 (6.6) を使用する。但し x = 0 の時は 0 を， x = π/2の時は 1 を返すようあらかじめ定数を設定しておく。
(b) π/2 < x ≤ π の時は， x′ := π − x として，sin x := sin x′ を計算する。
(c) π < x ≤ 2π の時は， x′ := x − π として，sin x := − sin x′ を計算する。
(d) x > 2π の時は x′ := x − 2π · ⌊x/2π⌋ として，(a)〜(c) を用いて sin x := sin x′ を計算する。
(e) x < 0 の時は sin(−x) の値を (a)〜(d) を使って求め，sin x := − sin(−x) とする。
但しこれでも不十分な部分がある。図 6.3 は sin x の値 (実線) と，Maclaurin 展開式 (6.6) を用いた近似値の相対誤差 (破線) をプロットしたものであるが，ちょうど sin x ≈ 0 となる部分で，相対誤差が著しく増大しているのが分かる。これは展開式の計算で桁落ちが発生していることを示している。従って，この部分の近似値を計算する際には特別の配慮をする必要がある。
問題 6.3.3

1. 三角関数 cos x を Maclaurin 展開を元に計算する場合，どのような配慮が必要か，sin x の例を元に述べよ。

2. sin 1.5, cos 1.5 を IEEE754 倍精度程度 (10 進 15〜16 桁) 求めるには最大何項まで足し込んでいく必要があるか？

第6章

64

初等関数の計算

sin(x)
相対誤差

sin(x)

11

相対
誤
差

1e-05

0.5

1e-100-0.5

1e-15-1-50x5

1e-20

図 6.3: sin x と Maclaurin 展開に基づく sin x の相対誤差6.3.4 log x の計算
log x の計算は，収束が遅い展開式を如何に改善させるか，という話題が出る所でよく取り上げられる。先に挙げた log(x + 1) の Maclaurin 展開式 (6.8) の項を，exp(x) や sin x のそれを比べてみると，分母の階乗がない分，高次項の係数は log(1 + x) の方がずっと大きい。そのため，収束するために必要となる項数も多くなる。また収束半径もごく限られた範囲に留まるため， (6.8) 式は実用性に乏しい。
そこで，収束を早める工夫が提案されている。高次項をぐっと小さくさせるために展開式を変更するのである。例えば

log x2=[展開式 2] log x=[展開式 1]!
!3!2n−1x−1
1 x−1
1x−1++ ··· ++ ···x+1
3 x+12n − 1 x + 1!
!2!nx−1
1 x−1
1 x−1++ ··· ++ ···x2xnx(6.11)(6.12)

等が有名である。収束範囲は展開式 1 の方が広く，収束も早いことは一目瞭然であろう。実際，これらの展開式を使って log(2) の近似値を求めてみると，展開式 1 の方がずっと早く収束していることが分かる (図 6.4)。
問題 6.3.4

1. a ∈ R, a > 0 に対して，a x を計算する方法を考えよ。
2. 展開式 1 が展開式 2 よりも収束が早い理由を説明せよ。
3. log 100 を展開式 1 を使って IEEE754 倍精度程度 (10 進 15〜16 桁) 求めたい時，何項目まで足し込めばよいか？

6.4. その他の関数

65

相1対誤
差

展開式21e-05

1e-10

展開式11e-15

1e-2005

項数1015

図 6.4: log(2) の近似値の相対誤差

6.4

その他の関数以上で取り上げた初等関数を初めとして，実用上必要な関数は数多くある。ここでは C 言語の標準規格である C99[14] と，フリーのコンパイラである GCC に規定されている関数群 (表 6.1, 6.2)のうち，あまり馴染みのない関数を簡潔に紹介する。

6.4.1

誤差関数

誤差関数 erf(x) は，次の式で定義される。

2erf(x) = √πZ xexp(−t2 )dt(6.13)0

これはちょうど，確率密度関数 f (t) が1f (t) = √ exp(−t2 )πであるような連続型確率変数 X における，確率 P(−x ≤ X ≤ x)
Z x1
P(−x ≤ X ≤ x) = √exp(−t2 )dt
π −xZ x2= √exp(−t2 )dt
π 0と等しい。これは Gauss が天文観測値の誤差の精度を推定するために導入した関数で [44]，そのためこのように命名されている。

第6章

66

初等関数の計算表 6.1: C99 に規定されている数学関数の一部
関数名C99 における型指定

変数の範囲

数式表記

逆三角関数

acos(x)x ∈ [−1, 1]

cos−1 (x)asin(x)atan(x)x ∈ [−1, 1]

sin−1 (x)tan−1 (x)三角関数

cos(x)
sin(x)
tan(x)

逆双曲線関数acosh(x)asinh(x)atanh(x)双曲線関数

cosh(x)sinh(x)tanh(x)cosh(x)sinh(x)tanh(x)

指数関数

exp(x)
exp2(x)exp(x) = e x
2x

対数関数

log(x)x ∈ (0, ∞)log x = loge x = ln x

log10(x)log2(x)x ∈ (0, ∞)x ∈ (0, ∞)平行根

sqrt(x)x ∈ [0, ∞)立方根

cbrt(x)log10 x = lg x
log2 x√x√3x = x1/3

べき乗

pow(x, y)

xy

誤差関数

erf(x)
erfc(x)erf(x)
1 − erf(x)ガンマ関数

tgamma(x)Γ(x)

対数ガンマ関数

lgamma(x)log |Γ(x)|cos(x)
sin(x)
tan(x)
cosh−1 (x)sinh−1 (x)tanh−1 (x)x ∈ [1, ∞)x ∈ [−1, 1]

表 6.2: GCC 3.2 で規定されている数学関数
関数名型指定数式表記

第一種 Bessel 関数

j0(x)
j1(x)

J0 (x)J1 (x)

jn(int n, x)

Jn (x)

y0(x)
y1(x)

Y0 (x)Y1 (x)yn(int n, x)

Yn (x)

第ニ種 Bessel 関数

6.4. その他の関数

67

誤差関数を無限級数で表現すると

2 X (−1)i x2i+1erf(x) = √π i=0 i!(2i + 1)∞となることが知られている [18]。

6.4.2

ガンマ関数

ガンマ関数 Γ(x) は，次の式で定義される。
Z ∞Γ(x) =t x−1 exp(−t)dt(6.14)

0

統計学では，この関数によって確率密度関数が規定されるガンマ分布を取り扱う。また，これを用いて定義されるベータ関数 B(x, y)B(x, y) =Γ(x)Γ(y)Γ(x + y)も，統計ではよく利用される (ベータ分布)。
ガンマ関数は，対数ガンマ関数 log Γ(x) を，次の漸近展開 (x → ∞ の時収束する無限級数)∞X1
1B2ilog Γ(x)〜(x − ) log x − x + log(2π) +222i(2i − 1)x2i−1i=1を用いて計算し，Γ(x) = exp(log Γ(x)) として求める [31]。 x は次の漸化式

Γ(x + 1) = xΓ(x)を用いて，なるべく大きくしてから漸近展開の計算を行う。ここで Bi は Bernoulli 数で，展開式∞
Xtti=Biexp(t) − 1 i=0 i!に基づいて定義される有理数1 である。

6.4.3

Bessel 関数

Bessel の常微分方程式
d2 ydy+ x + (x2 − n2 )y = 0 (n ∈ Z)
2dxdxの解として定義されるのが，Bessel 関数である。この解は複数あり，
x2lim Jn (x) = 0n→∞1 B = 1, B = −1/2, B = 1/6, B = 0, ...
0123(6.15)

第6章

68

初等関数の計算となる解 Jn (x) を，単なる Bessel 関数，あるいは第一種 Bessel 関数と呼ぶ。また

lim Yn (x) = ∞n→∞となる解 Yn (x) を，第二種 Bessel 関数，あるいは Weber 関数と呼ぶ。
第ー種，第二種 Bessel 関数はどちらもJn+1 (x) =Yn+1 (x) =2nJn (x) − Jn−1x2n
Yn (x) − Yn−1xという漸化式が成立し，これに基づいて計算を行う。但し，Jn (x) は n が大きい方から小さい方へ，Yn (x) は逆に n が小さい方から大きい方へと計算しなければならないことが知られている [28]。

演習問題
1. Newton 法を用いて，a1/3 を求めるための漸化式は次のようになる。
!
1axi+1 :=2xi + 23xi
これを使って 31/3 = 1.442249570307408... を求めよ。

2. ユーザが x の値を入力した時，cos x の値を求めるプログラムを作りたい。この時，次の問に答えよ。なお π は円周率とする。
(a) cos x の Maclaurin 展開はcos x =∞
Xx2i(−1)i(2i)!
i=0である。0 < x < π/4 である時，cos x の値に含まれる打ち切り誤差が 10−5 以下になるよ
うにするには，第何項まで計算する必要があるか？

(b) α = π/7 を，10 進 5 桁の浮動小数点数 eα に丸めた。近似値 eα を書け。またその近似値に含まれる絶対誤差と相対誤差もあわせて求めよ。

(c) cos eα を求めよ。
3. sinh x, cosh x, tanh x を計算するプログラムを作りたい。
(a) sinh x = (exp(x) − exp(−x))/2, cosh x = (exp(x) + exp(−x))/2, tanh x = sinh x/ cosh x という関係式に基づいて計算するプログラムを作れ。
(b) 展開式
sinh x=cosh x=に基づいてプログラムを作れ。

x2n+1
x3 x5++ ··· ++ ···3! 5!
(2n + 1)!
x2 x4x2n
1+++ ··· ++ ···
2! 4!
(2n)!
x+

6.4. その他の関数

694. 逆三角関数 sin−1 x, cos−1 x, tan−1 x を計算するプログラムを作れ。
5. 円周率 π を計算するプログラムを作れ。(Hint: π/4 = tan−1 1 である。) また，多倍長計算可能なソフトウェア (Mathematica 等) を使って，円周率 π を 1000 桁，10000 桁, 100000 桁・・・計算した際に要した時間を計測し，その増加率を求めよ。また，それだけの桁数を求めるた
めに必要となる項数はどの程度か？

6. erfc(x) は，次のように定義されることを示せ。
Z ∞2exp(−t2 )dt
erfc(x) = √π x
7. 漸化式 Γ(x + 1) = xΓ(x) が成立することを示せ。また，n が正の整数であるとき，Γ(n + 1) = n!となることも示せ。

8. 漸化式に基づいて n = 20 までの Jn (x) と Yn (x) を求めよ。

参考図書
初等関数は殆ど全ての数値計算で必要となるものだけに，その研究の歴史は長い。古い文献は沢山あるが，ある程度，理論的な内容をまとめたものとしては初等関数の数値計算一松信
教育出版

1974 年が優れている。また，長年この研究に携わってきた著者が著した近似式のプログラミング
浜田穂積
培風館

1995 年も，各種近似式の係数が掲載されており役に立つ。また，

Software Manual for the Elementary Functions
W.J.Cody and Jr., W. Waite
Prentice-Hall
1980 年は，各初等関数のテスト Fotran77 プログラムも掲載されていて，有用である。
なお，初等関数も含め，応用上よく利用される特殊関数についての性質，数表をまとめた

第6章

70

初等関数の計算

Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables
M. Abramowitz and I.A.Stegun 編
Dover
1965 年は，大判で太い本であるが，手元に一冊あると便利である。初版は古いが，着実に版を重ねているので，入手は容易い。