§E20.3安定な和・内積・多項式評価

最終更新

計算結果の誤差は、問題そのものが入力の摂動にどれだけ敏感であるかと、計算値が入力をどれだけ動かした問題の厳密な答えになっているかとの二つに分けて考えることができる。前者を表すのが条件数であり、後者を表すのが後退誤差である。

たとえば二次方程式z2+108z+1=0z^2+10^8z+1=0の絶対値の小さい根は、係数を成分ごとに相対的に摂動したときの条件数が約22であり、問題としては敏感でない。それでも解の公式(−108+1016−4)/2(-10^8+\sqrt{10^{16}-4})/2をそのまま基数22、53 桁の浮動小数点数系で計算すると、ほぼ等しい二数の差をとるために相対誤差は約0.250.25になる。一方、絶対値の大きい根を先に求めてその逆数をとると、相対誤差は単位丸め誤差と同じ程度に収まる。条件数は二つの計算に共通であるから、この違いは二つの計算値の後退誤差の違いとして説明される。

和は入力の成分について、内積は一方のベクトルを固定すれば他方の成分について、多項式の値は変数の値を固定すれば係数について線形な量であり、成分ごとの相対摂動に関する条件数と後退誤差を明示的に書き下すことができる。そのため、加算の順序の選び方や、加算で生じた丸め誤差を取り出して補正する方法を、後退誤差の上界によって比べることができる。これらの計算は、連立一次方程式や最小二乗法などの算法の内部で繰り返し用いられる。本記事では、和・内積・多項式の値の計算について、その誤差の評価と代表的な例を解説する。

1 加算木と和の誤差

定義 1.1.n∈N≥1n\in\NNとする。{1,…,n}\{1,\dots,n\}の空でない部分集合II上の 加算木 (summation tree) を、IIの元の個数に関して帰納的に次のように定める。

  1. I={i}I=\{i\}のとき、記号(i)(i)だけをII上の加算木とする。
  2. IIが二つ以上の元をもつとき、IIを交わらない空でない二つの部分集合I1,I2I_1,I_2の和集合に分け、T1T_1をI1I_1上の加算木、T2T_2をI2I_2上の加算木として、順序対T=(T1,T2)T=(T_1,T_2)をII上の加算木とする。

II上の加算木TTとi∈Ii\in Iに対して、T=(i)T=(i)のときdT(i):=0d_T(i):=0と置き、T=(T1,T2)T=(T_1,T_2)かつi∈Iki\in I_kのときdT(i):=1+dTk(i)d_T(i):=1+d_{T_k}(i)と置く。dT(i)d_T(i)をTTにおけるiiの 深さ (depth) という。

FFを浮動小数点数系、fl⁡\operatorname{fl}をFFの最近接丸めとし、x=(x1,…,xn)∈Fnx=(x_1,\dots,x_n)\in F^nとする。T=(i)T=(i)のときvT(x):=xiv_T(x):=x_iと置き、T=(T1,T2)T=(T_1,T_2)のときvT(x):=fl⁡(vT1(x)+vT2(x))v_T(x):=\operatorname{fl}(v_{T_1}(x)+v_{T_2}(x))と置く。vT(x)v_T(x)は、この帰納的な定義に現れる各実数vT1(x)+vT2(x)v_{T_1}(x)+v_{T_2}(x)(TTの加算の厳密な結果)の絶対値がNmax⁡N_{\max}以下であるときに定まる。vT(x)v_T(x)をTTによるxxの和の計算値という。

R1:=(1)R_1:=(1)、Rk:=(Rk−1,(k))R_k:=(R_{k-1},(k))(2≤k≤n2\le k\le n)と置き、vRn(x)v_{R_n}(x)をxxの 逐次和 (recursive summation) という。整数j≥1j\ge1とl≥0l\ge0に対して、Pj,0:=(j)P_{j,0}:=(j)、Pj,l:=(P2j−1,l−1,P2j,l−1)P_{j,l}:=(P_{2j-1,l-1},P_{2j,l-1})(l≥1l\ge1)と置く。llに関する帰納法により、n≥j2ln\ge j2^lのときPj,lP_{j,l}は{(j−1)2l+1,…,j2l}\{(j-1)2^l+1,\dots,j2^l\}上の加算木である。m:=⌈log⁡2n⌉m:=\lceil\log_2n\rceilと置き、xxに零を補ったx′:=(x1,…,xn,0,…,0)∈F2mx':=(x_1,\dots,x_n,0,\dots,0)\in F^{2^m}に対するvP1,m(x′)v_{P_{1,m}}(x')をxxの 対分割和 (pairwise summation) という。

例 1.2.F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)、fl⁡\operatorname{fl}を最近接偶数丸めとし、x=(1016,−1016,1)x=(10^{16},-10^{16},1)とする。253≤1016<2542^{53}\le10^{16}<2^{54}であり、§E20.1 補題 1.2 (1)によりF∩[253,254)F\cap[2^{53},2^{54})はs53=2s_{53}=2の倍数である偶数の全体であるから、±1016∈F\pm10^{16}\in Fである。逐次和はvR3(x)=fl⁡(fl⁡(1016−1016)+1)=fl⁡(0+1)=1v_{R_3}(x)=\operatorname{fl}(\operatorname{fl}(10^{16}-10^{16})+1)=\operatorname{fl}(0+1)=1であり、厳密な和に等しい。加算木T:=((1),((2),(3)))T:=((1),((2),(3)))では、−1016+1=−9999999999999999-10^{16}+1=-9999999999999999の最近接点は−1016-10^{16}と−9999999999999998-9999999999999998の二つである。1016=(5⋅1015)s5310^{16}=(5\cdot10^{15})s_{53}、9999999999999998=4999999999999999 s539999999999999998=4999999999999999\,s_{53}であり、整数5⋅10155\cdot10^{15}は偶数、49999999999999994999999999999999は奇数であるから、fl⁡(−1016+1)=−1016\operatorname{fl}(-10^{16}+1)=-10^{16}であり、vT(x)=fl⁡(1016−1016)=0v_T(x)=\operatorname{fl}(10^{16}-10^{16})=0である。−1016+1-10^{16}+1に対して−9999999999999998-9999999999999998を選ぶ最近接丸めでは、vT(x)=fl⁡(1016−9999999999999998)=2v_T(x)=\operatorname{fl}(10^{16}-9999999999999998)=2である。

定理 1.3.FFを浮動小数点数系、uuをその単位丸め誤差、fl⁡\operatorname{fl}をFFの最近接丸めとし、n∈N≥1n\in\NN、x=(x1,…,xn)∈Fnx=(x_1,\dots,x_n)\in F^n、TTを{1,…,n}\{1,\dots,n\}上の加算木とする。TTの各加算の厳密な結果の絶対値がNmax⁡N_{\max}以下であるとする。

  1. 各iiに対して、∣δi,j∣≤u|\delta_{i,j}|\le uを満たす実数δi,1,…,δi,dT(i)\delta_{i,1},\dots,\delta_{i,d_T(i)}が存在して vT(x)=∑i=1nxi∏j=1dT(i)(1+δi,j)v_T(x)=\sum_{i=1}^nx_i\prod_{j=1}^{d_T(i)}(1+\delta_{i,j}) が成り立つ。
  2. 整数d≥0d\ge0がdu<1du<1を満たし、すべてのiiについてdT(i)≤dd_T(i)\le dであるとする。0≤k≤d0\le k\le dに対してγk:=ku/(1−ku)\gamma_k:=ku/(1-ku)と置く。このとき∣θi∣≤γdT(i)≤γd|\theta_i|\le\gamma_{d_T(i)}\le\gamma_dを満たす実数θ1,…,θn\theta_1,\dots,\theta_nが存在してvT(x)=∑i=1nxi(1+θi)v_T(x)=\sum_{i=1}^nx_i(1+\theta_i)が成り立ち、 ∣vT(x)−∑i=1nxi∣≤γd∑i=1n∣xi∣\Bigl|v_T(x)-\sum_{i=1}^nx_i\Bigr|\le\gamma_d\sum_{i=1}^n|x_i| である。

証明.IIを{1,…,n}\{1,\dots,n\}の空でない部分集合、TTをII上の加算木とする。T=(i)T=(i)ならばvT(x)=xiv_T(x)=x_i、dT(i)=0d_T(i)=0であり、vT(x)=∑i∈Ixi∏j=1dT(i)(1+δi,j)v_T(x)=\sum_{i\in I}x_i\prod_{j=1}^{d_T(i)}(1+\delta_{i,j})は空積の表示として成り立つ。T=(T1,T2)T=(T_1,T_2)とし、k=1,2k=1,2について、∣δi,j∣≤u|\delta_{i,j}|\le uを満たす実数によりvTk(x)=∑i∈Ikxi∏j=1dTk(i)(1+δi,j)v_{T_k}(x)=\sum_{i\in I_k}x_i\prod_{j=1}^{d_{T_k}(i)}(1+\delta_{i,j})と書けるとする。§E20.1 系 3.3により∣δ∣≤u|\delta|\le uを満たすδ\deltaが存在して

vT(x)=(vT1(x)+vT2(x))(1+δ)=∑k=12∑i∈Ikxi(1+δ)∏j=1dTk(i)(1+δi,j)v_T(x)=\bigl(v_{T_1}(x)+v_{T_2}(x)\bigr)(1+\delta)=\sum_{k=1}^2\sum_{i\in I_k}x_i(1+\delta)\prod_{j=1}^{d_{T_k}(i)}(1+\delta_{i,j})

であり、i∈Iki\in I_kに対してdT(i)=1+dTk(i)d_T(i)=1+d_{T_k}(i)であるから、vT(x)v_T(x)も各iiの因子の個数がdT(i)d_T(i)である同じ形の表示をもつ。I1I_1、I2I_2の元の個数はIIより少ないので、IIの元の個数に関する帰納法により任意の加算木がこの表示をもち、I={1,…,n}I=\{1,\dots,n\}として(1)を得る。

(2)を示す。dT(i)≥1d_T(i)\ge1であるiiについては、dT(i)u≤du<1d_T(i)u\le du<1であるから、§E20.1 補題 4.1をk=dT(i)k=d_T(i)、εj=1\varepsilon_j=1として適用して、∏j=1dT(i)(1+δi,j)=1+θi\prod_{j=1}^{d_T(i)}(1+\delta_{i,j})=1+\theta_i、∣θi∣≤γdT(i)|\theta_i|\le\gamma_{d_T(i)}を満たすθi\theta_iを得る。dT(i)=0d_T(i)=0であるiiについてはθi:=0=γ0\theta_i:=0=\gamma_0と置く。k↦ku/(1−ku)k\mapsto ku/(1-ku)は0≤k≤d0\le k\le dで単調非減少であるからγdT(i)≤γd\gamma_{d_T(i)}\le\gamma_dであり、∣vT(x)−∑ixi∣=∣∑ixiθi∣≤γd∑i∣xi∣|v_T(x)-\sum_ix_i|=|\sum_ix_i\theta_i|\le\gamma_d\sum_i|x_i|が成り立つ。▨

