§E20.17有理近似と初等関数の計算

最終更新

指数関数の値を浮動小数点演算で求めるには、四則演算を有限回行って値を計算することができる関数で指数関数を近似する。Taylor の定理が与える多項式1+x+x2/21+x+x^2/2はその一つであり、指数関数の00における0,1,20,1,2次の係数から定まる。同じ三つの係数からは、多項式の商である有理関数(2+x)/(2−x)(2+x)/(2-x)も定まる。x=1/2x=1/2におけるe1/2e^{1/2}との差は、前者では0.0237…0.0237\ldots、後者では−0.0179…-0.0179\ldotsであり、後者の値は加減算二回と除算一回で計算される。

しかし、この近似の誤差の評価は00の近くに限られる。(2+x)/(2−x)(2+x)/(2-x)の誤差は∣x∣≤1\lvert x\rvert\le1の範囲で∣x∣3\lvert x\rvert^3の定数倍で押さえられるが、この有理関数はx=2x=2に極をもつ。そこで、r=2−sxr=2^{-s}xに対してex=(er)2se^x=(e^r)^{2^s}が成り立つことを用い、縮小した引数rrで近似してからss回二乗して値を戻す。二乗は近似の相対誤差を拡大し、厳密に二乗しても、正の相対誤差は2s2^s倍以上になる。

00における係数から有理関数を定める一般の方法が Padé 近似であり、分子と分母の次数を指定し、分母とffの積から分子を引いた差が00の近くで∣x∣\lvert x\rvertの高い冪で押さえられることを要求する。分母の係数が満たすべき条件は一次方程式として書くことができ、Padé 近似は存在すれば有理関数としてただ一つであるが、存在しない場合もある。Padé 近似と引数の縮小と繰り返す二乗を組み合わせると、有限区間の浮動小数点数に対する指数関数の計算について、浮動小数点数系と区間に課した条件の下で、近似の誤差と丸めの誤差を合わせた相対誤差の上界を得ることができる。本記事では、Padé 近似の基本的な性質と、それを用いた指数関数の計算の精度について解説する。

1 Padé 近似

定義 1.1.m,n∈N≥0m,n\in\Nとし、I⊂RI\subset\Rを00を含む開区間、f ⁣:I→Rf\colon I\to\RをCm+n+1C^{m+n+1}級の関数とする。実係数の多項式の組(p,q)(p,q)が次の三条件を満たすとき、有理関数p/qp/qをffの[m/n][m/n]型の Padé 近似 (Padé approximant) という。

  1. ppの次数はmm以下であり、qqの次数はnn以下である。
  2. q(0)=1q(0)=1である。
  3. あるδ>0\delta>0とM≥0M\ge0が存在して、(−δ,δ)⊂I(-\delta,\delta)\subset Iであり、∣x∣<δ\lvert x\rvert<\deltaを満たす任意の実数xxに対して ∣q(x)f(x)−p(x)∣≤M∣x∣m+n+1\lvert q(x)f(x)-p(x)\rvert\le M\lvert x\rvert^{m+n+1} が成り立つ。

補題 1.2.N∈N≥0N\in\Nとし、g(x)=∑i=0dgixig(x)=\sum_{i=0}^dg_ix^iを実係数の多項式とし、i>di>dに対してgi:=0g_i:=0と置く。あるδ>0\delta>0とM≥0M\ge0が存在して∣x∣<δ\lvert x\rvert<\deltaを満たす任意の実数xxに対して∣g(x)∣≤M∣x∣N\lvert g(x)\rvert\le M\lvert x\rvert^Nが成り立つことと、0≤i<N0\le i<Nを満たす任意の整数iiに対してgi=0g_i=0であることは同値である。

証明.0≤i<N0\le i<Nを満たす任意のiiに対してgi=0g_i=0であるとする。このときg(x)=xN∑i=Ndgixi−Ng(x)=x^N\sum_{i=N}^dg_ix^{i-N}であり、M:=∑i=Nd∣gi∣M:=\sum_{i=N}^d\lvert g_i\rvert、δ:=1\delta:=1と置けば、∣x∣<1\lvert x\rvert<1を満たす任意のxxに対して∣g(x)∣≤M∣x∣N\lvert g(x)\rvert\le M\lvert x\rvert^Nが成り立つ。

逆に、δ>0\delta>0とM≥0M\ge0が不等式を満たし、gi≠0g_i\ne0を満たす整数0≤i<N0\le i<Nが存在すると仮定する。そのようなiiの最小値をllとし、h(x):=∑i=l+1dgixi−l−1h(x):=\sum_{i=l+1}^dg_ix^{i-l-1}(l=dl=dならばh=0h=0)と置くと、g(x)=xl(gl+xh(x))g(x)=x^l\bigl(g_l+xh(x)\bigr)である。0<∣x∣<δ0<\lvert x\rvert<\deltaならば

∣gl+xh(x)∣=∣g(x)∣∣x∣l≤M∣x∣N−l\lvert g_l+xh(x)\rvert=\frac{\lvert g(x)\rvert}{\lvert x\rvert^l}\le M\lvert x\rvert^{N-l}

であり、N−l≥1N-l\ge1であるから右辺はx→0x\to0で00に収束する。多項式hhは連続であるから左辺はx→0x\to0で∣gl∣\lvert g_l\rvertに収束するので、∣gl∣≤0\lvert g_l\rvert\le0が成り立つ。∣gl∣≤0\lvert g_l\rvert\le0とgl≠0g_l\ne0は両立しない。▨

補題 1.3.m,n∈N≥0m,n\in\Nとし、I⊂RI\subset\Rを00を含む開区間、f ⁣:I→Rf\colon I\to\RをCm+n+1C^{m+n+1}級の関数とする。0≤i≤m+n0\le i\le m+nに対してci:=f(i)(0)/i!c_i:=f^{(i)}(0)/i!と置き、負の整数iiに対してci:=0c_i:=0と置く。

  1. p(x)=∑i=0maixip(x)=\sum_{i=0}^ma_ix^iとq(x)=∑j=0nbjxjq(x)=\sum_{j=0}^nb_jx^jをb0=1b_0=1を満たす実係数の多項式とする。(p,q)(p,q)が定義 1.1 条件 (c)を満たすことと、次の二つの等式系が成り立つことは同値である。 ai=∑j=0min⁡{i,n}bjci−j(0≤i≤m),∑j=1nci−jbj=−ci(m+1≤i≤m+n).a_i=\sum_{j=0}^{\min\{i,n\}}b_jc_{i-j}\quad(0\le i\le m),\qquad\sum_{j=1}^nc_{i-j}b_j=-c_i\quad(m+1\le i\le m+n).
  2. n=0n=0であるか、n≥1n\ge1かつnn次正方行列Cm,n:=(cm+i−j)1≤i,j≤nC_{m,n}:=(c_{m+i-j})_{1\le i,j\le n}が正則であるとする。このとき定義 1.1の三条件を満たす組(p,q)(p,q)はただ一つ存在する。

証明.(1)を示す。N:=m+n+1N:=m+n+1、T(x):=∑i=0N−1cixiT(x):=\sum_{i=0}^{N-1}c_ix^iと置き、[−δ0,δ0]⊂I[-\delta_0,\delta_0]\subset Iを満たすδ0>0\delta_0>0をとる。∣x∣≤δ0\lvert x\rvert\le\delta_0を満たす各xxに対して、§D1.16 定理 2.1をa=0a=0とn=Nn=Nに適用すると、∣ξ∣≤∣x∣\lvert\xi\rvert\le\lvert x\rvertを満たす実数ξ\xiが存在してf(x)−T(x)=f(N)(ξ)xN/N!f(x)-T(x)=f^{(N)}(\xi)x^N/N!が成り立つ。f(N)f^{(N)}とqqは[−δ0,δ0][-\delta_0,\delta_0]で連続であるから、L:=max⁡∣t∣≤δ0∣f(N)(t)∣L:=\max_{\lvert t\rvert\le\delta_0}\lvert f^{(N)}(t)\rvertとQ:=max⁡∣t∣≤δ0∣q(t)∣Q:=\max_{\lvert t\rvert\le\delta_0}\lvert q(t)\rvertが存在し、∣x∣≤δ0\lvert x\rvert\le\delta_0ならば

