1 加算木と和の誤差
定義 1.1.n∈N≥1とする。{1,…,n}の空でない部分集合I上の 加算木 (summation tree) を、Iの元の個数に関して帰納的に次のように定める。
- I={i}のとき、記号(i)だけをI上の加算木とする。
- Iが二つ以上の元をもつとき、Iを交わらない空でない二つの部分集合I1,I2の和集合に分け、T1をI1上の加算木、T2をI2上の加算木として、順序対T=(T1,T2)をI上の加算木とする。
I上の加算木Tとi∈Iに対して、T=(i)のときdT(i):=0と置き、T=(T1,T2)かつi∈IkのときdT(i):=1+dTk(i)と置く。dT(i)をTにおけるiの 深さ (depth) という。
Fを浮動小数点数系、flをFの最近接丸めとし、x=(x1,…,xn)∈Fnとする。T=(i)のときvT(x):=xiと置き、T=(T1,T2)のときvT(x):=fl(vT1(x)+vT2(x))と置く。vT(x)は、この帰納的な定義に現れる各実数vT1(x)+vT2(x)(Tの加算の厳密な結果)の絶対値がNmax以下であるときに定まる。vT(x)をTによるxの和の計算値という。
R1:=(1)、Rk:=(Rk−1,(k))(2≤k≤n)と置き、vRn(x)をxの 逐次和 (recursive summation) という。整数j≥1とl≥0に対して、Pj,0:=(j)、Pj,l:=(P2j−1,l−1,P2j,l−1)(l≥1)と置く。lに関する帰納法により、n≥j2lのときPj,lは{(j−1)2l+1,…,j2l}上の加算木である。m:=⌈log2n⌉と置き、xに零を補ったx′:=(x1,…,xn,0,…,0)∈F2mに対するvP1,m(x′)をxの 対分割和 (pairwise summation) という。
例 1.2.F=F(2,53,−1022,1023)、flを最近接偶数丸めとし、x=(1016,−1016,1)とする。253≤1016<254であり、§E20.1 補題 1.2 (1)によりF∩[253,254)はs53=2の倍数である偶数の全体であるから、±1016∈Fである。逐次和はvR3(x)=fl(fl(1016−1016)+1)=fl(0+1)=1であり、厳密な和に等しい。加算木T:=((1),((2),(3)))では、−1016+1=−9999999999999999の最近接点は−1016と−9999999999999998の二つである。1016=(5⋅1015)s53、9999999999999998=4999999999999999s53であり、整数5⋅1015は偶数、4999999999999999は奇数であるから、fl(−1016+1)=−1016であり、vT(x)=fl(1016−1016)=0である。−1016+1に対して−9999999999999998を選ぶ最近接丸めでは、vT(x)=fl(1016−9999999999999998)=2である。
定理 1.3.Fを浮動小数点数系、uをその単位丸め誤差、flをFの最近接丸めとし、n∈N≥1、x=(x1,…,xn)∈Fn、Tを{1,…,n}上の加算木とする。Tの各加算の厳密な結果の絶対値がNmax以下であるとする。
- 各iに対して、∣δi,j∣≤uを満たす実数δi,1,…,δi,dT(i)が存在して
vT(x)=i=1∑nxij=1∏dT(i)(1+δi,j)
が成り立つ。
- 整数d≥0がdu<1を満たし、すべてのiについてdT(i)≤dであるとする。0≤k≤dに対してγk:=ku/(1−ku)と置く。このとき∣θi∣≤γdT(i)≤γdを満たす実数θ1,…,θnが存在してvT(x)=∑i=1nxi(1+θi)が成り立ち、
vT(x)−i=1∑nxi≤γdi=1∑n∣xi∣
である。
証明.Iを{1,…,n}の空でない部分集合、TをI上の加算木とする。T=(i)ならばvT(x)=xi、dT(i)=0であり、vT(x)=∑i∈Ixi∏j=1dT(i)(1+δi,j)は空積の表示として成り立つ。T=(T1,T2)とし、k=1,2について、∣δi,j∣≤uを満たす実数によりvTk(x)=∑i∈Ikxi∏j=1dTk(i)(1+δi,j)と書けるとする。§E20.1 系 3.3により∣δ∣≤uを満たすδが存在して
vT(x)=(vT1(x)+vT2(x))(1+δ)=k=1∑2i∈Ik∑xi(1+δ)j=1∏dTk(i)(1+δi,j)であり、i∈Ikに対してdT(i)=1+dTk(i)であるから、vT(x)も各iの因子の個数がdT(i)である同じ形の表示をもつ。I1、I2の元の個数はIより少ないので、Iの元の個数に関する帰納法により任意の加算木がこの表示をもち、I={1,…,n}として(1)を得る。
(2)を示す。dT(i)≥1であるiについては、dT(i)u≤du<1であるから、§E20.1 補題 4.1をk=dT(i)、εj=1として適用して、∏j=1dT(i)(1+δi,j)=1+θi、∣θi∣≤γdT(i)を満たすθiを得る。dT(i)=0であるiについてはθi:=0=γ0と置く。k↦ku/(1−ku)は0≤k≤dで単調非減少であるからγdT(i)≤γdであり、∣vT(x)−∑ixi∣=∣∑ixiθi∣≤γd∑i∣xi∣が成り立つ。▨
系 1.4.Fを浮動小数点数系、uをその単位丸め誤差、flをFの最近接丸めとし、n∈N≥1、x∈Fnとする。整数k≥0がku<1を満たすときγk:=ku/(1−ku)と置く。
- dRn(1)=n−1であり、2≤i≤nに対してdRn(i)=n−i+1である。(n−1)u<1であり、逐次和の各加算の厳密な結果の絶対値がNmax以下であるならば、∣vRn(x)−∑ixi∣≤γn−1∑i∣xi∣が成り立つ。
- m:=⌈log2n⌉と置く。P1,mにおけるすべての元の深さはmである。mu<1であり、対分割和の各加算の厳密な結果の絶対値がNmax以下であるならば、対分割和s^は∣s^−∑ixi∣≤γm∑i∣xi∣を満たす。
証明.R1=(1)ではdR1(1)=0である。n≥2ならばRn=(Rn−1,(n))であるから、dRn(n)=1であり、i≤n−1に対してdRn(i)=1+dRn−1(i)である。nに関する帰納法により深さの式が成り立ち、その最大値はn−1である。Pj,0=(j)の元の深さは0であり、Pj,l=(P2j−1,l−1,P2j,l−1)の元の深さはP2j−1,l−1またはP2j,l−1における深さに1を加えたものであるから、lに関する帰納法によりPj,lの元の深さはすべてlである。零を補ったx′は∑ixi′=∑ixiと∑i∣xi′∣=∑i∣xi∣を満たす。定理 1.3 (2)を、逐次和にはd=n−1で、対分割和にはx′とd=mで適用する。▨
例 1.5.F=F(2,53,−1022,1023)、flを最近接偶数丸め、u=2−53とし、x:=(1,2−53,…,2−53)∈F8(2−53が七つ)とする。厳密な和は1+7uであり、全成分が正であるから∑i∣xi∣=1+7uである。1+2−53の最近接点は1=252s0と1+2−52=(252+1)s0であり、整数252が偶数であるからfl(1+2−53)=1である。したがって逐次和の各段の値は1であり、vR8(x)=1、絶対誤差は7u≈7.77×10−16である。系 1.4 (1)の上界はγ7(1+7u)=7u(1+7u)/(1−7u)であり、誤差との比は(1−7u)/(1+7u)>1−14uである。対分割和ではm=3であり、第一段の値はfl(1+2−53)=1と三つの2−53+2−53=2−52、第二段の値は1+2−52と2−51、第三段の値は1+2−52+2−51=1+3⋅2−52∈Fである。絶対誤差はu≈1.11×10−16であり、系 1.4 (2)の上界γ3(1+7u)≈3.33×10−16以下である。
2 和の条件数
補題 2.1.n∈N≥1、c∈Rnとし、L:Rn→RをL(x):=∑i=1ncixiで定める。x∈Rnに対してM:=∑i=1n∣cixi∣と置き、sgn(0):=0とし、δ∈Rnに対するx∘δを§E20.2 定義 5.4の成分ごとの積とする。
- L(x)=0ならばκcomp(L,x)=M/∣L(x)∣である。
- M>0ならば、任意のy^∈Rに対してηcomp(L,x,y^)=∣y^−L(x)∣/Mである。さらにδi:=sgn(cixi)(y^−L(x))/Mで定まるδ∈Rnは∥δ∥∞=∣y^−L(x)∣/MとL(x+x∘δ)=y^を満たす。
証明.Lの定義域はRnであるから§E20.2 定義 5.4のUxはRnであり、任意のδ∈Rnに対して
Gx(δ)−Gx(0)=i=1∑ncixiδi,∣Gx(δ)−Gx(0)∣≤M∥δ∥∞が成り立つ。
(1)を示す。M≥∣L(x)∣>0である。ε>0に対して、§E20.2 定義 3.1の上限s(ε)は上の不等式によりM以下である。δ:=ε(sgn(cixi))iは∥δ∥∞=εとGx(δ)−Gx(0)=εMを満たすから、s(ε)=Mである。したがってκabs(Gx,0)=Mである。
(2)を示す。L(x+x∘δ)=y^ならば、∣y^−L(x)∣=∣Gx(δ)−Gx(0)∣≤M∥δ∥∞であるから、∥δ∥∞≥∣y^−L(x)∣/Mである。主張のδは∑icixiδi=∑i∣cixi∣(y^−L(x))/M=y^−L(x)を満たし、M>0であるからcixi=0であるiが存在して∥δ∥∞=∣y^−L(x)∣/Mである。▨
系 2.2.n∈N≥1とし、σ:Rn→Rをσ(x):=∑i=1nxiで定める。
- σ(x)=0ならばκcomp(σ,x)=∑i∣xi∣/∣∑ixi∣である。
- F、u、fl、x∈Fn、T、d、γdが定理 1.3と定理 1.3 (2)の仮定を満たすならば、ηcomp(σ,x,vT(x))≤γdである。さらにσ(x)=0ならば
∣σ(x)∣∣vT(x)−σ(x)∣=κcomp(σ,x)ηcomp(σ,x,vT(x))≤γdκcomp(σ,x)
が成り立つ。
証明.(1)は補題 2.1 (1)をc=(1,…,1)として適用したものである。定理 1.3 (2)のθ=(θ1,…,θn)はσ(x+x∘θ)=vT(x)と∥θ∥∞≤γdを満たすので、ηcomp(σ,x,vT(x))≤γdである。σ(x)=0ならば∑i∣xi∣>0であり、補題 2.1 (2)によりηcomp(σ,x,vT(x))=∣vT(x)−σ(x)∣/∑i∣xi∣であるから、(1)と合わせて等式が成り立つ。▨
例 2.3.F=F(2,53,−1022,1023)、flを最近接偶数丸め、u=2−53とし、x:=(1016,1,1,−1016)とする。σ(x)=2、∑i∣xi∣=2⋅1016+2であり、系 2.2 (1)によりκcomp(σ,x)=1016+1である。1016+1の最近接点は1016と1016+2であり、1016=(5⋅1015)s53、1016+2=(5⋅1015+1)s53であり、整数5⋅1015が偶数であるからfl(1016+1)=1016である。したがって逐次和はfl(fl(1016+1)+1)=1016を経てvR4(x)=0である。対分割和は、例 1.2と同じくfl(1−1016)=−1016であるから、fl(1016−1016)=0である。どちらの計算値も相対誤差は1である。系 2.2 (2)の上界は、逐次和でγ3κcomp(σ,x)≈3.33、対分割和でγ2κcomp(σ,x)≈2.22であり、絶対誤差の上界はγ3(2⋅1016+2)≈6.66とγ2(2⋅1016+2)≈4.44である。
3 丸め誤差の回収
補題 3.1.F=F(β,p,emin,emax)を浮動小数点数系(非正規数を含む)、flをFの最近接丸めとする。x,y∈Fがy/2≤x≤2yを満たすならば、x−y∈Fであり、fl(x−y)=x−yである。
証明.y/2≤2yからy≥0であり、x≥y/2≥0である。条件y/2≤x≤2yはx/2≤y≤2xと同値であり、F=−Fであるから、x≥yとしてよい。x≤2yから0≤x−y≤yである。y=0ならばx=0であり、x−y=0∈Fである。y>0とし、§E20.1 補題 1.2の整数e(y)をeと置く。§E20.1 補題 1.3 (3)によりe(x)≥eであり、§E20.1 補題 1.3 (2)によりxとyはseの整数倍であって、y=mse、1≤m≤βp−1と書ける。したがってx−y=kseを満たす整数kがあり、0≤x−y≤yから0≤k≤m≤βp−1である。∣x−y∣≤y≤Nmaxであるから、§E20.1 補題 1.3 (1)によりx−y∈Fであり、§E20.1 系 3.2 (1)によりfl(x−y)=x−yである。▨
補題 3.3.F=F(β,p,emin,emax)を浮動小数点数系とし、a,b∈Fが∣a+b∣≤Nmaxを満たすとする。s∈Fがa+bの最近接点ならば、r:=a+b−sはFに属し、∣r∣≤min{∣a∣,∣b∣}を満たす。
証明.t:=a+bと置く。a,b∈Fであるから∣r∣=∣t−s∣≤∣t−a∣=∣b∣かつ∣r∣≤∣t−b∣=∣a∣である。主張はaとbについて対称であるから、∣b∣≤∣a∣としてよい。b=0ならばt=a∈Fであり、tとの距離が0であるFの元はtだけであるからs=t、r=0である。b=0とし、e:=e(∣b∣)と置く。§E20.1 補題 1.3 (3)によりe(∣a∣)≥eであり、§E20.1 補題 1.3 (2)によりaとbはseの整数倍であって、∣b∣=mse、1≤m≤βp−1と書ける。したがってtもseの整数倍である。sがseの整数倍でないと仮定する。このときs=0であり、§E20.1 補題 1.3 (2)によりe(∣s∣)<eである。e(⋅)の定義により∣s∣<βe(∣s∣)+1≤βeである。∣t∣≥βeならば、v:=sgn(t)βe=sgn(t)βp−1seはFの元であり、∣t−v∣=∣t∣−βe<∣t∣−∣s∣≤∣t−s∣を満たすので、sが最近接点であることに反する。よって∣t∣<βe=βp−1seであり、t=kse、∣k∣<βp−1を満たす整数kについて§E20.1 補題 1.3 (1)によりt∈Fである。このときs=tはseの整数倍であり、仮定に反する。したがってsはseの整数倍であり、r=kseを満たす整数kがある。∣r∣≤∣b∣から∣k∣≤m≤βp−1であり、∣r∣≤∣b∣≤Nmaxであるから、§E20.1 補題 1.3 (1)によりr∈Fである。▨
補題 3.4.F=F(2,p,emin,emax)を基数2の浮動小数点数系とし、x,y∈Fが∣x∣≥∣y∣と∣x+y∣≤Nmaxを満たすとする。s∈Fをx+yの最近接点とする。
- s−x∈Fである。
- x+y∈/Fならば∣s∣≥∣y∣である。
証明.t:=x+yと置く。F=−Fであり、−sは−tの最近接点であるから、(x,y,s)を(−x,−y,−s)に替えても仮定と結論は変わらない。したがってx≥0としてよい。t∈Fならば、tとの距離が0であるFの元はtだけであるからs=tであり、s−x=y∈Fである。このとき(2)の仮定は成り立たない。以下t∈/Fとする。このときy=0であり、x≥∣y∣>0である。v∈Fがv≤tを満たすならばv≤sである。実際、s<v≤tならば∣v−t∣<∣s−t∣となり、sが最近接点であることに反する。同様に、v∈Fがv≥tを満たすならばv≥sである。§E20.1 補題 1.3 (2)によりx=mse、e:=e(x)、1≤m≤2p−1と書く。
y>0の場合、x<t≤2xであり、v=xとしてs≥xである。2x≤Nmaxならば2x∈Fである。実際、m≤2p−1ならば2x=(2m)se、2m≤2pであり、m>2p−1ならばx≥2p−1se=2eかつ2x=mse+1であって、2e+1≤2x≤Nmax<2emax+1からe+1≤emaxであるので、どちらの場合も§E20.1 補題 1.3 (1)が適用される。このときv=2xとしてs≤2xである。2x>Nmaxならばs≤Nmax<2xである。したがってx/2≤x≤s≤2xであり、補題 3.1によりs−x∈Fである。また∣s∣=s≥x≥∣y∣である。
y<0の場合、0≤t<xである。−y≥x/2ならば、(−y)/2≤x≤2(−y)であるから補題 3.1によりt=x−(−y)∈Fとなり、t∈/Fに反する。よって−y<x/2であり、x/2<t<xである。e=eminならば、§E20.1 補題 1.3 (2)によりyもseminの整数倍であるからt=kseminを満たす整数kがあり、0<t<x=mseminから0<k<m≤2p−1であって、§E20.1 補題 1.3 (1)によりt∈Fとなり、t∈/Fに反する。よってe>eminであり、x/2=mse−1は§E20.1 補題 1.3 (1)によりFに属する。v=x/2とv=xとしてx/2≤s≤xであり、補題 3.1によりs−x∈Fである。また∣s∣=s≥x/2>−y=∣y∣である。▨
定義 3.5.Fを浮動小数点数系、flをFの最近接丸めとし、a,b∈Fとする。
s:=fl(a+b),z:=fl(s−a),a′:=fl(s−z),e:=fl(fl(a−a′)+fl(b−z))と置く。これら六つの演算の厳密な結果の絶対値がすべてNmax以下であるとき、対(s,e)を(a,b)の TwoSum (TwoSum) という。
定理 3.6.p≥1、emin≤emaxを整数とし、F=F(2,p,emin,emax)を基数2の浮動小数点数系、uをその単位丸め誤差、flをFの最近接丸めとする。a,b∈Fが∣a+b∣≤Nmaxと∣fl(a+b)−a∣≤Nmaxを満たすとし、s、z、a′、eを定義 3.5の値とする。
- s−z、a−a′、b−z、fl(a−a′)+fl(b−z)はすべてFに属する。したがって(a,b)の TwoSum(s,e)が定まり、a′=s−z、fl(a−a′)=a−a′、fl(b−z)=b−z、e=(a−a′)+(b−z)が成り立つ。
- a+b=s+eが成り立つ。
- ∣e∣≤min{∣a∣,∣b∣}かつ∣e∣≤u∣a+b∣である。
証明.t:=a+b、r:=t−sと置く。s=fl(t)はtの最近接点であるから、補題 3.3によりr∈F、∣r∣≤min{∣a∣,∣b∣}である。
∣a∣≥∣b∣の場合。補題 3.4 (1)を(x,y)=(a,b)に適用してs−a∈Fであるから、z=s−aである。s−z=a∈Fであるからa′=aであり、a−a′=0∈F、b−z=t−s=r∈Fである。(a−a′)+(b−z)=rである。
∣a∣<∣b∣かつt∈Fの場合。s=tであるからs−a=b∈Fであり、z=bである。s−z=a∈Fであるからa′=aであり、a−a′=0、b−z=0はともにFに属し、(a−a′)+(b−z)=0=rである。
∣a∣<∣b∣かつt∈/Fの場合。補題 3.4 (2)を(x,y)=(b,a)に適用して∣s∣≥∣a∣である。z=fl(s+(−a))はs+(−a)の最近接点であり、s,−a∈F、∣s∣≥∣−a∣、∣s−a∣≤Nmaxであるから、補題 3.4 (1)を(x,y)=(s,−a)に適用してz−s∈Fである。よってs−z∈Fであり、a′=s−zである。補題 3.3をsと−aに適用してρ:=(s−a)−z∈Fであり、a−a′=a−s+z=−ρ∈Fである。実数としてs−a=b−rであるから、zはb+(−r)の最近接点でもある。b,−r∈F、∣−r∣≤∣a∣<∣b∣、∣b−r∣=∣s−a∣≤Nmaxであるから、補題 3.4 (1)を(x,y)=(b,−r)に適用してz−b∈Fであり、b−z∈Fである。(a−a′)+(b−z)=(a−s+z)+(b−z)=t−s=rである。
三つの場合のいずれでもs−z、a−a′、b−zはFに属し、(a−a′)+(b−z)=r∈Fである。§E20.1 系 3.2 (1)により各演算の結果は厳密な結果に等しく、e=fl(r)=rであるから、(1)と(2)が成り立つ。
(3)を示す。e=rであるから∣e∣≤min{∣a∣,∣b∣}である。tが正規範囲にあるならば§E20.1 系 3.2 (2)により∣e∣=∣t−fl(t)∣≤u∣t∣であり、∣t∣<2eminならば§E20.1 系 3.2 (4)によりs=t、e=0である。▨
4 補償和
定義 4.1.F=F(2,p,emin,emax)を基数2の浮動小数点数系、flをFの最近接丸めとし、n≥2を整数、x=(x1,…,xn)∈Fnとする。s1:=x1と置き、i=2,…,nについて(si−1,xi)の TwoSum を(si,ei)とする。(e2,…,en)∈Fn−1の逐次和をcと置き、s^:=fl(sn+c)をxの 補償和 (compensated summation) という。s^は、各 TwoSum とcの計算の各加算と最後の加算の厳密な結果の絶対値がNmax以下であるときに定まる。
定理 4.2.F=F(2,p,emin,emax)を基数2の浮動小数点数系、uをその単位丸め誤差、flをFの最近接丸めとし、n≥2を整数、x∈Fnとする。各i=2,…,nについて(a,b)=(si−1,xi)が定理 3.6の仮定を満たし、cの計算の各加算とsn+cの厳密な結果の絶対値がNmax以下であるとする。整数k≥0がku<1を満たすときγk:=ku/(1−ku)と置く。
- ∑i=1nxi=sn+∑i=2neiが成り立つ。
- (n−2)u<1ならば、補償和s^は
s^−i=1∑nxi≤ui=1∑nxi+(1+u)γn−2i=2∑n∣ei∣
を満たす。
- (n−1)u<1ならば∑i=2n∣ei∣≤γn−1∑i=1n∣xi∣であり、
s^−i=1∑nxi≤ui=1∑nxi+(1+u)γn−2γn−1i=1∑n∣xi∣
が成り立つ。
証明.定理 3.6 (2)によりsi−1+xi=si+ei(2≤i≤n)であり、これらを加えてs1=x1を用いると(1)を得る。
(2)を示す。S:=∑ixi、E:=∑i=2neiと置く。cはn−1個の成分の逐次和であるから、系 1.4 (1)により∣c−E∣≤γn−2∑i∣ei∣である。§E20.1 系 3.3により∣δ∣≤uを満たすδが存在してs^=(sn+c)(1+δ)であり、(1)によりsn=S−Eであるから
s^−S=(S+(c−E))(1+δ)−S=δS+(1+δ)(c−E)である。したがって∣s^−S∣≤u∣S∣+(1+u)γn−2∑i∣ei∣である。
(3)を示す。定理 3.6 (3)により∣ei∣≤u∣si−1+xi∣である。si=fl(si−1+xi)であるから、si−1は(x1,…,xi−1)の逐次和であり、系 1.4 (1)により∣si−1−∑j<ixj∣≤γi−2∑j<i∣xj∣である。したがって
∣si−1+xi∣≤(1+γi−2)j≤i∑∣xj∣≤(1+γn−2)j=1∑n∣xj∣であり、i=2,…,nについて加えて
i=2∑n∣ei∣≤(n−1)u(1+γn−2)j∑∣xj∣=1−(n−2)u(n−1)uj∑∣xj∣≤1−(n−1)u(n−1)uj∑∣xj∣=γn−1j∑∣xj∣を得る。これを(2)に代入する。▨
例 4.3.F=F(2,53,−1022,1023)、flを最近接偶数丸めとし、例 2.3のx=(1016,1,1,−1016)の補償和を計算する。si=fl(si−1+xi)は逐次和の部分和であり、s2=s3=1016、s4=0である。定理 3.6 (2)によりei=si−1+xi−siであるから、e2=e3=1、e4=0である。c=fl(fl(1+1)+0)=2であり、s^=fl(0+2)=2=σ(x)である。
例 4.4.F=F(2,53,−1022,1023)、flを最近接偶数丸め、u=2−53とし、x:=(1,2−53,2−106,2−106)とその成分を並べ替えたx′:=(1,2−106,2−106,2−53)を考える。厳密な和はどちらもS:=1+2−53+2−105であり、Sのただ一つの最近接点は1+2−52である。例 1.5と同じくfl(1+2−53)=1である。1+2−106のただ一つの最近接点は1であり、fl(1+2−106)=1である。2−53+2−106の最近接点は2−53=252s−53と2−53+2−105=(252+1)s−53であり、整数252が偶数であるからfl(2−53+2−106)=2−53である。xの補償和では、s2=s3=s4=1、(e2,e3,e4)=(2−53,2−106,2−106)、c=fl(fl(2−53+2−106)+2−106)=2−53であり、s^=fl(1+2−53)=1である。誤差S−s^=2−53+2−105はu∣S∣=2−53+2−106+2−158より大きく、定理 4.2 (2)の上界u∣S∣+(1+u)γ2(2−53+2−105)以下である。x′の補償和では、s2=s3=s4=1、(e2,e3,e4)=(2−106,2−106,2−53)、fl(2−106+2−106)=2−105、c=fl(2−105+2−53)=2−53+2−105であり、s^=fl(1+2−53+2−105)=1+2−52である。xとx′の逐次和はどちらも1である。
5 内積
命題 5.1.Fを浮動小数点数系、uをその単位丸め誤差、flをFの最近接丸めとし、n∈N≥1、x,y∈Fn、Tを{1,…,n}上の加算木とする。整数d≥0が(d+1)u<1を満たし、すべてのiについてdT(i)≤dであるとする。各iについて厳密な積xiyiが0であるか正規範囲にあるとし、p^:=(fl(x1y1),…,fl(xnyn))についてTの各加算の厳密な結果の絶対値がNmax以下であるとする。γd+1:=(d+1)u/(1−(d+1)u)と置く。このとき∣θi∣≤γd+1を満たす実数θ1,…,θnが存在してvT(p^)=∑ixiyi(1+θi)が成り立ち、
vT(p^)−i=1∑nxiyi≤γd+1i=1∑n∣xiyi∣である。特にnu<1ならば、積を逐次和で加えた計算値vRn(p^)について、この評価がγnで成り立つ。
証明.xiyi=0ならば§E20.1 系 3.2 (1)によりfl(xiyi)=xiyiであり、xiyiが正規範囲にあるならば§E20.1 系 3.2 (2)により、いずれの場合も∣εi∣≤uを満たすεiが存在してfl(xiyi)=xiyi(1+εi)である。定理 1.3 (1)によりvT(p^)=∑ixiyi(1+εi)∏j=1dT(i)(1+δi,j)であり、各iの因子の個数dT(i)+1は1以上d+1以下である。§E20.1 補題 4.1をk=dT(i)+1として適用し、k↦ku/(1−ku)の単調性を用いてθiを得る。Rnの深さは系 1.4 (1)によりn−1以下であるから、d=n−1として最後の主張を得る。▨
例 5.2.F=F(2,53,−1022,1023)、flを最近接偶数丸め、u=2−53とし、x:=(1+2−30,−1)、y:=(1−2−30,1)とする。成分はすべて正規範囲のFの元である。x1y1=1−2−60の最近接点は、1との距離2−60が1−2−53との距離2−53−2−60より小さいので1だけであり、fl(x1y1)=1である。fl(x2y2)=−1であり、逐次和はfl(1−1)=0であって加算は厳密であるから、厳密な内積−2−60との誤差はすべて積x1y1の丸めから生じ、相対誤差は1である。∑i∣xiyi∣=2−2−60であり、命題 5.1の上界はγ2(2−2−60)≈4.44×10−16である。c=yとして補題 2.1 (1)を適用すると、yを固定してxの関数と見た内積の成分ごとの相対条件数は(2−2−60)/2−60=261−1である。
6 Horner 法
定義 6.1.Fを浮動小数点数系、flをFの最近接丸めとし、n∈N≥1、a0,…,an∈F、x∈Fとする。qn:=anと置き、k=n−1,…,0の順にqk:=fl(fl(qk+1x)+ak)と置く。q0を、多項式∑i=0naiziのxにおける値の Horner 法 (Horner's method) による計算値という。q0は各乗算と各加算の厳密な結果の絶対値がNmax以下であるときに定まる。
定理 6.2.Fを浮動小数点数系、uをその単位丸め誤差、flをFの最近接丸めとし、n∈N≥1、a0,…,an∈F、x∈F、p(x):=∑i=0naixiとする。2nu<1とし、γ2n:=2nu/(1−2nu)と置く。k=n−1,…,0の各段で、厳密な積qk+1xが0であるか正規範囲にあり、厳密な和fl(qk+1x)+akの絶対値がNmax以下であるとする。
- ∣θi∣≤γ2nを満たす実数θ0,…,θnが存在して、Horner 法による計算値はq0=∑i=0nai(1+θi)xiを満たす。
- ∣q0−p(x)∣≤γ2n∑i=0n∣ai∣∣x∣iが成り立つ。Px:Rn+1→RをPx(a0,…,an):=∑i=0naixiで定めると、p(x)=0ならばκcomp(Px,(a0,…,an))=∑i∣ai∣∣x∣i/∣p(x)∣であり、
∣p(x)∣∣q0−p(x)∣≤γ2nκcomp(Px,(a0,…,an))
が成り立つ。
証明.0≤k≤nとk≤i≤nに対して、mk,n:=2(n−k)、mk,i:=2(i−k)+1(i<n)と置く。qn=anは、mn,n=0であるからqn=∑i=nnaixi−n∏j=1mn,i(1+δn,i,j)の空積の表示をもつ。k<nとし、∣δk+1,i,j∣≤uを満たす実数によりqk+1=∑i=k+1naixi−k−1∏j=1mk+1,i(1+δk+1,i,j)と書けるとする。qk+1xが0ならば§E20.1 系 3.2 (1)により、正規範囲にあるならば§E20.1 系 3.2 (2)により、∣μ∣≤uを満たすμが存在してfl(qk+1x)=qk+1x(1+μ)である。§E20.1 系 3.3により∣α∣≤uを満たすαが存在して
qk=(qk+1x(1+μ)+ak)(1+α)=qk+1x(1+μ)(1+α)+ak(1+α)である。qk+1の表示を代入すると、i≥k+1の項の因子の個数はmk+1,i+2=mk,i、akの因子の個数は1=mk,kであり、qk=∑i=knaixi−k∏j=1mk,i(1+δk,i,j)、∣δk,i,j∣≤uと書ける。kの降順の帰納法によりk=0でこの表示が成り立つ。1≤m0,i≤2nであるから、§E20.1 補題 4.1を因子の個数m0,iに適用して得るθiは∣θi∣≤m0,iu/(1−m0,iu)≤γ2nを満たし、(1)が成り立つ。(1)から∣q0−p(x)∣=∣∑iaiθixi∣≤γ2n∑i∣ai∣∣x∣iである。Pxはc=(1,x,…,xn)に対する補題 2.1の線形写像であるから、補題 2.1 (1)により条件数の式が成り立ち、両辺を∣p(x)∣で割って最後の不等式を得る。▨
系 6.3.n∈N≥1とし、Rn+2にノルム∥⋅∥∞、Rに絶対値を入れ、Hn:Rn+2→RをHn(a0,…,an,x):=∑i=0naixiで定める。Snを、単位丸め誤差uΦが4nuΦ≤1を満たす浮動小数点数系Φの全体とし、各Φ∈Snに最近接丸めを一つ固定する。DΦを、(a0,…,an,x)∈Φn+2∖{0}であって定理 6.2の各段の仮定を満たすものの全体とし、F^Φ(a0,…,an,x)を Horner 法による計算値q0とする。このとき、nを添字とする族は定数(4n)nで後退安定である。
証明.Φ∈Snと(a,x)=(a0,…,an,x)∈DΦをとる。2nuΦ≤1/2<1であるから定理 6.2 (1)が適用され、∣θi∣≤γ2n=2nuΦ/(1−2nuΦ)≤4nuΦである。Δ:=(a0θ0,…,anθn,0)はHn((a,x)+Δ)=∑iai(1+θi)xi=q0と∥Δ∥∞≤4nuΦmaxi∣ai∣≤4nuΦ∥(a,x)∥∞を満たす。したがってq0の相対後退誤差は4nuΦ以下であり、§E20.2 定義 5.2 (2)の条件がcn=4nで成り立つ。▨
例 6.4.F=F(2,53,−1022,1023)、flを最近接偶数丸め、u=2−53とし、(z−1)3=z3−3z2+3z−1の係数(a0,a1,a2,a3)=(−1,3,−3,1)とx:=1+2−20をとる。p(x)=2−60である。
Horner 法の各段はq3=1、q2=fl(x−3)、q1=fl(fl(q2x)+3)である。x−3=−2+2−20、(−2+2−20)x=−2−2−20+2−40、−2−2−20+2−40+3=1−2−20+2−40は§E20.1 補題 1.2 (1)によりFに属するので、q2=−2+2−20、q1=1−2−20+2−40である。q1x=1+2−60の最近接点は1だけであり、q0=fl(1−1)=0である。すべての積は正規範囲にあり、誤差は積q1xの丸め一回から生じる。∑i∣ai∣∣x∣i=(1+x)3=(2+2−20)3であるから、補題 2.1 (2)によりq0の成分ごとの相対後退誤差は2−60/(2+2−20)3≈1.08×10−19≈9.8×10−4uであり、定理 6.2の上界γ6≈6u以下である。一方、相対前進誤差は1である。定理 6.2 (2)の条件数は(2+2−20)3/2−60≈9.22×1018であり、相対前進誤差は条件数と成分ごとの相対後退誤差の積に等しい。xが根1に近づくとp(x)→0、∑i∣ai∣∣x∣i→8であるから、条件数は+∞に発散する。根1とxからd:=fl(x−1)、fl(fl(d⋅d)⋅d)を計算すると、d=2−20と二回の乗算fl(d⋅d)=2−40、fl(2−40d)=2−60はすべて厳密であり、計算値は2−60=p(x)に等しい。
7 二次方程式の小さい根
命題 7.1.B>2を実数とし、U:={(a,b,c)∈R3∣a>0, b>0, c>0, b2−4ac>0}、R:U→RをR(a,b,c):=(−b+b2−4ac)/(2a)で定める。r:=R(1,B,1)と置く。
- (a,b,c)∈Uならば、az2+bz+cの二つの根はR(a,b,c)とc/(aR(a,b,c))であり、どちらも負であってR(a,b,c)>c/(aR(a,b,c))である。z^<0がaz2+bz+cの根でありz^2<c/aを満たすならば、z^=R(a,b,c)である。
- κcomp(R,(1,B,1))=2B/B2−4である。
- −1<z^<0ならばηcomp(R,(1,B,1),z^)=∣z^2+Bz^+1∣/(z^2+B∣z^∣+1)であり、この下限は最小値である。
証明.(1)を示す。D:=b2−4ac>0であるから、az2+bz+cの根はz±:=(−b±D)/(2a)の二つであり、z+=R(a,b,c)>z−である。z+z−=c/a>0かつz++z−=−b/a<0であるから、二つの根はどちらも負であり、z−=c/(az+)である。z^<0が根でありz^=z+ならばz^=z−<z+<0であるから∣z^∣>∣z+∣であり、z^2>∣z^∣∣z+∣=c/aである。
(2)を示す。x0:=(1,B,1)と置き、§E20.2 定義 5.4のUx0とG:=Gx0を用いる。Rは多項式、(0,∞)上の平方根、a>0での除算の合成であるからU上でC1級であり、G(δ)=R(x0+x0∘δ)は0を含む開集合Ux0上で全微分可能である。(1)により、任意のδ∈Ux0に対して(1+δ1)G(δ)2+B(1+δ2)G(δ)+(1+δ3)=0が成り立つ。この恒等式をδ=0で微分し、G(0)=rを用いると、任意のh∈R3に対して(2r+B)DG(0)h+r2h1+Brh2+h3=0を得る。2r+B=B2−4=0であるからDG(0)h=−(r2h1+Brh2+h3)/B2−4である。∥h∥∞=1ならば∣r2h1+Brh2+h3∣≤r2+B∣r∣+1であり、(1)によりr<0であるから、h=(1,−1,1)で等号が成り立つ。したがって∥DG(0)∥=(r2+B∣r∣+1)/B2−4であり、§E20.2 命題 3.2 (2)によりκabs(G,0)=(r2+B∣r∣+1)/B2−4である。r2+1=−Br=B∣r∣であるからκabs(G,0)=2B∣r∣/B2−4であり、κcomp(R,x0)=κabs(G,0)/∣r∣=2B/B2−4である。
(3)を示す。−1<z^<0とし、L(a,b,c):=az^2+bz^+c、M:=z^2+B∣z^∣+1、η∗:=∣L(x0)∣/Mと置く。L(x0)=z^2+Bz^+1である。δ∈Ux0がR(x0+x0∘δ)=z^を満たすならばL(x0+x0∘δ)=0であるから、補題 2.1 (2)により∥δ∥∞≥η∗である。(z^2,Bz^,1)の符号は(1,−1,1)であるから、補題 2.1 (2)のy^=0に対するδ∗は、τ:=L(x0)/Mと置いてδ∗=(−τ,τ,−τ)であり、∥δ∗∥∞=η∗とL(x0+x0∘δ∗)=0を満たす。z^=0であるから∣L(x0)∣=∣z^2+1−B∣z^∣∣<Mであり、∣τ∣<1である。(a′,b′,c′):=(1−τ,B(1+τ),1−τ)の成分は正であり、z^はa′z2+b′z+c′の実根であるからb′2−4a′c′≥0である。b′2−4a′c′=0ならばz^は重根であってz^2=c′/a′=1となり、∣z^∣<1に反する。よって(a′,b′,c′)∈Uであり、z^2<1=c′/a′であるから、(1)によりR(a′,b′,c′)=z^である。したがってδ∗∈Ux0は下限を達成する。▨
例 7.2.F=F(2,53,−1022,1023)、flを最近接偶数丸め、u=2−53とし、B:=108とする。z2+Bz+1の二つの根をr:=(−B+w)/2、ℓ:=(−B−w)/2(w:=B2−4)とすると、根と係数の関係によりrℓ=1である。平方根も正しく丸めてw^:=fl(fl(fl(B⋅B)−4))と置き、小さい根rの二つの計算値を
r^D:=fl(fl(−B+w^)/2),r^V:=fl(1/ℓ^),ℓ^:=fl(fl(−B−w^)/2)で定める。命題 7.1 (2)により、係数(1,B,1)から小さい根への写像Rの成分ごとの相対条件数は2B/B2−4=2/1−4⋅10−16であり、2と2+10−15の間にある。r^Dの分子は、写像ω↦ω−Bにwの近似値w^を入れた値である。§E20.2 例 3.4 (3)により、この写像のwにおける相対条件数はw/(B−w)=w(B+w)/4であり、4.99×1015より大きい。ℓ^の分子の絶対値は、写像ω↦ω+Bにw^を入れた値である。同じ例をa=−Bとして、この写像のwにおける相対条件数はw/(w+B)<1である。ω↦1/ωの相対条件数は§E20.2 系 3.3により1である。
例 7.4.注意 7.3のr^D=−2−27とr^Vはどちらも(−1,0)に属するので、命題 7.1 (3)により成分ごとの相対後退誤差は∣z^2+Bz^+1∣/(z^2+B∣z^∣+1)である。r^Dでは分子が1−108⋅2−27+2−54≈0.2549419403076172、分母が1+108⋅2−27+2−54≈1.7450580596923828であり、後退誤差は約0.14609367229452420である。r^Vの後退誤差は約3.9538719584935763×10−17≈0.356uである。条件数κ:=2B/B2−4と後退誤差ηの積は、r^Dで約0.2922、r^Vで約0.712uであり、注意 7.3の相対前進誤差0.2549と0.712uと同じ大きさである。