系 1.4.FFを浮動小数点数系、uuをその単位丸め誤差、fl⁡\operatorname{fl}をFFの最近接丸めとし、n∈N≥1n\in\NN、x∈Fnx\in F^nとする。整数k≥0k\ge0がku<1ku<1を満たすときγk:=ku/(1−ku)\gamma_k:=ku/(1-ku)と置く。

  1. dRn(1)=n−1d_{R_n}(1)=n-1であり、2≤i≤n2\le i\le nに対してdRn(i)=n−i+1d_{R_n}(i)=n-i+1である。(n−1)u<1(n-1)u<1であり、逐次和の各加算の厳密な結果の絶対値がNmax⁡N_{\max}以下であるならば、∣vRn(x)−∑ixi∣≤γn−1∑i∣xi∣|v_{R_n}(x)-\sum_ix_i|\le\gamma_{n-1}\sum_i|x_i|が成り立つ。
  2. m:=⌈log⁡2n⌉m:=\lceil\log_2n\rceilと置く。P1,mP_{1,m}におけるすべての元の深さはmmである。mu<1mu<1であり、対分割和の各加算の厳密な結果の絶対値がNmax⁡N_{\max}以下であるならば、対分割和s^\hat sは∣s^−∑ixi∣≤γm∑i∣xi∣|\hat s-\sum_ix_i|\le\gamma_m\sum_i|x_i|を満たす。

証明.R1=(1)R_1=(1)ではdR1(1)=0d_{R_1}(1)=0である。n≥2n\ge2ならばRn=(Rn−1,(n))R_n=(R_{n-1},(n))であるから、dRn(n)=1d_{R_n}(n)=1であり、i≤n−1i\le n-1に対してdRn(i)=1+dRn−1(i)d_{R_n}(i)=1+d_{R_{n-1}}(i)である。nnに関する帰納法により深さの式が成り立ち、その最大値はn−1n-1である。Pj,0=(j)P_{j,0}=(j)の元の深さは00であり、Pj,l=(P2j−1,l−1,P2j,l−1)P_{j,l}=(P_{2j-1,l-1},P_{2j,l-1})の元の深さはP2j−1,l−1P_{2j-1,l-1}またはP2j,l−1P_{2j,l-1}における深さに11を加えたものであるから、llに関する帰納法によりPj,lP_{j,l}の元の深さはすべてllである。零を補ったx′x'は∑ixi′=∑ixi\sum_ix'_i=\sum_ix_iと∑i∣xi′∣=∑i∣xi∣\sum_i|x'_i|=\sum_i|x_i|を満たす。定理 1.3 (2)を、逐次和にはd=n−1d=n-1で、対分割和にはx′x'とd=md=mで適用する。▨

例 1.5.F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)、fl⁡\operatorname{fl}を最近接偶数丸め、u=2−53u=2^{-53}とし、x:=(1,2−53,…,2−53)∈F8x:=(1,2^{-53},\dots,2^{-53})\in F^8(2−532^{-53}が七つ)とする。厳密な和は1+7u1+7uであり、全成分が正であるから∑i∣xi∣=1+7u\sum_i|x_i|=1+7uである。1+2−531+2^{-53}の最近接点は1=252s01=2^{52}s_0と1+2−52=(252+1)s01+2^{-52}=(2^{52}+1)s_0であり、整数2522^{52}が偶数であるからfl⁡(1+2−53)=1\operatorname{fl}(1+2^{-53})=1である。したがって逐次和の各段の値は11であり、vR8(x)=1v_{R_8}(x)=1、絶対誤差は7u≈7.77×10−167u\approx7.77\times10^{-16}である。系 1.4 (1)の上界はγ7(1+7u)=7u(1+7u)/(1−7u)\gamma_7(1+7u)=7u(1+7u)/(1-7u)であり、誤差との比は(1−7u)/(1+7u)>1−14u(1-7u)/(1+7u)>1-14uである。対分割和ではm=3m=3であり、第一段の値はfl⁡(1+2−53)=1\operatorname{fl}(1+2^{-53})=1と三つの2−53+2−53=2−522^{-53}+2^{-53}=2^{-52}、第二段の値は1+2−521+2^{-52}と2−512^{-51}、第三段の値は1+2−52+2−51=1+3⋅2−52∈F1+2^{-52}+2^{-51}=1+3\cdot2^{-52}\in Fである。絶対誤差はu≈1.11×10−16u\approx1.11\times10^{-16}であり、系 1.4 (2)の上界γ3(1+7u)≈3.33×10−16\gamma_3(1+7u)\approx3.33\times10^{-16}以下である。

2 和の条件数

補題 2.1.n∈N≥1n\in\NN、c∈Rnc\in\R^nとし、L ⁣:Rn→RL\colon\R^n\to\RをL(x):=∑i=1ncixiL(x):=\sum_{i=1}^nc_ix_iで定める。x∈Rnx\in\R^nに対してM:=∑i=1n∣cixi∣M:=\sum_{i=1}^n|c_ix_i|と置き、sgn⁡(0):=0\operatorname{sgn}(0):=0とし、δ∈Rn\delta\in\R^nに対するx∘δx\circ\deltaを§E20.2 定義 5.4の成分ごとの積とする。

  1. L(x)≠0L(x)\ne0ならばκcomp(L,x)=M/∣L(x)∣\kappa_{\mathrm{comp}}(L,x)=M/|L(x)|である。
  2. M>0M>0ならば、任意のy^∈R\hat y\in\Rに対してηcomp(L,x,y^)=∣y^−L(x)∣/M\eta_{\mathrm{comp}}(L,x,\hat y)=|\hat y-L(x)|/Mである。さらにδi:=sgn⁡(cixi)(y^−L(x))/M\delta_i:=\operatorname{sgn}(c_ix_i)(\hat y-L(x))/Mで定まるδ∈Rn\delta\in\R^nは∥δ∥∞=∣y^−L(x)∣/M\lVert\delta\rVert_\infty=|\hat y-L(x)|/MとL(x+x∘δ)=y^L(x+x\circ\delta)=\hat yを満たす。

証明.LLの定義域はRn\R^nであるから§E20.2 定義 5.4のUxU_xはRn\R^nであり、任意のδ∈Rn\delta\in\R^nに対して

Gx(δ)−Gx(0)=∑i=1ncixiδi,∣Gx(δ)−Gx(0)∣≤M∥δ∥∞G_x(\delta)-G_x(0)=\sum_{i=1}^nc_ix_i\delta_i,\qquad |G_x(\delta)-G_x(0)|\le M\lVert\delta\rVert_\infty

が成り立つ。

(1)を示す。M≥∣L(x)∣>0M\ge|L(x)|>0である。ε>0\varepsilon>0に対して、§E20.2 定義 3.1の上限s(ε)s(\varepsilon)は上の不等式によりMM以下である。δ:=ε(sgn⁡(cixi))i\delta:=\varepsilon(\operatorname{sgn}(c_ix_i))_iは∥δ∥∞=ε\lVert\delta\rVert_\infty=\varepsilonとGx(δ)−Gx(0)=εMG_x(\delta)-G_x(0)=\varepsilon Mを満たすから、s(ε)=Ms(\varepsilon)=Mである。したがってκabs(Gx,0)=M\kappa_{\mathrm{abs}}(G_x,0)=Mである。

(2)を示す。L(x+x∘δ)=y^L(x+x\circ\delta)=\hat yならば、∣y^−L(x)∣=∣Gx(δ)−Gx(0)∣≤M∥δ∥∞|\hat y-L(x)|=|G_x(\delta)-G_x(0)|\le M\lVert\delta\rVert_\inftyであるから、∥δ∥∞≥∣y^−L(x)∣/M\lVert\delta\rVert_\infty\ge|\hat y-L(x)|/Mである。主張のδ\deltaは∑icixiδi=∑i∣cixi∣(y^−L(x))/M=y^−L(x)\sum_ic_ix_i\delta_i=\sum_i|c_ix_i|(\hat y-L(x))/M=\hat y-L(x)を満たし、M>0M>0であるからcixi≠0c_ix_i\ne0であるiiが存在して∥δ∥∞=∣y^−L(x)∣/M\lVert\delta\rVert_\infty=|\hat y-L(x)|/Mである。▨

系 2.2.n∈N≥1n\in\NNとし、σ ⁣:Rn→R\sigma\colon\R^n\to\Rをσ(x):=∑i=1nxi\sigma(x):=\sum_{i=1}^nx_iで定める。

  1. σ(x)≠0\sigma(x)\ne0ならばκcomp(σ,x)=∑i∣xi∣/∣∑ixi∣\kappa_{\mathrm{comp}}(\sigma,x)=\sum_i|x_i|/|\sum_ix_i|である。
  2. FF、uu、fl⁡\operatorname{fl}、x∈Fnx\in F^n、TT、dd、γd\gamma_dが定理 1.3と定理 1.3 (2)の仮定を満たすならば、ηcomp(σ,x,vT(x))≤γd\eta_{\mathrm{comp}}(\sigma,x,v_T(x))\le\gamma_dである。さらにσ(x)≠0\sigma(x)\ne0ならば ∣vT(x)−σ(x)∣∣σ(x)∣=κcomp(σ,x) ηcomp(σ,x,vT(x))≤γd κcomp(σ,x)\frac{|v_T(x)-\sigma(x)|}{|\sigma(x)|}=\kappa_{\mathrm{comp}}(\sigma,x)\,\eta_{\mathrm{comp}}(\sigma,x,v_T(x))\le\gamma_d\,\kappa_{\mathrm{comp}}(\sigma,x) が成り立つ。