∣q(x)f(x)−p(x)−(q(x)T(x)−p(x))∣=∣q(x)∣ ∣f(x)−T(x)∣≤QLN!∣x∣N\bigl\lvert q(x)f(x)-p(x)-\bigl(q(x)T(x)-p(x)\bigr)\bigr\rvert=\lvert q(x)\rvert\,\lvert f(x)-T(x)\rvert\le\frac{QL}{N!}\lvert x\rvert^N

である。三角不等式により、(p,q)(p,q)が定義 1.1 条件 (c)を満たすことと、あるδ>0\delta>0とM≥0M\ge0が存在して∣x∣<δ\lvert x\rvert<\deltaならば∣q(x)T(x)−p(x)∣≤M∣x∣N\lvert q(x)T(x)-p(x)\rvert\le M\lvert x\rvert^Nが成り立つことは同値である。qT−pqT-pは多項式であるから、補題 1.2により、後者はqT−pqT-pの0≤i≤N−10\le i\le N-1次の係数がすべて00であることと同値である。0≤i≤N−10\le i\le N-1に対してqTqTのii次の係数は∑j=0min⁡{i,n}bjci−j\sum_{j=0}^{\min\{i,n\}}b_jc_{i-j}であり、ppのii次の係数はi≤mi\le mならばaia_i、m<i≤N−1m<i\le N-1ならば00である。m+1≤i≤m+nm+1\le i\le m+nに対しては、b0=1b_0=1であり、j>ij>iならばci−j=0c_{i-j}=0であるから、∑j=0min⁡{i,n}bjci−j=ci+∑j=1nci−jbj\sum_{j=0}^{\min\{i,n\}}b_jc_{i-j}=c_i+\sum_{j=1}^nc_{i-j}b_jである。したがって係数が00である条件は二つの等式系に一致する。

(2)を示す。二つ目の等式系は、iiをm+im+iに置き換えると∑j=1ncm+i−jbj=−cm+i\sum_{j=1}^nc_{m+i-j}b_j=-c_{m+i}(1≤i≤n1\le i\le n)、すなわちCm,n(b1,…,bn)T=−(cm+1,…,cm+n)TC_{m,n}(b_1,\dots,b_n)^{\mathsf T}=-(c_{m+1},\dots,c_{m+n})^{\mathsf T}である。n=0n=0ならばこの等式系は空であり、Cm,nC_{m,n}が正則ならば解(b1,…,bn)(b_1,\dots,b_n)はただ一つである。いずれの場合も、一つ目の等式系はa0,…,ama_0,\dots,a_mをただ一通りに定める。(1)により、定義 1.1 条件 (a)と定義 1.1 条件 (b)を満たす組のうち定義 1.1 条件 (c)を満たすものは、このようにして定まる組だけである。▨

定理 1.4.m,n∈N≥0m,n\in\Nとし、I⊂RI\subset\Rを00を含む開区間、f ⁣:I→Rf\colon I\to\RをCm+n+1C^{m+n+1}級の関数とする。(p1,q1)(p_1,q_1)と(p2,q2)(p_2,q_2)がどちらも定義 1.1の三条件を満たすならば、多項式としてp1q2=p2q1p_1q_2=p_2q_1が成り立つ。したがってp1/q1=p2/q2p_1/q_1=p_2/q_2が有理関数体R(x)\R(x)の元として成り立ち、q1(x)q2(x)≠0q_1(x)q_2(x)\ne0を満たす任意の実数xxに対してp1(x)/q1(x)=p2(x)/q2(x)p_1(x)/q_1(x)=p_2(x)/q_2(x)が成り立つ。ffの[m/n][m/n]型の Padé 近似が存在すれば、それはR(x)\R(x)の元としてただ一つである。

証明.k=1,2k=1,2について、定義 1.1 条件 (c)を満たすδk>0\delta_k>0とMk≥0M_k\ge0をとる。δ:=min⁡{δ1,δ2,1}\delta:=\min\{\delta_1,\delta_2,1\}、Qk:=max⁡∣t∣≤1∣qk(t)∣Q_k:=\max_{\lvert t\rvert\le1}\lvert q_k(t)\rvertと置く。g:=p1q2−p2q1g:=p_1q_2-p_2q_1はIIの上で

g(x)=q1(x)(q2(x)f(x)−p2(x))−q2(x)(q1(x)f(x)−p1(x))g(x)=q_1(x)\bigl(q_2(x)f(x)-p_2(x)\bigr)-q_2(x)\bigl(q_1(x)f(x)-p_1(x)\bigr)

を満たすので、∣x∣<δ\lvert x\rvert<\deltaならば∣g(x)∣≤(Q1M2+Q2M1)∣x∣m+n+1\lvert g(x)\rvert\le(Q_1M_2+Q_2M_1)\lvert x\rvert^{m+n+1}である。定義 1.1 条件 (a)によりggの次数はm+nm+n以下であるから、補題 1.2をN=m+n+1N=m+n+1として適用するとggのすべての係数は00であり、p1q2=p2q1p_1q_2=p_2q_1が成り立つ。q1(0)=q2(0)=1q_1(0)=q_2(0)=1であるからq1q_1とq2q_2は00でない多項式であり、両辺をq1q2q_1q_2で割ってR(x)\R(x)における等式と、q1(x)q2(x)≠0q_1(x)q_2(x)\ne0を満たすxxにおける値の等式を得る。▨

例 1.5.補題 1.3の記号を用いる。

  1. f(x)=1+x2f(x)=1+x^2、m=n=1m=n=1とする。c0=1c_0=1、c1=0c_1=0、c2=1c_2=1であり、補題 1.3 (1)の二つ目の等式系はc1b1=−c2c_1b_1=-c_2、すなわち0⋅b1=−10\cdot b_1=-1である。この等式を満たすb1b_1は存在しないので、ffの[1/1][1/1]型の Padé 近似は存在しない。定義 1.1 条件 (b)を外すと、p(x)=xp(x)=x、q(x)=xq(x)=xは定義 1.1 条件 (a)とq(x)f(x)−p(x)=x3q(x)f(x)-p(x)=x^3を満たす。この組のp/q=1p/q=1について1−f(x)=−x21-f(x)=-x^2の22次の係数は00でないので、補題 1.2により1−f(x)1-f(x)は∣x∣3\lvert x\rvert^3の定数倍で押さえられない。
  2. f(x)=1f(x)=1、m=n=1m=n=1とする。c0=1c_0=1、c1=c2=0c_1=c_2=0であり、C1,1=(c1)=(0)C_{1,1}=(c_1)=(0)は正則でない。二つ目の等式系0⋅b1=00\cdot b_1=0は任意のb1b_1で成り立ち、一つ目の等式系はa0=1a_0=1、a1=b1a_1=b_1を与える。したがって任意のb∈Rb\in\Rに対して組(1+bx,1+bx)(1+bx,1+bx)は三条件を満たし、組は一意でない。有理関数(1+bx)/(1+bx)(1+bx)/(1+bx)はbbによらず11である。

注意 1.6.Cm,nC_{m,n}が正則ならば、§E20.5 定理 2.5 (1)によりCm,nC_{m,n}の部分ピボット付き消去は厳密算術でどの段でも失敗せず、§E20.5 定理 2.5 (2)により前進代入と後退代入が補題 1.3 (1)の二つ目の等式系のただ一つの解(b1,…,bn)(b_1,\dots,b_n)を与える。

2 指数関数の[1/1][1/1]近似

命題 2.1.x≠2x\ne2に対してr1,1(x):=(1+x/2)/(1−x/2)=(2+x)/(2−x)r_{1,1}(x):=(1+x/2)/(1-x/2)=(2+x)/(2-x)と置く。0<ρ≤10<\rho\le1に対して

K(ρ):=(1+ρ)eρ6(2−ρ)K(\rho):=\frac{(1+\rho)e^\rho}{6(2-\rho)}

と置く。

  1. 組(1+x/2, 1−x/2)(1+x/2,\,1-x/2)は、f=exp⁡f=\exp、m=n=1m=n=1について定義 1.1の三条件を満たすただ一つの組である。特にr1,1r_{1,1}は指数関数の[1/1][1/1]型の Padé 近似である。
  2. 0<ρ≤10<\rho\le1かつ∣x∣≤ρ\lvert x\rvert\le\rhoならば、1−x/2≥1−ρ/2≥1/21-x/2\ge1-\rho/2\ge1/2である。
  3. 0<ρ≤10<\rho\le1かつ∣x∣≤ρ\lvert x\rvert\le\rhoならば、∣ξ∣≤∣x∣\lvert\xi\rvert\le\lvert x\rvertと∣ξ′∣≤∣x∣\lvert\xi'\rvert\le\lvert x\rvertを満たす実数ξ,ξ′\xi,\xi'が存在して ex−r1,1(x)=−(1+ξ)eξ12(1−x/2)x3,r1,1(x)e−x−1=(1+ξ′)eξ′12(1−x/2)x3e^x-r_{1,1}(x)=-\frac{(1+\xi)e^{\xi}}{12(1-x/2)}x^3,\qquad r_{1,1}(x)e^{-x}-1=\frac{(1+\xi')e^{\xi'}}{12(1-x/2)}x^3 が成り立ち、 ∣ex−r1,1(x)∣≤K(ρ)∣x∣3,∣r1,1(x)e−x−1∣≤K(ρ)∣x∣3\lvert e^x-r_{1,1}(x)\rvert\le K(\rho)\lvert x\rvert^3,\qquad\lvert r_{1,1}(x)e^{-x}-1\rvert\le K(\rho)\lvert x\rvert^3 である。特にK(1/2)=e1/2/6<0.275K(1/2)=e^{1/2}/6<0.275である。

証明.(1)を示す。f=exp⁡f=\expについてc0=c1=1c_0=c_1=1、c2=1/2c_2=1/2であり、C1,1=(c1)=(1)C_{1,1}=(c_1)=(1)は正則である。補題 1.3 (1)の二つ目の等式系c1b1=−c2c_1b_1=-c_2からb1=−1/2b_1=-1/2であり、一つ目の等式系からa0=c0=1a_0=c_0=1、a1=c1+b1c0=1/2a_1=c_1+b_1c_0=1/2である。補題 1.3 (2)により、三条件を満たす組は(1+x/2, 1−x/2)(1+x/2,\,1-x/2)だけである。

∣x∣≤ρ≤1\lvert x\rvert\le\rho\le1ならばx/2≤ρ/2≤1/2x/2\le\rho/2\le1/2であるから、(2)が成り立つ。

(3)を示す。実数yyに対してg(y):=(1−y/2)ey−(1+y/2)g(y):=(1-y/2)e^y-(1+y/2)と置く。

g′(y)=1−y2ey−12,g′′(y)=−y2ey,g′′′(y)=−1+y2eyg'(y)=\frac{1-y}2e^y-\frac12,\qquad g''(y)=-\frac y2e^y,\qquad g'''(y)=-\frac{1+y}2e^y

であり、g(0)=g′(0)=g′′(0)=0g(0)=g'(0)=g''(0)=0である。§D1.16 定理 2.1をa=0a=0とn=3n=3に適用すると、各yyに対して∣η∣≤∣y∣\lvert\eta\rvert\le\lvert y\rvertを満たす実数η\etaが存在してg(y)=g′′′(η)y3/6=−(1+η)eηy3/12g(y)=g'''(\eta)y^3/6=-(1+\eta)e^\eta y^3/12が成り立つ。(2)により1−x/2>01-x/2>0であり、

ex−r1,1(x)=g(x)1−x/2,r1,1(x)e−x−1=(1+x/2)e−x−(1−x/2)1−x/2=g(−x)1−x/2e^x-r_{1,1}(x)=\frac{g(x)}{1-x/2},\qquad r_{1,1}(x)e^{-x}-1=\frac{(1+x/2)e^{-x}-(1-x/2)}{1-x/2}=\frac{g(-x)}{1-x/2}