証明.(1)は補題 2.1 (1)をc=(1,…,1)c=(1,\dots,1)として適用したものである。定理 1.3 (2)のθ=(θ1,…,θn)\theta=(\theta_1,\dots,\theta_n)はσ(x+x∘θ)=vT(x)\sigma(x+x\circ\theta)=v_T(x)と∥θ∥∞≤γd\lVert\theta\rVert_\infty\le\gamma_dを満たすので、ηcomp(σ,x,vT(x))≤γd\eta_{\mathrm{comp}}(\sigma,x,v_T(x))\le\gamma_dである。σ(x)≠0\sigma(x)\ne0ならば∑i∣xi∣>0\sum_i|x_i|>0であり、補題 2.1 (2)によりηcomp(σ,x,vT(x))=∣vT(x)−σ(x)∣/∑i∣xi∣\eta_{\mathrm{comp}}(\sigma,x,v_T(x))=|v_T(x)-\sigma(x)|/\sum_i|x_i|であるから、(1)と合わせて等式が成り立つ。▨

例 2.3.F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)、fl⁡\operatorname{fl}を最近接偶数丸め、u=2−53u=2^{-53}とし、x:=(1016,1,1,−1016)x:=(10^{16},1,1,-10^{16})とする。σ(x)=2\sigma(x)=2、∑i∣xi∣=2⋅1016+2\sum_i|x_i|=2\cdot10^{16}+2であり、系 2.2 (1)によりκcomp(σ,x)=1016+1\kappa_{\mathrm{comp}}(\sigma,x)=10^{16}+1である。1016+110^{16}+1の最近接点は101610^{16}と1016+210^{16}+2であり、1016=(5⋅1015)s5310^{16}=(5\cdot10^{15})s_{53}、1016+2=(5⋅1015+1)s5310^{16}+2=(5\cdot10^{15}+1)s_{53}であり、整数5⋅10155\cdot10^{15}が偶数であるからfl⁡(1016+1)=1016\operatorname{fl}(10^{16}+1)=10^{16}である。したがって逐次和はfl⁡(fl⁡(1016+1)+1)=1016\operatorname{fl}(\operatorname{fl}(10^{16}+1)+1)=10^{16}を経てvR4(x)=0v_{R_4}(x)=0である。対分割和は、例 1.2と同じくfl⁡(1−1016)=−1016\operatorname{fl}(1-10^{16})=-10^{16}であるから、fl⁡(1016−1016)=0\operatorname{fl}(10^{16}-10^{16})=0である。どちらの計算値も相対誤差は11である。系 2.2 (2)の上界は、逐次和でγ3κcomp(σ,x)≈3.33\gamma_3\kappa_{\mathrm{comp}}(\sigma,x)\approx3.33、対分割和でγ2κcomp(σ,x)≈2.22\gamma_2\kappa_{\mathrm{comp}}(\sigma,x)\approx2.22であり、絶対誤差の上界はγ3(2⋅1016+2)≈6.66\gamma_3(2\cdot10^{16}+2)\approx6.66とγ2(2⋅1016+2)≈4.44\gamma_2(2\cdot10^{16}+2)\approx4.44である。

3 丸め誤差の回収

補題 3.1.F=F(β,p,emin⁡,emax⁡)F=F(\beta,p,e_{\min},e_{\max})を浮動小数点数系(非正規数を含む)、fl⁡\operatorname{fl}をFFの最近接丸めとする。x,y∈Fx,y\in Fがy/2≤x≤2yy/2\le x\le2yを満たすならば、x−y∈Fx-y\in Fであり、fl⁡(x−y)=x−y\operatorname{fl}(x-y)=x-yである。

証明.y/2≤2yy/2\le2yからy≥0y\ge0であり、x≥y/2≥0x\ge y/2\ge0である。条件y/2≤x≤2yy/2\le x\le2yはx/2≤y≤2xx/2\le y\le2xと同値であり、F=−FF=-Fであるから、x≥yx\ge yとしてよい。x≤2yx\le2yから0≤x−y≤y0\le x-y\le yである。y=0y=0ならばx=0x=0であり、x−y=0∈Fx-y=0\in Fである。y>0y>0とし、§E20.1 補題 1.2の整数e(y)e(y)をeeと置く。§E20.1 補題 1.3 (3)によりe(x)≥ee(x)\ge eであり、§E20.1 補題 1.3 (2)によりxxとyyはses_eの整数倍であって、y=msey=ms_e、1≤m≤βp−11\le m\le\beta^p-1と書ける。したがってx−y=ksex-y=ks_eを満たす整数kkがあり、0≤x−y≤y0\le x-y\le yから0≤k≤m≤βp−10\le k\le m\le\beta^p-1である。∣x−y∣≤y≤Nmax⁡|x-y|\le y\le N_{\max}であるから、§E20.1 補題 1.3 (1)によりx−y∈Fx-y\in Fであり、§E20.1 系 3.2 (1)によりfl⁡(x−y)=x−y\operatorname{fl}(x-y)=x-yである。▨

注意 3.2.p≥2p\ge2とし、F=F(β,p,emin⁡,emax⁡)F=F(\beta,p,e_{\min},e_{\max})から非正規数を除いた集合をF′F'とする。x:=(βp−1+1)semin⁡x:=(\beta^{p-1}+1)s_{e_{\min}}、y:=βp−1semin⁡=βemin⁡y:=\beta^{p-1}s_{e_{\min}}=\beta^{e_{\min}}はF′F'の元であり、y/2≤x≤2yy/2\le x\le2yを満たす。x−y=semin⁡x-y=s_{e_{\min}}は非正規数であるからF′F'に属さず、補題 3.1の結論はF′F'では成り立たない。

補題 3.3.F=F(β,p,emin⁡,emax⁡)F=F(\beta,p,e_{\min},e_{\max})を浮動小数点数系とし、a,b∈Fa,b\in Fが∣a+b∣≤Nmax⁡|a+b|\le N_{\max}を満たすとする。s∈Fs\in Fがa+ba+bの最近接点ならば、r:=a+b−sr:=a+b-sはFFに属し、∣r∣≤min⁡{∣a∣,∣b∣}|r|\le\min\{|a|,|b|\}を満たす。

証明.t:=a+bt:=a+bと置く。a,b∈Fa,b\in Fであるから∣r∣=∣t−s∣≤∣t−a∣=∣b∣|r|=|t-s|\le|t-a|=|b|かつ∣r∣≤∣t−b∣=∣a∣|r|\le|t-b|=|a|である。主張はaaとbbについて対称であるから、∣b∣≤∣a∣|b|\le|a|としてよい。b=0b=0ならばt=a∈Ft=a\in Fであり、ttとの距離が00であるFFの元はttだけであるからs=ts=t、r=0r=0である。b≠0b\ne0とし、e:=e(∣b∣)e:=e(|b|)と置く。§E20.1 補題 1.3 (3)によりe(∣a∣)≥ee(|a|)\ge eであり、§E20.1 補題 1.3 (2)によりaaとbbはses_eの整数倍であって、∣b∣=mse|b|=ms_e、1≤m≤βp−11\le m\le\beta^p-1と書ける。したがってttもses_eの整数倍である。ssがses_eの整数倍でないと仮定する。このときs≠0s\ne0であり、§E20.1 補題 1.3 (2)によりe(∣s∣)<ee(|s|)<eである。e(⋅)e(\cdot)の定義により∣s∣<βe(∣s∣)+1≤βe|s|<\beta^{e(|s|)+1}\le\beta^eである。∣t∣≥βe|t|\ge\beta^eならば、v:=sgn⁡(t)βe=sgn⁡(t)βp−1sev:=\operatorname{sgn}(t)\beta^e=\operatorname{sgn}(t)\beta^{p-1}s_eはFFの元であり、∣t−v∣=∣t∣−βe<∣t∣−∣s∣≤∣t−s∣|t-v|=|t|-\beta^e<|t|-|s|\le|t-s|を満たすので、ssが最近接点であることに反する。よって∣t∣<βe=βp−1se|t|<\beta^e=\beta^{p-1}s_eであり、t=kset=ks_e、∣k∣<βp−1|k|<\beta^{p-1}を満たす整数kkについて§E20.1 補題 1.3 (1)によりt∈Ft\in Fである。このときs=ts=tはses_eの整数倍であり、仮定に反する。したがってssはses_eの整数倍であり、r=kser=ks_eを満たす整数kkがある。∣r∣≤∣b∣|r|\le|b|から∣k∣≤m≤βp−1|k|\le m\le\beta^p-1であり、∣r∣≤∣b∣≤Nmax⁡|r|\le|b|\le N_{\max}であるから、§E20.1 補題 1.3 (1)によりr∈Fr\in Fである。▨

補題 3.4.F=F(2,p,emin⁡,emax⁡)F=F(2,p,e_{\min},e_{\max})を基数22の浮動小数点数系とし、x,y∈Fx,y\in Fが∣x∣≥∣y∣|x|\ge|y|と∣x+y∣≤Nmax⁡|x+y|\le N_{\max}を満たすとする。s∈Fs\in Fをx+yx+yの最近接点とする。

  1. s−x∈Fs-x\in Fである。
  2. x+y∉Fx+y\notin Fならば∣s∣≥∣y∣|s|\ge|y|である。

証明.t:=x+yt:=x+yと置く。F=−FF=-Fであり、−s-sは−t-tの最近接点であるから、(x,y,s)(x,y,s)を(−x,−y,−s)(-x,-y,-s)に替えても仮定と結論は変わらない。したがってx≥0x\ge0としてよい。t∈Ft\in Fならば、ttとの距離が00であるFFの元はttだけであるからs=ts=tであり、s−x=y∈Fs-x=y\in Fである。このとき(2)の仮定は成り立たない。以下t∉Ft\notin Fとする。このときy≠0y\ne0であり、x≥∣y∣>0x\ge|y|>0である。v∈Fv\in Fがv≤tv\le tを満たすならばv≤sv\le sである。実際、s<v≤ts<v\le tならば∣v−t∣<∣s−t∣|v-t|<|s-t|となり、ssが最近接点であることに反する。同様に、v∈Fv\in Fがv≥tv\ge tを満たすならばv≥sv\ge sである。§E20.1 補題 1.3 (2)によりx=msex=ms_e、e:=e(x)e:=e(x)、1≤m≤2p−11\le m\le2^p-1と書く。

y>0y>0の場合、x<t≤2xx<t\le2xであり、v=xv=xとしてs≥xs\ge xである。2x≤Nmax⁡2x\le N_{\max}ならば2x∈F2x\in Fである。実際、m≤2p−1m\le2^{p-1}ならば2x=(2m)se2x=(2m)s_e、2m≤2p2m\le2^pであり、m>2p−1m>2^{p-1}ならばx≥2p−1se=2ex\ge2^{p-1}s_e=2^eかつ2x=mse+12x=ms_{e+1}であって、2e+1≤2x≤Nmax⁡<2emax⁡+12^{e+1}\le2x\le N_{\max}<2^{e_{\max}+1}からe+1≤emax⁡e+1\le e_{\max}であるので、どちらの場合も§E20.1 補題 1.3 (1)が適用される。このときv=2xv=2xとしてs≤2xs\le2xである。2x>Nmax⁡2x>N_{\max}ならばs≤Nmax⁡<2xs\le N_{\max}<2xである。したがってx/2≤x≤s≤2xx/2\le x\le s\le2xであり、補題 3.1によりs−x∈Fs-x\in Fである。また∣s∣=s≥x≥∣y∣|s|=s\ge x\ge|y|である。