である。y=xy=xに対するη\etaをξ\xi、y=−xy=-xに対するη\etaをξ′\xi'とすると、g(−x)=(1+ξ′)eξ′x3/12g(-x)=(1+\xi')e^{\xi'}x^3/12であるから、二つの等式が成り立つ。関数t↦(1+t)ett\mapsto(1+t)e^tの導関数は(2+t)et(2+t)e^tであり、[−1,∞)[-1,\infty)で正であるから、この関数は[−1,∞)[-1,\infty)で増加し、−1≤t≤ρ-1\le t\le\rhoならば0≤(1+t)et≤(1+ρ)eρ0\le(1+t)e^t\le(1+\rho)e^\rhoである。∣ξ∣,∣ξ′∣≤ρ≤1\lvert\xi\rvert,\lvert\xi'\rvert\le\rho\le1と(2)により、二つの誤差の絶対値はどちらも

(1+ρ)eρ12(1−ρ/2)∣x∣3=K(ρ)∣x∣3\frac{(1+\rho)e^\rho}{12(1-\rho/2)}\lvert x\rvert^3=K(\rho)\lvert x\rvert^3

以下である。K(1/2)=(3/2)e1/2/(6⋅3/2)=e1/2/6K(1/2)=(3/2)e^{1/2}/(6\cdot3/2)=e^{1/2}/6であり、e1/2<1.65e^{1/2}<1.65からK(1/2)<0.275K(1/2)<0.275である。▨

例 2.2.r1,1(x)=(2+x)/(2−x)r_{1,1}(x)=(2+x)/(2-x)とT2(x):=1+x+x2/2T_2(x):=1+x+x^2/2はどちらも指数関数の係数ci=1/i!c_i=1/i!(i=0,1,2i=0,1,2)から定まり、補題 1.3 (2)によりT2T_2は指数関数の[2/0][2/0]型の Padé 近似である。r1,1(x)r_{1,1}(x)の評価は加減算二回と除算一回を、T2(x)=1+x(1+x⋅(1/2))T_2(x)=1+x(1+x\cdot(1/2))の評価は乗算二回と加算二回を行う。§D1.16 定理 2.1によりex−T2(x)=eξx3/6e^x-T_2(x)=e^{\xi}x^3/6(∣ξ∣≤∣x∣\lvert\xi\rvert\le\lvert x\rvert)であるから、x→0x\to0のときx3x^3の係数は、ex−T2(x)e^x-T_2(x)では1/61/6に、命題 2.1 (3)の表示のex−r1,1(x)e^x-r_{1,1}(x)では−1/12-1/12に収束する。x=1/2x=1/2ではe1/2=1.6487212…e^{1/2}=1.6487212\ldots、r1,1(1/2)=5/3r_{1,1}(1/2)=5/3、T2(1/2)=1.625T_2(1/2)=1.625であり、e1/2−r1,1(1/2)=−0.017945…e^{1/2}-r_{1,1}(1/2)=-0.017945\ldots、e1/2−T2(1/2)=0.023721…e^{1/2}-T_2(1/2)=0.023721\ldotsである。前者の絶対値は命題 2.1 (3)の上界K(1/2)/8=0.034348…K(1/2)/8=0.034348\ldots以下である。x=−1/2x=-1/2ではe−1/2=0.6065306…e^{-1/2}=0.6065306\ldots、r1,1(−1/2)=3/5r_{1,1}(-1/2)=3/5、T2(−1/2)=0.625T_2(-1/2)=0.625であり、e−1/2−r1,1(−1/2)=0.0065306…e^{-1/2}-r_{1,1}(-1/2)=0.0065306\ldots、e−1/2−T2(−1/2)=−0.018469…e^{-1/2}-T_2(-1/2)=-0.018469\ldotsである。r1,1r_{1,1}はx=2x=2に極をもつ。x→2−0x\to2-0のときr1,1(x)→+∞r_{1,1}(x)\to+\inftyであり、x>2x>2ならばr1,1(x)<0<exr_{1,1}(x)<0<e^xである。x=1.9x=1.9ではr1,1(1.9)=39r_{1,1}(1.9)=39、e1.9=6.6858944…e^{1.9}=6.6858944\ldots、T2(1.9)=4.705T_2(1.9)=4.705である。

3 引数の縮小と繰り返す二乗

定義 3.1.x∈Rx\in\R、s∈N≥0s\in\Nとし、r:=2−sxr:=2^{-s}xと置く。このときex=(er)2se^x=(e^r)^{2^s}が成り立つ。xxをrrに置き換えることを 引数の縮小 (argument reduction) といい、ssをその 縮小回数 (number of halvings) という。FFを浮動小数点数系、fl⁡\operatorname{fl}をFFの最近接丸めとし、y0∈Fy_0\in Fとする。0≤j<s0\le j<sについて順にyj+1:=fl⁡(yj⋅yj)y_{j+1}:=\operatorname{fl}(y_j\cdot y_j)と置いてy1,…,ysy_1,\dots,y_sを定めることを、y0y_0からのss回の 繰り返す二乗 (repeated squaring) という。yj+1y_{j+1}は、厳密な積yj⋅yjy_j\cdot y_jがNmax⁡N_{\max}以下であるときに定まる。

補題 3.2.k∈N≥1k\in\NNとし、t1,…,tk∈[0,1]t_1,\dots,t_k\in[0,1]とする。このとき

∏i=1k(1−ti)≥1−∑i=1kti\prod_{i=1}^k(1-t_i)\ge1-\sum_{i=1}^kt_i

が成り立つ。

証明. 演習とする(問題 5.1)。▨

定理 3.3.FFを浮動小数点数系、uuをその単位丸め誤差、fl⁡\operatorname{fl}をFFの最近接丸めとする。r∈Rr\in\R、s∈N≥0s\in\N、y0∈Fy_0\in Fとし、y1,…,ysy_1,\dots,y_sをy0y_0からのss回の繰り返す二乗とする。0≤j<s0\le j<sを満たす各jjについて、厳密な積yj⋅yjy_j\cdot y_jは正規範囲にあるとする。0≤j≤s0\le j\le sに対して実数ηj\eta_jをyj=e2jr(1+ηj)y_j=e^{2^jr}(1+\eta_j)により定め、ε0:=η0\varepsilon_0:=\eta_0と置く。

  1. z0:=y0z_0:=y_0、zj+1:=zj2z_{j+1}:=z_j^2(0≤j<s0\le j<s)と置くと、zs=e2sr(1+ε0)2sz_s=e^{2^sr}(1+\varepsilon_0)^{2^s}が成り立つ。すなわち、厳密な二乗によるzsz_sのe2sre^{2^sr}に対する相対誤差は(1+ε0)2s−1(1+\varepsilon_0)^{2^s}-1である。
  2. 0≤j<s0\le j<sを満たす各jjに対して、∣δj∣≤u\lvert\delta_j\rvert\le uを満たす実数δj\delta_jが存在して 1+ηj+1=(1+ηj)2(1+δj)1+\eta_{j+1}=(1+\eta_j)^2(1+\delta_j) が成り立つ。したがって 1+ηs=(1+ε0)2s∏j=0s−1(1+δj)2s−1−j1+\eta_s=(1+\varepsilon_0)^{2^s}\prod_{j=0}^{s-1}(1+\delta_j)^{2^{s-1-j}} である。
  3. 0≤b<10\le b<1が1−b≤1+ε0≤(1−b)−11-b\le1+\varepsilon_0\le(1-b)^{-1}を満たすとし、0≤j≤s0\le j\le sに対してσj:=2jb+(2j−1)u\sigma_j:=2^jb+(2^j-1)uと置き、σs<1\sigma_s<1とする。このとき0≤j≤s0\le j\le sを満たす各jjに対して 1−σj≤1+ηj≤(1−σj)−1,∣ηj∣≤σj1−σj1-\sigma_j\le1+\eta_j\le(1-\sigma_j)^{-1},\qquad\lvert\eta_j\rvert\le\frac{\sigma_j}{1-\sigma_j} が成り立つ。∣ε0∣≤b\lvert\varepsilon_0\rvert\le bならば1−b≤1+ε0≤(1−b)−11-b\le1+\varepsilon_0\le(1-b)^{-1}が成り立つ。

証明.(1)を示す。j=0j=0ではz0=y0=er(1+ε0)z_0=y_0=e^r(1+\varepsilon_0)である。zj=e2jr(1+ε0)2jz_j=e^{2^jr}(1+\varepsilon_0)^{2^j}ならばzj+1=zj2=e2j+1r(1+ε0)2j+1z_{j+1}=z_j^2=e^{2^{j+1}r}(1+\varepsilon_0)^{2^{j+1}}であり、jjに関する帰納法により主張が成り立つ。

(2)を示す。yj∈Fy_j\in Fであり、厳密な積yj⋅yjy_j\cdot y_jは正規範囲にあるから、§E20.1 系 3.2 (2)により∣δj∣≤u\lvert\delta_j\rvert\le uを満たすδj\delta_jが存在してyj+1=yj2(1+δj)y_{j+1}=y_j^2(1+\delta_j)である。yj=e2jr(1+ηj)y_j=e^{2^jr}(1+\eta_j)を代入するとe2j+1r(1+ηj+1)=e2j+1r(1+ηj)2(1+δj)e^{2^{j+1}r}(1+\eta_{j+1})=e^{2^{j+1}r}(1+\eta_j)^2(1+\delta_j)であり、両辺をe2j+1r>0e^{2^{j+1}r}>0で割って漸化式を得る。1+ηj=(1+ε0)2j∏i=0j−1(1+δi)2j−1−i1+\eta_j=(1+\varepsilon_0)^{2^j}\prod_{i=0}^{j-1}(1+\delta_i)^{2^{j-1-i}}ならば、漸化式により

1+ηj+1=(1+ε0)2j+1∏i=0j−1(1+δi)2j−i⋅(1+δj)=(1+ε0)2j+1∏i=0j(1+δi)2j−i1+\eta_{j+1}=(1+\varepsilon_0)^{2^{j+1}}\prod_{i=0}^{j-1}(1+\delta_i)^{2^{j-i}}\cdot(1+\delta_j)=(1+\varepsilon_0)^{2^{j+1}}\prod_{i=0}^{j}(1+\delta_i)^{2^{j-i}}

であり、j=0j=0の場合は空積の表示として成り立つので、jjに関する帰納法によりj=sj=sの表示を得る。

(3)を示す。σ0=b\sigma_0=bであり、σj+1=2σj+u\sigma_{j+1}=2\sigma_j+uであるから、0≤σ0≤σ1≤⋯≤σs<10\le\sigma_0\le\sigma_1\le\dots\le\sigma_s<1である。u=β1−p/2≤1/2u=\beta^{1-p}/2\le1/2であり、∣δj∣≤u\lvert\delta_j\rvert\le uから1−u≤1+δj≤1+u≤(1−u)−11-u\le1+\delta_j\le1+u\le(1-u)^{-1}である。最後の不等式は(1+u)(1−u)=1−u2≤1(1+u)(1-u)=1-u^2\le1による。j=0j=0では仮定により1−σ0≤1+η0≤(1−σ0)−11-\sigma_0\le1+\eta_0\le(1-\sigma_0)^{-1}である。j<sj<sとし、1−σj≤1+ηj≤(1−σj)−11-\sigma_j\le1+\eta_j\le(1-\sigma_j)^{-1}が成り立つとする。1−σj>01-\sigma_j>0であるから(1−σj)2≤(1+ηj)2≤(1−σj)−2(1-\sigma_j)^2\le(1+\eta_j)^2\le(1-\sigma_j)^{-2}であり、(2)の漸化式により

(1−σj)2(1−u)≤1+ηj+1≤((1−σj)2(1−u))−1(1-\sigma_j)^2(1-u)\le1+\eta_{j+1}\le\bigl((1-\sigma_j)^2(1-u)\bigr)^{-1}

である。補題 3.2をt1=t2=σjt_1=t_2=\sigma_j、t3=ut_3=uに適用すると(1−σj)2(1−u)≥1−2σj−u=1−σj+1>0(1-\sigma_j)^2(1-u)\ge1-2\sigma_j-u=1-\sigma_{j+1}>0であるから、1−σj+1≤1+ηj+1≤(1−σj+1)−11-\sigma_{j+1}\le1+\eta_{j+1}\le(1-\sigma_{j+1})^{-1}が成り立つ。jjに関する帰納法により、0≤j≤s0\le j\le sを満たす各jjで一つ目の不等式が成り立つ。このときηj≤(1−σj)−1−1=σj/(1−σj)\eta_j\le(1-\sigma_j)^{-1}-1=\sigma_j/(1-\sigma_j)であり、−ηj≤σj≤σj/(1−σj)-\eta_j\le\sigma_j\le\sigma_j/(1-\sigma_j)である。∣ε0∣≤b\lvert\varepsilon_0\rvert\le bならば1−b≤1+ε0≤1+b1-b\le1+\varepsilon_0\le1+bであり、(1+b)(1−b)=1−b2≤1(1+b)(1-b)=1-b^2\le1から1+b≤(1−b)−11+b\le(1-b)^{-1}である。▨

注意 3.4.定理 3.3 (1)においてε0>0\varepsilon_0>0ならば、(1+ε0)2s(1+\varepsilon_0)^{2^s}の二項展開の各項は非負であるから(1+ε0)2s−1≥2sε0(1+\varepsilon_0)^{2^s}-1\ge2^s\varepsilon_0である。x∈Rx\in\Rを固定し、r=2−sxr=2^{-s}xとする。ssによらない定数K>0K>0とc≥0c\ge0について∣ε0∣≤b:=K∣r∣3+cu\lvert\varepsilon_0\rvert\le b:=K\lvert r\rvert^3+cuであるとき、定理 3.3 (3)のσs\sigma_sは

σs=K∣x∣3 4−s+c 2su+(2s−1)u\sigma_s=K\lvert x\rvert^3\,4^{-s}+c\,2^su+(2^s-1)u

である。第1項は縮小回数ssについて減少し、第2項と第3項はssについて増加する。

4 有限区間での計算精度

補題 4.1.F=F(2,p,emin⁡,emax⁡)F=F(2,p,e_{\min},e_{\max})をp≥3p\ge3、emin⁡≤−2e_{\min}\le-2、emax⁡≥2e_{\max}\ge2を満たす浮動小数点数系とし、u=2−pu=2^{-p}をその単位丸め誤差、fl⁡\operatorname{fl}をFFの最近接丸めとする。r∈Fr\in Fが∣r∣≤1\lvert r\rvert\le1を満たすとする。このとき2∈F2\in Fであり、n^:=fl⁡(2+r)\hat n:=\operatorname{fl}(2+r)、d^:=fl⁡(2−r)\hat d:=\operatorname{fl}(2-r)、y^:=fl⁡(n^/d^)\hat y:=\operatorname{fl}(\hat n/\hat d)は定まり、厳密な商n^/d^\hat n/\hat dは[1/4,4][1/4,4]に属する。さらに、(1−u)3≤1+θ≤(1−u)−3(1-u)^3\le1+\theta\le(1-u)^{-3}を満たす実数θ\thetaが存在して

y^=2+r2−r(1+θ)\hat y=\frac{2+r}{2-r}(1+\theta)

が成り立つ。

証明.s1=22−ps_1=2^{2-p}であるから2=2p−1s12=2^{p-1}s_1であり、1≤emax⁡1\le e_{\max}から22は指数11の正規化数である。Nmax⁡=2emax⁡+1(1−2−p)≥8⋅7/8=7N_{\max}=2^{e_{\max}+1}(1-2^{-p})\ge8\cdot7/8=7であり、∣r∣≤1\lvert r\rvert\le1から1≤2±r≤3≤Nmax⁡1\le2\pm r\le3\le N_{\max}である。F=−FF=-Fであるから−r∈F-r\in Fであり、§E20.1 系 3.3により∣α1∣,∣α2∣≤u\lvert\alpha_1\rvert,\lvert\alpha_2\rvert\le uを満たすα1,α2\alpha_1,\alpha_2が存在してn^=(2+r)(1+α1)\hat n=(2+r)(1+\alpha_1)、d^=(2−r)(1+α2)\hat d=(2-r)(1+\alpha_2)である。u≤1/8u\le1/8であるからd^≥1−u>0\hat d\ge1-u>0であり、(1−u)/(1+u)≥7/9(1-u)/(1+u)\ge7/9から

14<727≤13⋅1−u1+u≤n^d^≤3⋅1+u1−u≤277<4\frac14<\frac{7}{27}\le\frac13\cdot\frac{1-u}{1+u}\le\frac{\hat n}{\hat d}\le3\cdot\frac{1+u}{1-u}\le\frac{27}{7}<4

である。2emin⁡≤1/42^{e_{\min}}\le1/4かつ4<Nmax⁡4<N_{\max}であるからn^/d^\hat n/\hat dは正規範囲にあり、§E20.1 系 3.2 (2)により∣α3∣≤u\lvert\alpha_3\rvert\le uを満たすα3\alpha_3が存在してy^=(n^/d^)(1+α3)\hat y=(\hat n/\hat d)(1+\alpha_3)である。したがって

y^=2+r2−r⋅(1+α1)(1+α3)1+α2\hat y=\frac{2+r}{2-r}\cdot\frac{(1+\alpha_1)(1+\alpha_3)}{1+\alpha_2}

である。i=1,2,3i=1,2,3について1−u≤1+αi≤1+u≤(1−u)−11-u\le1+\alpha_i\le1+u\le(1-u)^{-1}であり、(1+u)−1≤(1+α2)−1≤(1−u)−1(1+u)^{-1}\le(1+\alpha_2)^{-1}\le(1-u)^{-1}と(1+u)−1≥1−u(1+u)^{-1}\ge1-uにより1−u≤(1+α2)−1≤(1−u)−11-u\le(1+\alpha_2)^{-1}\le(1-u)^{-1}である。1+θ:=(1+α1)(1+α3)(1+α2)−11+\theta:=(1+\alpha_1)(1+\alpha_3)(1+\alpha_2)^{-1}と置けば(1−u)3≤1+θ≤(1−u)−3(1-u)^3\le1+\theta\le(1-u)^{-3}である。▨

定義 4.2.F=F(2,p,emin⁡,emax⁡)F=F(2,p,e_{\min},e_{\max})をp≥3p\ge3、emin⁡≤−2e_{\min}\le-2、emax⁡≥2e_{\max}\ge2を満たす浮動小数点数系、fl⁡\operatorname{fl}をFFの最近接丸めとし、k∈N≥0k\in\Nはk≤−emin⁡−1k\le-e_{\min}-1を満たすとする。x∈Fx\in Fに対して、∣x∣≤2s−k\lvert x\rvert\le2^{s-k}を満たす最小のs∈N≥0s\in\Nをs(x)s(x)と置き、r(x):=2−s(x)xr(x):=2^{-s(x)}xと置く。r(x)∈Fr(x)\in Fのとき、y^0:=fl⁡(fl⁡(2+r(x))/fl⁡(2−r(x)))\hat y_0:=\operatorname{fl}\bigl(\operatorname{fl}(2+r(x))/\operatorname{fl}(2-r(x))\bigr)と置き、y^1,…,y^s(x)\hat y_1,\dots,\hat y_{s(x)}をy^0\hat y_0からのs(x)s(x)回の繰り返す二乗とする。これらの値がすべて定まるとき、Ek(x):=y^s(x)E_k(x):=\hat y_{s(x)}と置く。

定理 4.3.F=F(2,p,emin⁡,emax⁡)F=F(2,p,e_{\min},e_{\max})、fl⁡\operatorname{fl}、kkを定義 4.2のとおりとし、u=2−pu=2^{-p}をFFの単位丸め誤差とする。

ρ:=2−k,K:=(1+ρ)eρ6(2−ρ)\rho:=2^{-k},\qquad K:=\frac{(1+\rho)e^\rho}{6(2-\rho)}

と置く。X>0X>0とし、X≤2S−kX\le2^{S-k}を満たす最小のS∈N≥0S\in\Nをとって

σ:=2SKρ3+(2S+2−1)u\sigma:=2^SK\rho^3+(2^{S+2}-1)u

と置く。次の条件を仮定する。

  1. σ<1\sigma<1である。
  2. 2emin⁡≤e−X(1−σ)22^{e_{\min}}\le e^{-X}(1-\sigma)^2かつeX(1−σ)−2≤Nmax⁡e^X(1-\sigma)^{-2}\le N_{\max}である。

このとき、∣x∣≤X\lvert x\rvert\le Xを満たす任意のx∈Fx\in Fに対してEk(x)E_k(x)は定まり、

1−σ≤Ek(x)e−x≤(1−σ)−1,∣Ek(x)−ex∣ex≤σ1−σ1-\sigma\le E_k(x)e^{-x}\le(1-\sigma)^{-1},\qquad\frac{\lvert E_k(x)-e^x\rvert}{e^x}\le\frac{\sigma}{1-\sigma}

が成り立つ。

証明.x∈Fx\in F、∣x∣≤X\lvert x\rvert\le Xをとり、s:=s(x)s:=s(x)、r:=r(x)=2−sxr:=r(x)=2^{-s}xと置く。∣x∣≤X≤2S−k\lvert x\rvert\le X\le2^{S-k}であるから、s(x)s(x)の最小性によりs≤Ss\le Sである。s(x)s(x)の定義により∣r∣=2−s∣x∣≤2−k=ρ≤1\lvert r\rvert=2^{-s}\lvert x\rvert\le2^{-k}=\rho\le1である。

s=0s=0ならばr=x∈Fr=x\in Fである。s≥1s\ge1ならば、s(x)s(x)の最小性により∣x∣>2s−1−k\lvert x\rvert>2^{s-1-k}であるから、∣r∣>2−k−1≥2emin⁡\lvert r\rvert>2^{-k-1}\ge2^{e_{\min}}である。このとき∣x∣≥∣r∣≥2emin⁡\lvert x\rvert\ge\lvert r\rvert\ge2^{e_{\min}}であり、ℓ:=e(∣x∣)\ell:=e(\lvert x\rvert)を§E20.1 補題 1.2の整数とすると2ℓ≤∣x∣<2ℓ+12^\ell\le\lvert x\rvert<2^{\ell+1}である。§E20.1 補題 1.3 (2)により1≤∣m∣≤2p−11\le\lvert m\rvert\le2^p-1を満たす整数mmが存在してx=msℓx=ms_\ellであり、r=msℓ−sr=ms_{\ell-s}である。2emin⁡<∣r∣<2ℓ+1−s2^{e_{\min}}<\lvert r\rvert<2^{\ell+1-s}からemin⁡≤ℓ−s≤emax⁡e_{\min}\le\ell-s\le e_{\max}であり、∣r∣≤1≤Nmax⁡\lvert r\rvert\le1\le N_{\max}であるから、§E20.1 補題 1.3 (1)によりr∈Fr\in Fである。

補題 4.1によりy^0\hat y_0は定まり、(1−u)3≤1+θ≤(1−u)−3(1-u)^3\le1+\theta\le(1-u)^{-3}を満たすθ\thetaが存在してy^0=r1,1(r)(1+θ)\hat y_0=r_{1,1}(r)(1+\theta)である。ここでr1,1(t):=(2+t)/(2−t)r_{1,1}(t):=(2+t)/(2-t)である。a:=∣r1,1(r)e−r−1∣a:=\lvert r_{1,1}(r)e^{-r}-1\rvertと置くと、命題 2.1 (3)によりa≤K∣r∣3≤Kρ3a\le K\lvert r\rvert^3\le K\rho^3である。2sa≤2SKρ3≤σ<12^sa\le2^SK\rho^3\le\sigma<1であるからa<1a<1であり、1−a≤r1,1(r)e−r≤1+a≤(1−a)−11-a\le r_{1,1}(r)e^{-r}\le1+a\le(1-a)^{-1}である。ε0:=y^0e−r−1\varepsilon_0:=\hat y_0e^{-r}-1と置くと

(1−a)(1−u)3≤1+ε0≤((1−a)(1−u)3)−1(1-a)(1-u)^3\le1+\varepsilon_0\le\bigl((1-a)(1-u)^3\bigr)^{-1}

であり、補題 3.2により(1−a)(1−u)3≥1−a−3u(1-a)(1-u)^3\ge1-a-3uである。b:=a+3ub:=a+3uと置き、0≤j≤s0\le j\le sに対して

σj:=2jb+(2j−1)u=2ja+(2j+2−1)u\sigma_j:=2^jb+(2^j-1)u=2^ja+(2^{j+2}-1)u

と置く。σj≤σs≤2SKρ3+(2S+2−1)u=σ<1\sigma_j\le\sigma_s\le2^SK\rho^3+(2^{S+2}-1)u=\sigma<1であり、特にb=σ0<1b=\sigma_0<1であるから1−b≤1+ε0≤(1−b)−11-b\le1+\varepsilon_0\le(1-b)^{-1}が成り立つ。

0≤j≤s0\le j\le sに対して、y^0,…,y^j\hat y_0,\dots,\hat y_jが定まり、0≤i<j0\le i<jについて厳密な積y^i⋅y^i\hat y_i\cdot\hat y_iが正規範囲にあり、1−σj≤y^je−2jr≤(1−σj)−11-\sigma_j\le\hat y_je^{-2^jr}\le(1-\sigma_j)^{-1}が成り立つことをP(j)P(j)と書く。P(0)P(0)は前段で示した。j<sj<sとし、P(j)P(j)が成り立つとする。∣2j+1r∣≤2s∣r∣=∣x∣≤X\lvert2^{j+1}r\rvert\le2^s\lvert r\rvert=\lvert x\rvert\le Xとσj≤σ\sigma_j\le\sigmaにより、

y^j⋅y^j=e2j+1r(y^je−2jr)2∈[e−X(1−σ)2,  eX(1−σ)−2]\hat y_j\cdot\hat y_j=e^{2^{j+1}r}\bigl(\hat y_je^{-2^jr}\bigr)^2\in\bigl[e^{-X}(1-\sigma)^2,\;e^X(1-\sigma)^{-2}\bigr]

であり、条件 (b)によりこの区間は正規範囲に含まれる。したがってy^j+1\hat y_{j+1}は定まる。y^0\hat y_0からのj+1j+1回の繰り返す二乗に定理 3.3 (3)をssをj+1j+1として適用すると、1−σj+1≤y^j+1e−2j+1r≤(1−σj+1)−11-\sigma_{j+1}\le\hat y_{j+1}e^{-2^{j+1}r}\le(1-\sigma_{j+1})^{-1}であり、P(j+1)P(j+1)が成り立つ。jjに関する帰納法によりP(s)P(s)が成り立つ。2sr=x2^sr=xであるから、P(s)P(s)によりEk(x)=y^sE_k(x)=\hat y_sは定まり、

1−σ≤1−σs≤Ek(x)e−x≤(1−σs)−1≤(1−σ)−11-\sigma\le1-\sigma_s\le E_k(x)e^{-x}\le(1-\sigma_s)^{-1}\le(1-\sigma)^{-1}

である。したがってEk(x)e−x−1≤σ/(1−σ)E_k(x)e^{-x}-1\le\sigma/(1-\sigma)かつ1−Ek(x)e−x≤σ≤σ/(1−σ)1-E_k(x)e^{-x}\le\sigma\le\sigma/(1-\sigma)であり、両辺にex>0e^x>0を掛けて相対誤差の不等式を得る。▨

例 4.4.F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)、u=2−53u=2^{-53}とし、fl⁡\operatorname{fl}を最近接偶数丸めとする。定理 4.3のσ\sigmaと上界σ/(1−σ)\sigma/(1-\sigma)を 60 桁の十進演算で計算し、CPython の float で実行したEk(x)E_k(x)の相対誤差(Ek(x)−ex)/ex(E_k(x)-e^x)/e^xと並べる。以下に挙げるkkのうちσ<1\sigma<1を満たすものについて、定理 4.3 条件 (b)も成り立つ。

  1. X=1X=1とするとS=kS=kであり、σ=K4−k+(2k+2−1)u\sigma=K4^{-k}+(2^{k+2}-1)uである。x=1x=1ではs(1)=ks(1)=k、r(1)=2−kr(1)=2^{-k}である。

    kk σ/(1−σ)\sigma/(1-\sigma) Ek(1)E_k(1)の相対誤差(観察)
    00 9.6499.649 1.036×10−11.036\times10^{-1}
    22 9.646×10−39.646\times10^{-3} 5.272×10−35.272\times10^{-3}
    44 3.802×10−43.802\times10^{-4} 3.258×10−43.258\times10^{-4}
    88 1.284×10−61.284\times10^{-6} 1.272×10−61.272\times10^{-6}
    1212 4.972×10−94.972\times10^{-9} 4.967×10−94.967\times10^{-9}
    1515 9.217×10−119.217\times10^{-11} 7.782×10−117.782\times10^{-11}
    1616 4.851×10−114.851\times10^{-11} 1.950×10−111.950\times10^{-11}
    1717 6.306×10−116.306\times10^{-11} 1.950×10−111.950\times10^{-11}
    2020 4.657×10−104.657\times10^{-10} −9.668×10−12-9.668\times10^{-12}
    2525 1.490×10−81.490\times10^{-8} −9.668×10−12-9.668\times10^{-12}

    k≤12k\le12では上界の大部分は第1項K4−kK4^{-k}であり、k≥17k\ge17では第2項(2k+2−1)u(2^{k+2}-1)uである。x=1x=1の計算値の相対誤差は18≤k≤2518\le k\le25で−9.668×10−12-9.668\times10^{-12}のまま変わらず、上界だけが増加する。

  2. X=10X=10とするとS=k+4S=k+4であり、σ=16K4−k+(2k+6−1)u\sigma=16K4^{-k}+(2^{k+6}-1)uである。k=0,1k=0,1ではσ=14.50, 1.099\sigma=14.50,\ 1.099であり、定理 4.3 条件 (a)は成り立たない。k=2,8,16,20,25k=2,8,16,20,25の上界σ/(1−σ)\sigma/(1-\sigma)は順に1.804×10−11.804\times10^{-1}、2.054×10−52.054\times10^{-5}、7.761×10−107.761\times10^{-10}、7.452×10−97.452\times10^{-9}、2.384×10−72.384\times10^{-7}であり、Ek(10)E_k(10)の相対誤差(観察)は順に2.063×10−22.063\times10^{-2}、4.967×10−64.967\times10^{-6}、7.362×10−117.362\times10^{-11}、−1.585×10−10-1.585\times10^{-10}、4.454×10−84.454\times10^{-8}である。x=10x=10ではk=23,24,25k=23,24,25で計算値の相対誤差の絶対値が10−810^{-8}を超える。

  3. X=1X=1で相対誤差10−1010^{-10}を要求する。(1+ρ)eρ>1(1+\rho)e^\rho>1かつ6(2−ρ)<126(2-\rho)<12であるからK>1/12K>1/12である。k≤14k\le14ならばσ≥4−14/12>3.1×10−10\sigma\ge4^{-14}/12>3.1\times10^{-10}であり、k≥18k\ge18ならばσ≥(220−1)u>1.16×10−10\sigma\ge(2^{20}-1)u>1.16\times10^{-10}である。k=15,16,17k=15,16,17の上界は10−1010^{-10}以下であるから、定理 4.3が要求を保証する縮小の指定はk∈{15,16,17}k\in\{15,16,17\}だけである。相対誤差10−1110^{-11}を要求すると、k≤15k\le15ならばσ≥4−15/12>7.7×10−11\sigma\ge4^{-15}/12>7.7\times10^{-11}、k=16k=16ならばσ>4.8×10−11\sigma>4.8\times10^{-11}、k≥17k\ge17ならばσ≥(219−1)u>5.8×10−11\sigma\ge(2^{19}-1)u>5.8\times10^{-11}であるから、どのkkでもσ>10−11\sigma>10^{-11}であり、定理 4.3は相対誤差10−1110^{-11}を保証しない。他方、(1)で観察した18≤k≤2518\le k\le25におけるEk(1)E_k(1)の相対誤差の絶対値9.668×10−129.668\times10^{-12}は10−1110^{-11}未満である。σ\sigmaの第1項2SKρ32^SK\rho^3の因子ρ3\rho^3は命題 2.1 (3)の局所誤差K∣x∣3K\lvert x\rvert^3の次数33から来るので、上界による保証を得るには近似の次数を上げることが一つの方針である。

注意 4.5.定理 4.3はx∈Fx\in Fを厳密な入力としてEk(x)E_k(x)をexe^xと比べる。実数ξ\xiが正規範囲にありx=fl⁡(ξ)=ξ(1+δ)x=\operatorname{fl}(\xi)=\xi(1+\delta)、∣δ∣≤u\lvert\delta\rvert\le uであるとき、ex/eξ=eξδe^x/e^\xi=e^{\xi\delta}であるから∣ex/eξ−1∣≤e∣ξ∣u−1\lvert e^x/e^\xi-1\rvert\le e^{\lvert\xi\rvert u}-1である。§E20.2 系 3.3により、ξ≠0\xi\ne0における指数関数の相対条件数は∣ξ∣\lvert\xi\rvertである。F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)で∣ξ∣≤10\lvert\xi\rvert\le10ならばe10u−1<1.2×10−15e^{10u}-1<1.2\times10^{-15}であり、この値は例 4.4 (2)のk=16k=16の上界7.761×10−107.761\times10^{-10}より小さい。

5 演習

問題 5.1.補題 3.2の証明を完成させよ。

解答.

1≤j≤k1\le j\le kに対してSj:=∑i=1jtiS_j:=\sum_{i=1}^jt_i、Pj:=∏i=1j(1−ti)P_j:=\prod_{i=1}^j(1-t_i)と置く。P1=1−t1=1−S1P_1=1-t_1=1-S_1である。1≤j<k1\le j<kとし、Pj≥1−SjP_j\ge1-S_jと仮定する。不等式Pj≥1−SjP_j\ge1-S_jの両辺に1−tj+1≥01-t_{j+1}\ge0を掛けると、Sj≥0S_j\ge0とtj+1≥0t_{j+1}\ge0により