y<0y<0の場合、0≤t<x0\le t<xである。−y≥x/2-y\ge x/2ならば、(−y)/2≤x≤2(−y)(-y)/2\le x\le2(-y)であるから補題 3.1によりt=x−(−y)∈Ft=x-(-y)\in Fとなり、t∉Ft\notin Fに反する。よって−y<x/2-y<x/2であり、x/2<t<xx/2<t<xである。e=emin⁡e=e_{\min}ならば、§E20.1 補題 1.3 (2)によりyyもsemin⁡s_{e_{\min}}の整数倍であるからt=ksemin⁡t=ks_{e_{\min}}を満たす整数kkがあり、0<t<x=msemin⁡0<t<x=ms_{e_{\min}}から0<k<m≤2p−10<k<m\le2^p-1であって、§E20.1 補題 1.3 (1)によりt∈Ft\in Fとなり、t∉Ft\notin Fに反する。よってe>emin⁡e>e_{\min}であり、x/2=mse−1x/2=ms_{e-1}は§E20.1 補題 1.3 (1)によりFFに属する。v=x/2v=x/2とv=xv=xとしてx/2≤s≤xx/2\le s\le xであり、補題 3.1によりs−x∈Fs-x\in Fである。また∣s∣=s≥x/2>−y=∣y∣|s|=s\ge x/2>-y=|y|である。▨

定義 3.5.FFを浮動小数点数系、fl⁡\operatorname{fl}をFFの最近接丸めとし、a,b∈Fa,b\in Fとする。

s:=fl⁡(a+b),z:=fl⁡(s−a),a′:=fl⁡(s−z),e:=fl⁡(fl⁡(a−a′)+fl⁡(b−z))s:=\operatorname{fl}(a+b),\quad z:=\operatorname{fl}(s-a),\quad a':=\operatorname{fl}(s-z),\quad e:=\operatorname{fl}\bigl(\operatorname{fl}(a-a')+\operatorname{fl}(b-z)\bigr)

と置く。これら六つの演算の厳密な結果の絶対値がすべてNmax⁡N_{\max}以下であるとき、対(s,e)(s,e)を(a,b)(a,b)の TwoSum (TwoSum) という。

定理 3.6.p≥1p\ge1、emin⁡≤emax⁡e_{\min}\le e_{\max}を整数とし、F=F(2,p,emin⁡,emax⁡)F=F(2,p,e_{\min},e_{\max})を基数22の浮動小数点数系、uuをその単位丸め誤差、fl⁡\operatorname{fl}をFFの最近接丸めとする。a,b∈Fa,b\in Fが∣a+b∣≤Nmax⁡|a+b|\le N_{\max}と∣fl⁡(a+b)−a∣≤Nmax⁡|\operatorname{fl}(a+b)-a|\le N_{\max}を満たすとし、ss、zz、a′a'、eeを定義 3.5の値とする。

  1. s−zs-z、a−a′a-a'、b−zb-z、fl⁡(a−a′)+fl⁡(b−z)\operatorname{fl}(a-a')+\operatorname{fl}(b-z)はすべてFFに属する。したがって(a,b)(a,b)の TwoSum(s,e)(s,e)が定まり、a′=s−za'=s-z、fl⁡(a−a′)=a−a′\operatorname{fl}(a-a')=a-a'、fl⁡(b−z)=b−z\operatorname{fl}(b-z)=b-z、e=(a−a′)+(b−z)e=(a-a')+(b-z)が成り立つ。
  2. a+b=s+ea+b=s+eが成り立つ。
  3. ∣e∣≤min⁡{∣a∣,∣b∣}|e|\le\min\{|a|,|b|\}かつ∣e∣≤u∣a+b∣|e|\le u|a+b|である。

証明.t:=a+bt:=a+b、r:=t−sr:=t-sと置く。s=fl⁡(t)s=\operatorname{fl}(t)はttの最近接点であるから、補題 3.3によりr∈Fr\in F、∣r∣≤min⁡{∣a∣,∣b∣}|r|\le\min\{|a|,|b|\}である。

∣a∣≥∣b∣|a|\ge|b|の場合。補題 3.4 (1)を(x,y)=(a,b)(x,y)=(a,b)に適用してs−a∈Fs-a\in Fであるから、z=s−az=s-aである。s−z=a∈Fs-z=a\in Fであるからa′=aa'=aであり、a−a′=0∈Fa-a'=0\in F、b−z=t−s=r∈Fb-z=t-s=r\in Fである。(a−a′)+(b−z)=r(a-a')+(b-z)=rである。

∣a∣<∣b∣|a|<|b|かつt∈Ft\in Fの場合。s=ts=tであるからs−a=b∈Fs-a=b\in Fであり、z=bz=bである。s−z=a∈Fs-z=a\in Fであるからa′=aa'=aであり、a−a′=0a-a'=0、b−z=0b-z=0はともにFFに属し、(a−a′)+(b−z)=0=r(a-a')+(b-z)=0=rである。

∣a∣<∣b∣|a|<|b|かつt∉Ft\notin Fの場合。補題 3.4 (2)を(x,y)=(b,a)(x,y)=(b,a)に適用して∣s∣≥∣a∣|s|\ge|a|である。z=fl⁡(s+(−a))z=\operatorname{fl}(s+(-a))はs+(−a)s+(-a)の最近接点であり、s,−a∈Fs,-a\in F、∣s∣≥∣−a∣|s|\ge|-a|、∣s−a∣≤Nmax⁡|s-a|\le N_{\max}であるから、補題 3.4 (1)を(x,y)=(s,−a)(x,y)=(s,-a)に適用してz−s∈Fz-s\in Fである。よってs−z∈Fs-z\in Fであり、a′=s−za'=s-zである。補題 3.3をssと−a-aに適用してρ:=(s−a)−z∈F\rho:=(s-a)-z\in Fであり、a−a′=a−s+z=−ρ∈Fa-a'=a-s+z=-\rho\in Fである。実数としてs−a=b−rs-a=b-rであるから、zzはb+(−r)b+(-r)の最近接点でもある。b,−r∈Fb,-r\in F、∣−r∣≤∣a∣<∣b∣|-r|\le|a|<|b|、∣b−r∣=∣s−a∣≤Nmax⁡|b-r|=|s-a|\le N_{\max}であるから、補題 3.4 (1)を(x,y)=(b,−r)(x,y)=(b,-r)に適用してz−b∈Fz-b\in Fであり、b−z∈Fb-z\in Fである。(a−a′)+(b−z)=(a−s+z)+(b−z)=t−s=r(a-a')+(b-z)=(a-s+z)+(b-z)=t-s=rである。

三つの場合のいずれでもs−zs-z、a−a′a-a'、b−zb-zはFFに属し、(a−a′)+(b−z)=r∈F(a-a')+(b-z)=r\in Fである。§E20.1 系 3.2 (1)により各演算の結果は厳密な結果に等しく、e=fl⁡(r)=re=\operatorname{fl}(r)=rであるから、(1)と(2)が成り立つ。

(3)を示す。e=re=rであるから∣e∣≤min⁡{∣a∣,∣b∣}|e|\le\min\{|a|,|b|\}である。ttが正規範囲にあるならば§E20.1 系 3.2 (2)により∣e∣=∣t−fl⁡(t)∣≤u∣t∣|e|=|t-\operatorname{fl}(t)|\le u|t|であり、∣t∣<2emin⁡|t|<2^{e_{\min}}ならば§E20.1 系 3.2 (4)によりs=ts=t、e=0e=0である。▨

注意 3.7.F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)とし、b:=Nmax⁡=(253−1)2971b:=N_{\max}=(2^{53}-1)2^{971}、a:=−2970a:=-2^{970}とする。a+b=Nmax⁡−2970a+b=N_{\max}-2^{970}は∣a+b∣≤Nmax⁡|a+b|\le N_{\max}を満たし、§E20.1 補題 1.2 (1)によりFFの隣接する二元Nmax⁡−2971N_{\max}-2^{971}とNmax⁡N_{\max}の中点である。fl⁡(Nmax⁡−2970)=Nmax⁡\operatorname{fl}(N_{\max}-2^{970})=N_{\max}を満たす最近接丸めfl⁡\operatorname{fl}では、fl⁡(a+b)−a=Nmax⁡+2970>Nmax⁡\operatorname{fl}(a+b)-a=N_{\max}+2^{970}>N_{\max}であり、定理 3.6の第二の仮定は第一の仮定から従わない。最近接偶数丸めでは、Nmax⁡=(253−1)s1023N_{\max}=(2^{53}-1)s_{1023}の整数253−12^{53}-1が奇数、Nmax⁡−2971=(253−2)s1023N_{\max}-2^{971}=(2^{53}-2)s_{1023}の整数253−22^{53}-2が偶数であるからfl⁡(a+b)=Nmax⁡−2971\operatorname{fl}(a+b)=N_{\max}-2^{971}であり、fl⁡(a+b)−a=Nmax⁡−2970\operatorname{fl}(a+b)-a=N_{\max}-2^{970}は第二の仮定を満たす。

4 補償和

定義 4.1.F=F(2,p,emin⁡,emax⁡)F=F(2,p,e_{\min},e_{\max})を基数22の浮動小数点数系、fl⁡\operatorname{fl}をFFの最近接丸めとし、n≥2n\ge2を整数、x=(x1,…,xn)∈Fnx=(x_1,\dots,x_n)\in F^nとする。s1:=x1s_1:=x_1と置き、i=2,…,ni=2,\dots,nについて(si−1,xi)(s_{i-1},x_i)の TwoSum を(si,ei)(s_i,e_i)とする。(e2,…,en)∈Fn−1(e_2,\dots,e_n)\in F^{n-1}の逐次和をccと置き、s^:=fl⁡(sn+c)\hat s:=\operatorname{fl}(s_n+c)をxxの 補償和 (compensated summation) という。s^\hat sは、各 TwoSum とccの計算の各加算と最後の加算の厳密な結果の絶対値がNmax⁡N_{\max}以下であるときに定まる。

定理 4.2.F=F(2,p,emin⁡,emax⁡)F=F(2,p,e_{\min},e_{\max})を基数22の浮動小数点数系、uuをその単位丸め誤差、fl⁡\operatorname{fl}をFFの最近接丸めとし、n≥2n\ge2を整数、x∈Fnx\in F^nとする。各i=2,…,ni=2,\dots,nについて(a,b)=(si−1,xi)(a,b)=(s_{i-1},x_i)が定理 3.6の仮定を満たし、ccの計算の各加算とsn+cs_n+cの厳密な結果の絶対値がNmax⁡N_{\max}以下であるとする。整数k≥0k\ge0がku<1ku<1を満たすときγk:=ku/(1−ku)\gamma_k:=ku/(1-ku)と置く。

  1. ∑i=1nxi=sn+∑i=2nei\sum_{i=1}^nx_i=s_n+\sum_{i=2}^ne_iが成り立つ。
  2. (n−2)u<1(n-2)u<1ならば、補償和s^\hat sは ∣s^−∑i=1nxi∣≤u∣∑i=1nxi∣+(1+u)γn−2∑i=2n∣ei∣\Bigl|\hat s-\sum_{i=1}^nx_i\Bigr|\le u\Bigl|\sum_{i=1}^nx_i\Bigr|+(1+u)\gamma_{n-2}\sum_{i=2}^n|e_i| を満たす。
  3. (n−1)u<1(n-1)u<1ならば∑i=2n∣ei∣≤γn−1∑i=1n∣xi∣\sum_{i=2}^n|e_i|\le\gamma_{n-1}\sum_{i=1}^n|x_i|であり、 ∣s^−∑i=1nxi∣≤u∣∑i=1nxi∣+(1+u)γn−2γn−1∑i=1n∣xi∣\Bigl|\hat s-\sum_{i=1}^nx_i\Bigr|\le u\Bigl|\sum_{i=1}^nx_i\Bigr|+(1+u)\gamma_{n-2}\gamma_{n-1}\sum_{i=1}^n|x_i| が成り立つ。

証明.定理 3.6 (2)によりsi−1+xi=si+eis_{i-1}+x_i=s_i+e_i(2≤i≤n2\le i\le n)であり、これらを加えてs1=x1s_1=x_1を用いると(1)を得る。

(2)を示す。S:=∑ixiS:=\sum_ix_i、E:=∑i=2neiE:=\sum_{i=2}^ne_iと置く。ccはn−1n-1個の成分の逐次和であるから、系 1.4 (1)により∣c−E∣≤γn−2∑i∣ei∣|c-E|\le\gamma_{n-2}\sum_i|e_i|である。§E20.1 系 3.3により∣δ∣≤u|\delta|\le uを満たすδ\deltaが存在してs^=(sn+c)(1+δ)\hat s=(s_n+c)(1+\delta)であり、(1)によりsn=S−Es_n=S-Eであるから

s^−S=(S+(c−E))(1+δ)−S=δS+(1+δ)(c−E)\hat s-S=(S+(c-E))(1+\delta)-S=\delta S+(1+\delta)(c-E)

である。したがって∣s^−S∣≤u∣S∣+(1+u)γn−2∑i∣ei∣|\hat s-S|\le u|S|+(1+u)\gamma_{n-2}\sum_i|e_i|である。

(3)を示す。定理 3.6 (3)により∣ei∣≤u∣si−1+xi∣|e_i|\le u|s_{i-1}+x_i|である。si=fl⁡(si−1+xi)s_i=\operatorname{fl}(s_{i-1}+x_i)であるから、si−1s_{i-1}は(x1,…,xi−1)(x_1,\dots,x_{i-1})の逐次和であり、系 1.4 (1)により∣si−1−∑j<ixj∣≤γi−2∑j<i∣xj∣|s_{i-1}-\sum_{j<i}x_j|\le\gamma_{i-2}\sum_{j<i}|x_j|である。したがって

∣si−1+xi∣≤(1+γi−2)∑j≤i∣xj∣≤(1+γn−2)∑j=1n∣xj∣|s_{i-1}+x_i|\le(1+\gamma_{i-2})\sum_{j\le i}|x_j|\le(1+\gamma_{n-2})\sum_{j=1}^n|x_j|

であり、i=2,…,ni=2,\dots,nについて加えて

∑i=2n∣ei∣≤(n−1)u(1+γn−2)∑j∣xj∣=(n−1)u1−(n−2)u∑j∣xj∣≤(n−1)u1−(n−1)u∑j∣xj∣=γn−1∑j∣xj∣\sum_{i=2}^n|e_i|\le(n-1)u(1+\gamma_{n-2})\sum_j|x_j|=\frac{(n-1)u}{1-(n-2)u}\sum_j|x_j|\le\frac{(n-1)u}{1-(n-1)u}\sum_j|x_j|=\gamma_{n-1}\sum_j|x_j|

を得る。これを(2)に代入する。▨

例 4.3.F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)、fl⁡\operatorname{fl}を最近接偶数丸めとし、例 2.3のx=(1016,1,1,−1016)x=(10^{16},1,1,-10^{16})の補償和を計算する。si=fl⁡(si−1+xi)s_i=\operatorname{fl}(s_{i-1}+x_i)は逐次和の部分和であり、s2=s3=1016s_2=s_3=10^{16}、s4=0s_4=0である。定理 3.6 (2)によりei=si−1+xi−sie_i=s_{i-1}+x_i-s_iであるから、e2=e3=1e_2=e_3=1、e4=0e_4=0である。c=fl⁡(fl⁡(1+1)+0)=2c=\operatorname{fl}(\operatorname{fl}(1+1)+0)=2であり、s^=fl⁡(0+2)=2=σ(x)\hat s=\operatorname{fl}(0+2)=2=\sigma(x)である。