Pj+1=Pj(1−tj+1)≥(1−Sj)(1−tj+1)=1−Sj+1+Sjtj+1≥1−Sj+1P_{j+1}=P_j(1-t_{j+1})\ge(1-S_j)(1-t_{j+1})=1-S_{j+1}+S_jt_{j+1}\ge1-S_{j+1}

である。jjに関する帰納法によりPk≥1−SkP_k\ge1-S_kが成り立つ。▨

問題 5.2. 指数関数の[2/2][2/2]型の Padé 近似p/qp/qを補題 1.3により求め、C2,2C_{2,2}が正則であることを確かめよ。さらに、任意の実数xxに対してq(x)>0q(x)>0であり、∣x∣≤1\lvert x\rvert\le1ならばq(x)≥7/12q(x)\ge7/12であることを示せ。

解答.

ci=1/i!c_i=1/i!であるから

C2,2=(c2c1c3c2)=(1/211/61/2),det⁡C2,2=14−16=112≠0C_{2,2}=\begin{pmatrix}c_2&c_1\\c_3&c_2\end{pmatrix}=\begin{pmatrix}1/2&1\\1/6&1/2\end{pmatrix},\qquad\det C_{2,2}=\frac14-\frac16=\frac1{12}\ne0

であり、C2,2C_{2,2}は正則である。補題 1.3 (1)の二つ目の等式系は

12b1+b2=−16,16b1+12b2=−124\frac12b_1+b_2=-\frac16,\qquad\frac16b_1+\frac12b_2=-\frac1{24}

である。一つ目の式からb2=−1/6−b1/2b_2=-1/6-b_1/2であり、二つ目の式に代入すると−b1/12−1/12=−1/24-b_1/12-1/12=-1/24であるからb1=−1/2b_1=-1/2、b2=1/12b_2=1/12である。一つ目の等式系からa0=c0=1a_0=c_0=1、a1=c1+b1c0=1/2a_1=c_1+b_1c_0=1/2、a2=c2+b1c1+b2c0=1/12a_2=c_2+b_1c_1+b_2c_0=1/12である。補題 1.3 (2)により、三条件を満たす組は

p(x)=1+x2+x212,q(x)=1−x2+x212p(x)=1+\frac x2+\frac{x^2}{12},\qquad q(x)=1-\frac x2+\frac{x^2}{12}

だけであり、p/qp/qが指数関数の[2/2][2/2]型の Padé 近似である。q(x)=(x−3)2/12+1/4q(x)=(x-3)^2/12+1/4であるから、任意の実数xxに対してq(x)≥1/4>0q(x)\ge1/4>0である。∣x∣≤1\lvert x\rvert\le1ならば(x−3)2≥4(x-3)^2\ge4であるからq(x)≥4/12+1/4=7/12q(x)\ge4/12+1/4=7/12であり、x=1x=1で等号が成り立つ。▨