例 4.4.F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)、fl⁡\operatorname{fl}を最近接偶数丸め、u=2−53u=2^{-53}とし、x:=(1,2−53,2−106,2−106)x:=(1,2^{-53},2^{-106},2^{-106})とその成分を並べ替えたx′:=(1,2−106,2−106,2−53)x':=(1,2^{-106},2^{-106},2^{-53})を考える。厳密な和はどちらもS:=1+2−53+2−105S:=1+2^{-53}+2^{-105}であり、SSのただ一つの最近接点は1+2−521+2^{-52}である。例 1.5と同じくfl⁡(1+2−53)=1\operatorname{fl}(1+2^{-53})=1である。1+2−1061+2^{-106}のただ一つの最近接点は11であり、fl⁡(1+2−106)=1\operatorname{fl}(1+2^{-106})=1である。2−53+2−1062^{-53}+2^{-106}の最近接点は2−53=252s−532^{-53}=2^{52}s_{-53}と2−53+2−105=(252+1)s−532^{-53}+2^{-105}=(2^{52}+1)s_{-53}であり、整数2522^{52}が偶数であるからfl⁡(2−53+2−106)=2−53\operatorname{fl}(2^{-53}+2^{-106})=2^{-53}である。xxの補償和では、s2=s3=s4=1s_2=s_3=s_4=1、(e2,e3,e4)=(2−53,2−106,2−106)(e_2,e_3,e_4)=(2^{-53},2^{-106},2^{-106})、c=fl⁡(fl⁡(2−53+2−106)+2−106)=2−53c=\operatorname{fl}(\operatorname{fl}(2^{-53}+2^{-106})+2^{-106})=2^{-53}であり、s^=fl⁡(1+2−53)=1\hat s=\operatorname{fl}(1+2^{-53})=1である。誤差S−s^=2−53+2−105S-\hat s=2^{-53}+2^{-105}はu∣S∣=2−53+2−106+2−158u|S|=2^{-53}+2^{-106}+2^{-158}より大きく、定理 4.2 (2)の上界u∣S∣+(1+u)γ2(2−53+2−105)u|S|+(1+u)\gamma_2(2^{-53}+2^{-105})以下である。x′x'の補償和では、s2=s3=s4=1s_2=s_3=s_4=1、(e2,e3,e4)=(2−106,2−106,2−53)(e_2,e_3,e_4)=(2^{-106},2^{-106},2^{-53})、fl⁡(2−106+2−106)=2−105\operatorname{fl}(2^{-106}+2^{-106})=2^{-105}、c=fl⁡(2−105+2−53)=2−53+2−105c=\operatorname{fl}(2^{-105}+2^{-53})=2^{-53}+2^{-105}であり、s^=fl⁡(1+2−53+2−105)=1+2−52\hat s=\operatorname{fl}(1+2^{-53}+2^{-105})=1+2^{-52}である。xxとx′x'の逐次和はどちらも11である。

注意 4.5. 成分の並べ方と加算木を固定すると、逐次和・対分割和・補償和は、それぞれ定義される入力の上でFnF^nからFFへの写像である。例 4.4のxxとx′x'は成分の置換で移り合い、補償和の値が異なるので、補償和は成分の置換で不変な写像ではない。対分割和も、例 2.3の(1016,1,1,−1016)(10^{16},1,1,-10^{16})では00であり、その並べ替え(1016,−1016,1,1)(10^{16},-10^{16},1,1)ではfl⁡(fl⁡(1016−1016)+fl⁡(1+1))=2\operatorname{fl}(\operatorname{fl}(10^{16}-10^{16})+\operatorname{fl}(1+1))=2であるので、成分の置換で不変な写像ではない。定理 4.2 (3)と系 1.4の上界は∑i∣xi∣\sum_i|x_i|と∣∑ixi∣|\sum_ix_i|とnnだけで表されるので成分の置換で不変であり、計算値そのものの置換不変性は与えない。

5 内積

命題 5.1.FFを浮動小数点数系、uuをその単位丸め誤差、fl⁡\operatorname{fl}をFFの最近接丸めとし、n∈N≥1n\in\NN、x,y∈Fnx,y\in F^n、TTを{1,…,n}\{1,\dots,n\}上の加算木とする。整数d≥0d\ge0が(d+1)u<1(d+1)u<1を満たし、すべてのiiについてdT(i)≤dd_T(i)\le dであるとする。各iiについて厳密な積xiyix_iy_iが00であるか正規範囲にあるとし、p^:=(fl⁡(x1y1),…,fl⁡(xnyn))\hat p:=(\operatorname{fl}(x_1y_1),\dots,\operatorname{fl}(x_ny_n))についてTTの各加算の厳密な結果の絶対値がNmax⁡N_{\max}以下であるとする。γd+1:=(d+1)u/(1−(d+1)u)\gamma_{d+1}:=(d+1)u/(1-(d+1)u)と置く。このとき∣θi∣≤γd+1|\theta_i|\le\gamma_{d+1}を満たす実数θ1,…,θn\theta_1,\dots,\theta_nが存在してvT(p^)=∑ixiyi(1+θi)v_T(\hat p)=\sum_ix_iy_i(1+\theta_i)が成り立ち、

∣vT(p^)−∑i=1nxiyi∣≤γd+1∑i=1n∣xiyi∣\Bigl|v_T(\hat p)-\sum_{i=1}^nx_iy_i\Bigr|\le\gamma_{d+1}\sum_{i=1}^n|x_iy_i|

である。特にnu<1nu<1ならば、積を逐次和で加えた計算値vRn(p^)v_{R_n}(\hat p)について、この評価がγn\gamma_nで成り立つ。

証明.xiyi=0x_iy_i=0ならば§E20.1 系 3.2 (1)によりfl⁡(xiyi)=xiyi\operatorname{fl}(x_iy_i)=x_iy_iであり、xiyix_iy_iが正規範囲にあるならば§E20.1 系 3.2 (2)により、いずれの場合も∣εi∣≤u|\varepsilon_i|\le uを満たすεi\varepsilon_iが存在してfl⁡(xiyi)=xiyi(1+εi)\operatorname{fl}(x_iy_i)=x_iy_i(1+\varepsilon_i)である。定理 1.3 (1)によりvT(p^)=∑ixiyi(1+εi)∏j=1dT(i)(1+δi,j)v_T(\hat p)=\sum_ix_iy_i(1+\varepsilon_i)\prod_{j=1}^{d_T(i)}(1+\delta_{i,j})であり、各iiの因子の個数dT(i)+1d_T(i)+1は11以上d+1d+1以下である。§E20.1 補題 4.1をk=dT(i)+1k=d_T(i)+1として適用し、k↦ku/(1−ku)k\mapsto ku/(1-ku)の単調性を用いてθi\theta_iを得る。RnR_nの深さは系 1.4 (1)によりn−1n-1以下であるから、d=n−1d=n-1として最後の主張を得る。▨

例 5.2.F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)、fl⁡\operatorname{fl}を最近接偶数丸め、u=2−53u=2^{-53}とし、x:=(1+2−30,−1)x:=(1+2^{-30},-1)、y:=(1−2−30,1)y:=(1-2^{-30},1)とする。成分はすべて正規範囲のFFの元である。x1y1=1−2−60x_1y_1=1-2^{-60}の最近接点は、11との距離2−602^{-60}が1−2−531-2^{-53}との距離2−53−2−602^{-53}-2^{-60}より小さいので11だけであり、fl⁡(x1y1)=1\operatorname{fl}(x_1y_1)=1である。fl⁡(x2y2)=−1\operatorname{fl}(x_2y_2)=-1であり、逐次和はfl⁡(1−1)=0\operatorname{fl}(1-1)=0であって加算は厳密であるから、厳密な内積−2−60-2^{-60}との誤差はすべて積x1y1x_1y_1の丸めから生じ、相対誤差は11である。∑i∣xiyi∣=2−2−60\sum_i|x_iy_i|=2-2^{-60}であり、命題 5.1の上界はγ2(2−2−60)≈4.44×10−16\gamma_2(2-2^{-60})\approx4.44\times10^{-16}である。c=yc=yとして補題 2.1 (1)を適用すると、yyを固定してxxの関数と見た内積の成分ごとの相対条件数は(2−2−60)/2−60=261−1(2-2^{-60})/2^{-60}=2^{61}-1である。

注意 5.3.F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)、fl⁡\operatorname{fl}を最近接偶数丸め、u=2−53u=2^{-53}とする。x=(Nmax⁡,Nmax⁡,−Nmax⁡)x=(N_{\max},N_{\max},-N_{\max})の逐次和では、最初の加算の厳密な結果2Nmax⁡2N_{\max}がNmax⁡N_{\max}を超え、定理 1.3の仮定は満たされない。CPython の float の計算では、この加算の結果は inf であり、続く加算の結果も inf である。加算木((1),((2),(3)))((1),((2),(3)))ではfl⁡(Nmax⁡−Nmax⁡)=0\operatorname{fl}(N_{\max}-N_{\max})=0、fl⁡(Nmax⁡+0)=Nmax⁡\operatorname{fl}(N_{\max}+0)=N_{\max}であり、厳密な和に等しい。x=y=(2−600)∈F1x=y=(2^{-600})\in F^1では、厳密な積2−12002^{-1200}はアンダーフローの範囲にあり、2−1200<2−1075=s−1022/22^{-1200}<2^{-1075}=s_{-1022}/2であるから、2−12002^{-1200}のただ一つの最近接点は00であり、fl⁡(2−1200)=0\operatorname{fl}(2^{-1200})=0である。相対誤差11は命題 5.1の上界γ1=u/(1−u)\gamma_1=u/(1-u)を超える。

6 Horner 法