問題 5.3.p(x):=1+x/2+x2/12p(x):=1+x/2+x^2/12、q(x):=1−x/2+x2/12q(x):=1-x/2+x^2/12と置く。任意の実数xxに対してq(x)>0q(x)>0であることを示し、r2,2(x):=p(x)/q(x)r_{2,2}(x):=p(x)/q(x)と置く。g(x):=exq(x)−p(x)g(x):=e^xq(x)-p(x)の00における Taylor 展開の最初の零でない項を求め、∣x∣≤1\lvert x\rvert\le1ならば∣g(x)∣≤(7e/1440)∣x∣5\lvert g(x)\rvert\le(7e/1440)\lvert x\rvert^5であることを示せ。さらに、実数xxを固定し、2s≥∣x∣2^s\ge\lvert x\rvertを満たすs∈N≥0s\in\Nに対してts:=2−sxt_s:=2^{-s}xと置くとき、ssによらない定数L≥0L\ge0が存在して

∣r2,2(ts)2s−ex∣≤L 2−4s\bigl\lvert r_{2,2}(t_s)^{2^s}-e^x\bigr\rvert\le L\,2^{-4s}

が成り立つことを示せ。

解答.

q(x)=((x−3)2+3)/12≥1/4q(x)=\bigl((x-3)^2+3\bigr)/12\ge1/4であるから、任意の実数xxに対してq(x)>0q(x)>0であり、r2,2r_{2,2}はR\Rの全体で定まる。∣x∣≤1\lvert x\rvert\le1ならば(x−3)2≥4(x-3)^2\ge4であるからq(x)≥7/12q(x)\ge7/12である。

q′(y)=−1/2+y/6q'(y)=-1/2+y/6、q′′(y)=1/6q''(y)=1/6であり、i≥3i\ge3ならばq(i)=0q^{(i)}=0、p(i)=0p^{(i)}=0である。Leibniz の公式により、i≥0i\ge0に対して

g(i)(y)=ey(q(y)+iq′(y)+i(i−1)2q′′(y))−p(i)(y)g^{(i)}(y)=e^y\Bigl(q(y)+iq'(y)+\frac{i(i-1)}2q''(y)\Bigr)-p^{(i)}(y)

である。q(0)=1q(0)=1、q′(0)=−1/2q'(0)=-1/2、q′′(0)=1/6q''(0)=1/6、p(0)=1p(0)=1、p′(0)=1/2p'(0)=1/2、p′′(0)=1/6p''(0)=1/6であるから

g(i)(0)=1−i2+i(i−1)12−p(i)(0)g^{(i)}(0)=1-\frac i2+\frac{i(i-1)}{12}-p^{(i)}(0)

である。i=0,1,2,3,4i=0,1,2,3,4に対してこの値は順に1−11-1、1−1/2−1/21-1/2-1/2、1−1+1/6−1/61-1+1/6-1/6、1−3/2+1/21-3/2+1/2、1−2+11-2+1であり、いずれも00である。i=5i=5に対しては1−5/2+5/3=1/61-5/2+5/3=1/6である。したがって最初の零でない項はg(5)(0)x5/5!=x5/720g^{(5)}(0)x^5/5!=x^5/720である。i=5i=5の式は

g(5)(y)=ey(1−y2+y212−52+5y6+106)=ey(y2+4y+2)12g^{(5)}(y)=e^y\Bigl(1-\frac y2+\frac{y^2}{12}-\frac52+\frac{5y}6+\frac{10}6\Bigr)=\frac{e^y(y^2+4y+2)}{12}

である。§D1.16 定理 2.1をa=0a=0とn=5n=5に適用すると、各xxに対して∣ξ∣≤∣x∣\lvert\xi\rvert\le\lvert x\rvertを満たす実数ξ\xiが存在して

g(x)=g(5)(ξ)120x5=eξ(ξ2+4ξ+2)1440x5g(x)=\frac{g^{(5)}(\xi)}{120}x^5=\frac{e^\xi(\xi^2+4\xi+2)}{1440}x^5

が成り立つ。∣x∣≤1\lvert x\rvert\le1ならば∣ξ∣≤1\lvert\xi\rvert\le1であり、∣ξ2+4ξ+2∣≤1+4+2=7\lvert\xi^2+4\xi+2\rvert\le1+4+2=7、eξ≤ee^\xi\le eであるから、∣g(x)∣≤(7e/1440)∣x∣5\lvert g(x)\rvert\le(7e/1440)\lvert x\rvert^5である。

∣t∣≤1\lvert t\rvert\le1を満たす実数ttに対してε(t):=r2,2(t)e−t−1\varepsilon(t):=r_{2,2}(t)e^{-t}-1と置くと、ε(t)=−g(t)e−t/q(t)\varepsilon(t)=-g(t)e^{-t}/q(t)である。前段の表示をx=tx=tに適用したときのξ\xiは∣ξ∣≤∣t∣≤1\lvert\xi\rvert\le\lvert t\rvert\le1を満たすから∣g(t)∣≤7e∣t∣∣t∣5/1440\lvert g(t)\rvert\le7e^{\lvert t\rvert}\lvert t\rvert^5/1440であり、e−t≤e∣t∣e^{-t}\le e^{\lvert t\rvert}とq(t)≥7/12q(t)\ge7/12により

∣ε(t)∣≤7e2∣t∣∣t∣51440⋅127=e2∣t∣120∣t∣5≤C∣t∣5,C:=e2120\lvert\varepsilon(t)\rvert\le\frac{7e^{2\lvert t\rvert}\lvert t\rvert^5}{1440}\cdot\frac{12}7=\frac{e^{2\lvert t\rvert}}{120}\lvert t\rvert^5\le C\lvert t\rvert^5,\qquad C:=\frac{e^2}{120}

である。2s≥∣x∣2^s\ge\lvert x\rvertならば∣ts∣≤1\lvert t_s\rvert\le1であり、N:=2sN:=2^sと置くと、Nts=xNt_s=xから

r2,2(ts)N=eNts(1+ε(ts))N=ex(1+ε(ts))Nr_{2,2}(t_s)^N=e^{Nt_s}\bigl(1+\varepsilon(t_s)\bigr)^N=e^x\bigl(1+\varepsilon(t_s)\bigr)^N

である。二項展開と三角不等式により∣(1+ε)N−1∣≤(1+∣ε∣)N−1\lvert(1+\varepsilon)^N-1\rvert\le(1+\lvert\varepsilon\rvert)^N-1であり、1+∣ε∣≤e∣ε∣1+\lvert\varepsilon\rvert\le e^{\lvert\varepsilon\rvert}から(1+∣ε∣)N−1≤eN∣ε∣−1(1+\lvert\varepsilon\rvert)^N-1\le e^{N\lvert\varepsilon\rvert}-1である。v≥0v\ge0に対して、平均値の定理により0≤ζ≤v0\le\zeta\le vを満たすζ\zetaが存在してev−1=veζ≤veve^v-1=ve^\zeta\le ve^vである。v:=N∣ε(ts)∣v:=N\lvert\varepsilon(t_s)\rvertと置くと

v≤2sC∣x∣52−5s=C∣x∣52−4s=C∣x∣(∣x∣2s)4≤C∣x∣v\le2^sC\lvert x\rvert^52^{-5s}=C\lvert x\rvert^52^{-4s}=C\lvert x\rvert\Bigl(\frac{\lvert x\rvert}{2^s}\Bigr)^4\le C\lvert x\rvert

であるから

∣r2,2(ts)2s−ex∣=ex∣(1+ε(ts))N−1∣≤exvev≤exC∣x∣5eC∣x∣ 2−4s\bigl\lvert r_{2,2}(t_s)^{2^s}-e^x\bigr\rvert=e^x\bigl\lvert(1+\varepsilon(t_s))^N-1\bigr\rvert\le e^xve^v\le e^xC\lvert x\rvert^5e^{C\lvert x\rvert}\,2^{-4s}

である。L:=exC∣x∣5eC∣x∣L:=e^xC\lvert x\rvert^5e^{C\lvert x\rvert}はssによらない。▨

前提記事