定義 6.1.FFを浮動小数点数系、fl⁡\operatorname{fl}をFFの最近接丸めとし、n∈N≥1n\in\NN、a0,…,an∈Fa_0,\dots,a_n\in F、x∈Fx\in Fとする。qn:=anq_n:=a_nと置き、k=n−1,…,0k=n-1,\dots,0の順にqk:=fl⁡(fl⁡(qk+1x)+ak)q_k:=\operatorname{fl}(\operatorname{fl}(q_{k+1}x)+a_k)と置く。q0q_0を、多項式∑i=0naizi\sum_{i=0}^na_iz^iのxxにおける値の Horner 法 (Horner's method) による計算値という。q0q_0は各乗算と各加算の厳密な結果の絶対値がNmax⁡N_{\max}以下であるときに定まる。

定理 6.2.FFを浮動小数点数系、uuをその単位丸め誤差、fl⁡\operatorname{fl}をFFの最近接丸めとし、n∈N≥1n\in\NN、a0,…,an∈Fa_0,\dots,a_n\in F、x∈Fx\in F、p(x):=∑i=0naixip(x):=\sum_{i=0}^na_ix^iとする。2nu<12nu<1とし、γ2n:=2nu/(1−2nu)\gamma_{2n}:=2nu/(1-2nu)と置く。k=n−1,…,0k=n-1,\dots,0の各段で、厳密な積qk+1xq_{k+1}xが00であるか正規範囲にあり、厳密な和fl⁡(qk+1x)+ak\operatorname{fl}(q_{k+1}x)+a_kの絶対値がNmax⁡N_{\max}以下であるとする。

  1. ∣θi∣≤γ2n|\theta_i|\le\gamma_{2n}を満たす実数θ0,…,θn\theta_0,\dots,\theta_nが存在して、Horner 法による計算値はq0=∑i=0nai(1+θi)xiq_0=\sum_{i=0}^na_i(1+\theta_i)x^iを満たす。
  2. ∣q0−p(x)∣≤γ2n∑i=0n∣ai∣∣x∣i|q_0-p(x)|\le\gamma_{2n}\sum_{i=0}^n|a_i||x|^iが成り立つ。Px ⁣:Rn+1→RP_x\colon\R^{n+1}\to\RをPx(a0,…,an):=∑i=0naixiP_x(a_0,\dots,a_n):=\sum_{i=0}^na_ix^iで定めると、p(x)≠0p(x)\ne0ならばκcomp(Px,(a0,…,an))=∑i∣ai∣∣x∣i/∣p(x)∣\kappa_{\mathrm{comp}}(P_x,(a_0,\dots,a_n))=\sum_i|a_i||x|^i/|p(x)|であり、 ∣q0−p(x)∣∣p(x)∣≤γ2n κcomp(Px,(a0,…,an))\frac{|q_0-p(x)|}{|p(x)|}\le\gamma_{2n}\,\kappa_{\mathrm{comp}}(P_x,(a_0,\dots,a_n)) が成り立つ。

証明.0≤k≤n0\le k\le nとk≤i≤nk\le i\le nに対して、mk,n:=2(n−k)m_{k,n}:=2(n-k)、mk,i:=2(i−k)+1m_{k,i}:=2(i-k)+1(i<ni<n)と置く。qn=anq_n=a_nは、mn,n=0m_{n,n}=0であるからqn=∑i=nnaixi−n∏j=1mn,i(1+δn,i,j)q_n=\sum_{i=n}^na_ix^{i-n}\prod_{j=1}^{m_{n,i}}(1+\delta_{n,i,j})の空積の表示をもつ。k<nk<nとし、∣δk+1,i,j∣≤u|\delta_{k+1,i,j}|\le uを満たす実数によりqk+1=∑i=k+1naixi−k−1∏j=1mk+1,i(1+δk+1,i,j)q_{k+1}=\sum_{i=k+1}^na_ix^{i-k-1}\prod_{j=1}^{m_{k+1,i}}(1+\delta_{k+1,i,j})と書けるとする。qk+1xq_{k+1}xが00ならば§E20.1 系 3.2 (1)により、正規範囲にあるならば§E20.1 系 3.2 (2)により、∣μ∣≤u|\mu|\le uを満たすμ\muが存在してfl⁡(qk+1x)=qk+1x(1+μ)\operatorname{fl}(q_{k+1}x)=q_{k+1}x(1+\mu)である。§E20.1 系 3.3により∣α∣≤u|\alpha|\le uを満たすα\alphaが存在して

qk=(qk+1x(1+μ)+ak)(1+α)=qk+1x(1+μ)(1+α)+ak(1+α)q_k=\bigl(q_{k+1}x(1+\mu)+a_k\bigr)(1+\alpha)=q_{k+1}x(1+\mu)(1+\alpha)+a_k(1+\alpha)

である。qk+1q_{k+1}の表示を代入すると、i≥k+1i\ge k+1の項の因子の個数はmk+1,i+2=mk,im_{k+1,i}+2=m_{k,i}、aka_kの因子の個数は1=mk,k1=m_{k,k}であり、qk=∑i=knaixi−k∏j=1mk,i(1+δk,i,j)q_k=\sum_{i=k}^na_ix^{i-k}\prod_{j=1}^{m_{k,i}}(1+\delta_{k,i,j})、∣δk,i,j∣≤u|\delta_{k,i,j}|\le uと書ける。kkの降順の帰納法によりk=0k=0でこの表示が成り立つ。1≤m0,i≤2n1\le m_{0,i}\le2nであるから、§E20.1 補題 4.1を因子の個数m0,im_{0,i}に適用して得るθi\theta_iは∣θi∣≤m0,iu/(1−m0,iu)≤γ2n|\theta_i|\le m_{0,i}u/(1-m_{0,i}u)\le\gamma_{2n}を満たし、(1)が成り立つ。(1)から∣q0−p(x)∣=∣∑iaiθixi∣≤γ2n∑i∣ai∣∣x∣i|q_0-p(x)|=|\sum_ia_i\theta_ix^i|\le\gamma_{2n}\sum_i|a_i||x|^iである。PxP_xはc=(1,x,…,xn)c=(1,x,\dots,x^n)に対する補題 2.1の線形写像であるから、補題 2.1 (1)により条件数の式が成り立ち、両辺を∣p(x)∣|p(x)|で割って最後の不等式を得る。▨

系 6.3.n∈N≥1n\in\NNとし、Rn+2\R^{n+2}にノルム∥⋅∥∞\lVert\cdot\rVert_\infty、R\Rに絶対値を入れ、Hn ⁣:Rn+2→RH_n\colon\R^{n+2}\to\RをHn(a0,…,an,x):=∑i=0naixiH_n(a_0,\dots,a_n,x):=\sum_{i=0}^na_ix^iで定める。Sn\mathfrak S_nを、単位丸め誤差uΦu_\Phiが4nuΦ≤14nu_\Phi\le1を満たす浮動小数点数系Φ\Phiの全体とし、各Φ∈Sn\Phi\in\mathfrak S_nに最近接丸めを一つ固定する。DΦD_\Phiを、(a0,…,an,x)∈Φn+2∖{0}(a_0,\dots,a_n,x)\in\Phi^{n+2}\setminus\{0\}であって定理 6.2の各段の仮定を満たすものの全体とし、F^Φ(a0,…,an,x)\hat F_\Phi(a_0,\dots,a_n,x)を Horner 法による計算値q0q_0とする。このとき、nnを添字とする族は定数(4n)n(4n)_nで後退安定である。

証明.Φ∈Sn\Phi\in\mathfrak S_nと(a,x)=(a0,…,an,x)∈DΦ(a,x)=(a_0,\dots,a_n,x)\in D_\Phiをとる。2nuΦ≤1/2<12nu_\Phi\le1/2<1であるから定理 6.2 (1)が適用され、∣θi∣≤γ2n=2nuΦ/(1−2nuΦ)≤4nuΦ|\theta_i|\le\gamma_{2n}=2nu_\Phi/(1-2nu_\Phi)\le4nu_\Phiである。Δ:=(a0θ0,…,anθn,0)\Delta:=(a_0\theta_0,\dots,a_n\theta_n,0)はHn((a,x)+Δ)=∑iai(1+θi)xi=q0H_n((a,x)+\Delta)=\sum_ia_i(1+\theta_i)x^i=q_0と∥Δ∥∞≤4nuΦmax⁡i∣ai∣≤4nuΦ∥(a,x)∥∞\lVert\Delta\rVert_\infty\le4nu_\Phi\max_i|a_i|\le4nu_\Phi\lVert(a,x)\rVert_\inftyを満たす。したがってq0q_0の相対後退誤差は4nuΦ4nu_\Phi以下であり、§E20.2 定義 5.2 (2)の条件がcn=4nc_n=4nで成り立つ。▨

例 6.4.F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)、fl⁡\operatorname{fl}を最近接偶数丸め、u=2−53u=2^{-53}とし、(z−1)3=z3−3z2+3z−1(z-1)^3=z^3-3z^2+3z-1の係数(a0,a1,a2,a3)=(−1,3,−3,1)(a_0,a_1,a_2,a_3)=(-1,3,-3,1)とx:=1+2−20x:=1+2^{-20}をとる。p(x)=2−60p(x)=2^{-60}である。 Horner 法の各段はq3=1q_3=1、q2=fl⁡(x−3)q_2=\operatorname{fl}(x-3)、q1=fl⁡(fl⁡(q2x)+3)q_1=\operatorname{fl}(\operatorname{fl}(q_2x)+3)である。x−3=−2+2−20x-3=-2+2^{-20}、(−2+2−20)x=−2−2−20+2−40(-2+2^{-20})x=-2-2^{-20}+2^{-40}、−2−2−20+2−40+3=1−2−20+2−40-2-2^{-20}+2^{-40}+3=1-2^{-20}+2^{-40}は§E20.1 補題 1.2 (1)によりFFに属するので、q2=−2+2−20q_2=-2+2^{-20}、q1=1−2−20+2−40q_1=1-2^{-20}+2^{-40}である。q1x=1+2−60q_1x=1+2^{-60}の最近接点は11だけであり、q0=fl⁡(1−1)=0q_0=\operatorname{fl}(1-1)=0である。すべての積は正規範囲にあり、誤差は積q1xq_1xの丸め一回から生じる。∑i∣ai∣∣x∣i=(1+x)3=(2+2−20)3\sum_i|a_i||x|^i=(1+x)^3=(2+2^{-20})^3であるから、補題 2.1 (2)によりq0q_0の成分ごとの相対後退誤差は2−60/(2+2−20)3≈1.08×10−19≈9.8×10−4u2^{-60}/(2+2^{-20})^3\approx1.08\times10^{-19}\approx9.8\times10^{-4}uであり、定理 6.2の上界γ6≈6u\gamma_6\approx6u以下である。一方、相対前進誤差は11である。定理 6.2 (2)の条件数は(2+2−20)3/2−60≈9.22×1018(2+2^{-20})^3/2^{-60}\approx9.22\times10^{18}であり、相対前進誤差は条件数と成分ごとの相対後退誤差の積に等しい。xxが根11に近づくとp(x)→0p(x)\to0、∑i∣ai∣∣x∣i→8\sum_i|a_i||x|^i\to8であるから、条件数は+∞+\inftyに発散する。根11とxxからd:=fl⁡(x−1)d:=\operatorname{fl}(x-1)、fl⁡(fl⁡(d⋅d)⋅d)\operatorname{fl}(\operatorname{fl}(d\cdot d)\cdot d)を計算すると、d=2−20d=2^{-20}と二回の乗算fl⁡(d⋅d)=2−40\operatorname{fl}(d\cdot d)=2^{-40}、fl⁡(2−40d)=2−60\operatorname{fl}(2^{-40}d)=2^{-60}はすべて厳密であり、計算値は2−60=p(x)2^{-60}=p(x)に等しい。

7 二次方程式の小さい根

命題 7.1.B>2B>2を実数とし、U:={(a,b,c)∈R3∣a>0, b>0, c>0, b2−4ac>0}U:=\{(a,b,c)\in\R^3\mid a>0,\ b>0,\ c>0,\ b^2-4ac>0\}、R ⁣:U→RR\colon U\to\RをR(a,b,c):=(−b+b2−4ac)/(2a)R(a,b,c):=(-b+\sqrt{b^2-4ac})/(2a)で定める。r:=R(1,B,1)r:=R(1,B,1)と置く。

  1. (a,b,c)∈U(a,b,c)\in Uならば、az2+bz+caz^2+bz+cの二つの根はR(a,b,c)R(a,b,c)とc/(aR(a,b,c))c/(aR(a,b,c))であり、どちらも負であってR(a,b,c)>c/(aR(a,b,c))R(a,b,c)>c/(aR(a,b,c))である。z^<0\hat z<0がaz2+bz+caz^2+bz+cの根でありz^2<c/a\hat z^2<c/aを満たすならば、z^=R(a,b,c)\hat z=R(a,b,c)である。
  2. κcomp(R,(1,B,1))=2B/B2−4\kappa_{\mathrm{comp}}(R,(1,B,1))=2B/\sqrt{B^2-4}である。
  3. −1<z^<0-1<\hat z<0ならばηcomp(R,(1,B,1),z^)=∣z^2+Bz^+1∣/(z^2+B∣z^∣+1)\eta_{\mathrm{comp}}(R,(1,B,1),\hat z)=|\hat z^2+B\hat z+1|/(\hat z^2+B|\hat z|+1)であり、この下限は最小値である。

証明.(1)を示す。D:=b2−4ac>0D:=b^2-4ac>0であるから、az2+bz+caz^2+bz+cの根はz±:=(−b±D)/(2a)z_\pm:=(-b\pm\sqrt D)/(2a)の二つであり、z+=R(a,b,c)>z−z_+=R(a,b,c)>z_-である。z+z−=c/a>0z_+z_-=c/a>0かつz++z−=−b/a<0z_++z_-=-b/a<0であるから、二つの根はどちらも負であり、z−=c/(az+)z_-=c/(az_+)である。z^<0\hat z<0が根でありz^≠z+\hat z\ne z_+ならばz^=z−<z+<0\hat z=z_-<z_+<0であるから∣z^∣>∣z+∣|\hat z|>|z_+|であり、z^2>∣z^∣∣z+∣=c/a\hat z^2>|\hat z||z_+|=c/aである。

(2)を示す。x0:=(1,B,1)x_0:=(1,B,1)と置き、§E20.2 定義 5.4のUx0U_{x_0}とG:=Gx0G:=G_{x_0}を用いる。RRは多項式、(0,∞)(0,\infty)上の平方根、a>0a>0での除算の合成であるからUU上でC1C^1級であり、G(δ)=R(x0+x0∘δ)G(\delta)=R(x_0+x_0\circ\delta)は00を含む開集合Ux0U_{x_0}上で全微分可能である。(1)により、任意のδ∈Ux0\delta\in U_{x_0}に対して(1+δ1)G(δ)2+B(1+δ2)G(δ)+(1+δ3)=0(1+\delta_1)G(\delta)^2+B(1+\delta_2)G(\delta)+(1+\delta_3)=0が成り立つ。この恒等式をδ=0\delta=0で微分し、G(0)=rG(0)=rを用いると、任意のh∈R3h\in\R^3に対して(2r+B) DG(0)h+r2h1+Brh2+h3=0(2r+B)\,DG(0)h+r^2h_1+Brh_2+h_3=0を得る。2r+B=B2−4≠02r+B=\sqrt{B^2-4}\ne0であるからDG(0)h=−(r2h1+Brh2+h3)/B2−4DG(0)h=-(r^2h_1+Brh_2+h_3)/\sqrt{B^2-4}である。∥h∥∞=1\lVert h\rVert_\infty=1ならば∣r2h1+Brh2+h3∣≤r2+B∣r∣+1|r^2h_1+Brh_2+h_3|\le r^2+B|r|+1であり、(1)によりr<0r<0であるから、h=(1,−1,1)h=(1,-1,1)で等号が成り立つ。したがって∥DG(0)∥=(r2+B∣r∣+1)/B2−4\lVert DG(0)\rVert=(r^2+B|r|+1)/\sqrt{B^2-4}であり、§E20.2 命題 3.2 (2)によりκabs(G,0)=(r2+B∣r∣+1)/B2−4\kappa_{\mathrm{abs}}(G,0)=(r^2+B|r|+1)/\sqrt{B^2-4}である。r2+1=−Br=B∣r∣r^2+1=-Br=B|r|であるからκabs(G,0)=2B∣r∣/B2−4\kappa_{\mathrm{abs}}(G,0)=2B|r|/\sqrt{B^2-4}であり、κcomp(R,x0)=κabs(G,0)/∣r∣=2B/B2−4\kappa_{\mathrm{comp}}(R,x_0)=\kappa_{\mathrm{abs}}(G,0)/|r|=2B/\sqrt{B^2-4}である。

(3)を示す。−1<z^<0-1<\hat z<0とし、L(a,b,c):=az^2+bz^+cL(a,b,c):=a\hat z^2+b\hat z+c、M:=z^2+B∣z^∣+1M:=\hat z^2+B|\hat z|+1、η∗:=∣L(x0)∣/M\eta^*:=|L(x_0)|/Mと置く。L(x0)=z^2+Bz^+1L(x_0)=\hat z^2+B\hat z+1である。δ∈Ux0\delta\in U_{x_0}がR(x0+x0∘δ)=z^R(x_0+x_0\circ\delta)=\hat zを満たすならばL(x0+x0∘δ)=0L(x_0+x_0\circ\delta)=0であるから、補題 2.1 (2)により∥δ∥∞≥η∗\lVert\delta\rVert_\infty\ge\eta^*である。(z^2,Bz^,1)(\hat z^2,B\hat z,1)の符号は(1,−1,1)(1,-1,1)であるから、補題 2.1 (2)のy^=0\hat y=0に対するδ∗\delta^*は、τ:=L(x0)/M\tau:=L(x_0)/Mと置いてδ∗=(−τ,τ,−τ)\delta^*=(-\tau,\tau,-\tau)であり、∥δ∗∥∞=η∗\lVert\delta^*\rVert_\infty=\eta^*とL(x0+x0∘δ∗)=0L(x_0+x_0\circ\delta^*)=0を満たす。z^≠0\hat z\ne0であるから∣L(x0)∣=∣z^2+1−B∣z^∣∣<M|L(x_0)|=|\hat z^2+1-B|\hat z||<Mであり、∣τ∣<1|\tau|<1である。(a′,b′,c′):=(1−τ,B(1+τ),1−τ)(a',b',c'):=(1-\tau,B(1+\tau),1-\tau)の成分は正であり、z^\hat zはa′z2+b′z+c′a'z^2+b'z+c'の実根であるからb′2−4a′c′≥0b'^2-4a'c'\ge0である。b′2−4a′c′=0b'^2-4a'c'=0ならばz^\hat zは重根であってz^2=c′/a′=1\hat z^2=c'/a'=1となり、∣z^∣<1|\hat z|<1に反する。よって(a′,b′,c′)∈U(a',b',c')\in Uであり、z^2<1=c′/a′\hat z^2<1=c'/a'であるから、(1)によりR(a′,b′,c′)=z^R(a',b',c')=\hat zである。したがってδ∗∈Ux0\delta^*\in U_{x_0}は下限を達成する。▨

例 7.2.F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)、fl⁡\operatorname{fl}を最近接偶数丸め、u=2−53u=2^{-53}とし、B:=108B:=10^8とする。z2+Bz+1z^2+Bz+1の二つの根をr:=(−B+w)/2r:=(-B+w)/2、ℓ:=(−B−w)/2\ell:=(-B-w)/2(w:=B2−4w:=\sqrt{B^2-4})とすると、根と係数の関係によりrℓ=1r\ell=1である。平方根も正しく丸めてw^:=fl⁡(fl⁡(fl⁡(B⋅B)−4))\hat w:=\operatorname{fl}(\sqrt{\operatorname{fl}(\operatorname{fl}(B\cdot B)-4)})と置き、小さい根rrの二つの計算値を

r^D:=fl⁡(fl⁡(−B+w^)/2),r^V:=fl⁡(1/ℓ^),ℓ^:=fl⁡(fl⁡(−B−w^)/2)\hat r_{\mathrm D}:=\operatorname{fl}\bigl(\operatorname{fl}(-B+\hat w)/2\bigr),\qquad\hat r_{\mathrm V}:=\operatorname{fl}(1/\hat\ell),\quad\hat\ell:=\operatorname{fl}\bigl(\operatorname{fl}(-B-\hat w)/2\bigr)

で定める。命題 7.1 (2)により、係数(1,B,1)(1,B,1)から小さい根への写像RRの成分ごとの相対条件数は2B/B2−4=2/1−4⋅10−162B/\sqrt{B^2-4}=2/\sqrt{1-4\cdot10^{-16}}であり、22と2+10−152+10^{-15}の間にある。r^D\hat r_{\mathrm D}の分子は、写像ω↦ω−B\omega\mapsto\omega-Bにwwの近似値w^\hat wを入れた値である。§E20.2 例 3.4 (3)により、この写像のwwにおける相対条件数はw/(B−w)=w(B+w)/4w/(B-w)=w(B+w)/4であり、4.99×10154.99\times10^{15}より大きい。ℓ^\hat\ellの分子の絶対値は、写像ω↦ω+B\omega\mapsto\omega+Bにw^\hat wを入れた値である。同じ例をa=−Ba=-Bとして、この写像のwwにおける相対条件数はw/(w+B)<1w/(w+B)<1である。ω↦1/ω\omega\mapsto1/\omegaの相対条件数は§E20.2 系 3.3により11である。

注意 7.3.例 7.2において、101610^{16}と1016−410^{16}-4は[253,254)[2^{53},2^{54})に属する偶数であるからFFに属し、fl⁡(B⋅B)=1016\operatorname{fl}(B\cdot B)=10^{16}、fl⁡(1016−4)=1016−4\operatorname{fl}(10^{16}-4)=10^{16}-4である。B−w=4/(B+w)B-w=4/(B+w)は2⋅10−82\cdot10^{-8}より大きく2.0000001⋅10−82.0000001\cdot10^{-8}より小さい。F∩[226,227)F\cap[2^{26},2^{27})はs26=2−26s_{26}=2^{-26}の倍数からなり、wwに近い元BB、B−2−26B-2^{-26}、B−2−25B-2^{-25}とwwとの距離はそれぞれ約2.0×10−82.0\times10^{-8}、5.1×10−95.1\times10^{-9}、9.8×10−99.8\times10^{-9}であるから、w^=B−2−26\hat w=B-2^{-26}である。w^\hat wの相対誤差は約0.459u0.459uである。w^/2≤B≤2w^\hat w/2\le B\le2\hat wであるから、補題 3.1によりfl⁡(−B+w^)=−2−26\operatorname{fl}(-B+\hat w)=-2^{-26}であり、−2−27∈F-2^{-27}\in Fであるからr^D=−2−27≈−7.450580596923828×10−9\hat r_{\mathrm D}=-2^{-27}\approx-7.450580596923828\times10^{-9}である。r^D=(w^−B)/2\hat r_{\mathrm D}=(\hat w-B)/2とr=(w−B)/2r=(w-B)/2から

r^D−rr=w^−ww⋅ww−B\frac{\hat r_{\mathrm D}-r}{r}=\frac{\hat w-w}{w}\cdot\frac{w}{w-B}

であり、減算は厳密に行われ、相対誤差の絶対値はw^\hat wの相対誤差の絶対値のw/(B−w)w/(B-w)倍である。−B−w^=−(2⋅108−2−26)-B-\hat w=-(2\cdot10^8-2^{-26})は、s27=2−25s_{27}=2^{-25}の倍数からなるF∩(−228,−227]F\cap(-2^{28},-2^{27}]の隣接する二元−(2⋅108−2−25)-(2\cdot10^8-2^{-25})と−2⋅108-2\cdot10^8の中点である。2⋅108=(2⋅108⋅225)s272\cdot10^8=(2\cdot10^8\cdot2^{25})s_{27}の整数2⋅108⋅2252\cdot10^8\cdot2^{25}は偶数であり、他方の元の整数はそれより11小さい奇数であるからfl⁡(−B−w^)=−2⋅108\operatorname{fl}(-B-\hat w)=-2\cdot10^8であり、ℓ^=−108\hat\ell=-10^8、r^V=fl⁡(−10−8)≈−1.00000000000000002×10−8\hat r_{\mathrm V}=\operatorname{fl}(-10^{-8})\approx-1.00000000000000002\times10^{-8}である。参照値r=−2/(B+w)=−1.0000000000000001000000000000000200…×10−8r=-2/(B+w)=-1.0000000000000001000000000000000200\ldots\times10^{-8}を 60 桁の十進演算で求めると、相対前進誤差はr^D\hat r_{\mathrm D}で0.25494194030761726…0.25494194030761726\ldots、r^V\hat r_{\mathrm V}で7.9077439169871539…×10−17≈0.712u7.9077439169871539\ldots\times10^{-17}\approx0.712uである。CPython の float と math.sqrt による計算値はr^D\hat r_{\mathrm D}、r^V\hat r_{\mathrm V}と一致する。

例 7.4.注意 7.3のr^D=−2−27\hat r_{\mathrm D}=-2^{-27}とr^V\hat r_{\mathrm V}はどちらも(−1,0)(-1,0)に属するので、命題 7.1 (3)により成分ごとの相対後退誤差は∣z^2+Bz^+1∣/(z^2+B∣z^∣+1)|\hat z^2+B\hat z+1|/(\hat z^2+B|\hat z|+1)である。r^D\hat r_{\mathrm D}では分子が1−108⋅2−27+2−54≈0.25494194030761721-10^8\cdot2^{-27}+2^{-54}\approx0.2549419403076172、分母が1+108⋅2−27+2−54≈1.74505805969238281+10^8\cdot2^{-27}+2^{-54}\approx1.7450580596923828であり、後退誤差は約0.146093672294524200.14609367229452420である。r^V\hat r_{\mathrm V}の後退誤差は約3.9538719584935763×10−17≈0.356u3.9538719584935763\times10^{-17}\approx0.356uである。条件数κ:=2B/B2−4\kappa:=2B/\sqrt{B^2-4}と後退誤差η\etaの積は、r^D\hat r_{\mathrm D}で約0.29220.2922、r^V\hat r_{\mathrm V}で約0.712u0.712uであり、注意 7.3の相対前進誤差0.25490.2549と0.712u0.712uと同じ大きさである。

前提記事