数値解析日記

数値解析関係のネタ帳。C++,pythonなど

分数区間を用いた、浮動小数点数の区間包含についてのアイデア

分数区間を用いた、浮動小数点数の区間包含について

python 3.5.1で動作を確認
区間演算の実装について(2)python上解決する場合のアイデア
いろいろと検証が足りていないが公開はする。
jupyter notebook用のコードもアップロードした。

import pint as pn
from pint import roundfloat as rf
from pint import roundmode as rdm
import fractions
def frac_interval_proto(a):
    aH, aL = rf.split(a)
    if aL ==0:
        answer = pn.interval(a)
    else:
        aS = rf.succ(a)
        aP = rf.pred(a)
        aS_c, aS_p = aS.as_integer_ratio()
        aP_c, aP_p = aP.as_integer_ratio()
        a_c, a_p = a.as_integer_ratio()
        bS_c = a_c * aS_p + aS_c * a_p
        bS_p = a_p * aS_p * 2
        bP_c = a_c  * aP_p + aP_c * a_p
        bP_p = a_p * aP_p * 2
        answer = pn.interval(0.)
        answer.inf = fractions.Fraction(bP_c, bP_p)
        answer.sup = fractions.Fraction(bS_c, bS_p)
    return answer

今回は、代入演算子を用いて区間型を生成した際に、初期値を包含しない区間が発生するのを防ぐ分数区間を生成する。
すなわち

a = 0.1
itv_a = pn.interval(a)
format(itv_a.inf, '.17g')
'0.10000000000000001'
format(itv_a.sup, '.17g')
'0.10000000000000001'

のような区間が発生することを防ぐ。
なお、今回以外の解決方法としては

itv_a.inf = rf.pred(a)
itv_a.sup = rf.succ(a)
format(itv_a.inf, '.17g')
'0.099999999999999992'
format(itv_a.sup,'.17g')
'0.10000000000000002'

のようにsucc predを用いることで解決可能である。
しかし、今回のアルゴリズムのほうが、より区間幅を縮小できる。

まず、上位ビットと下位ビットを分離する

aH, aL = rf.split(a)

もし、下位ビットが存在する場合、else文以降の通りにする。
(下位ビットが存在しない浮動小数点数で点区間を生成した場合、元の初期値を包含する区間になる可能性が高い
元の初期値の下位ビットが全て1で、近似の際にくりあがりで下位ビットが消去された場合、おかしくなるかも)

if aL ==0:
    print(True)
    answer = pn.interval(a)
else:
    print(False)
False

以下はelseの続き

次にaをsucc predし、as_intger_ratio()で近似分数にする

aS = rf.succ(a)
aP = rf.pred(a)
aS_c, aS_p = aS.as_integer_ratio()
aP_c, aP_p = aP.as_integer_ratio()
a_c, a_p = a.as_integer_ratio()

次に、aとsuccしたaの中点と、aとpredしたaの中点を求める。

bS_c = a_c * aS_p + aS_c * a_p
bS_p = a_p * aS_p * 2
bP_c = a_c  * aP_p + aP_c * a_p
bP_p = a_p * aP_p * 2

こうして得た、分数bSとbPとaは以下の性質を満たす

浮動小数点数として評価した場合、同値である。

2016-02-01 例外がたくさん発生した
正しくは bP_c / bP_p <= a <= bS_c / bS_p

bS_c / bS_p == a == bP_c / bP_p
True

分数として評価するとbP <= a <= bSが成立する
(以下は、分母と通分するためbS_pやら、a_pがかけられている。今回は bP < a < bSとなっていることをわかりやすくするため、等号の場合は考慮していない)

a <= bS

a_c * bS_p < bS_c * a_p 
True

bP <= a

bP_c * a_p < a_c * bP_p
True

よって、分数で評価すると大小関係が存在するが、浮動小数点数で評価すると値が一致する分数区間が得られた

なお、succ predで得られたitv_aと比較すると

rf.pred(a) < bP_c / bP_p
True
bS_c / bS_p < rf.succ(a)
True

となり、succ predよりも狭い区間幅が得られた。

answer = pn.interval(0.)
answer.inf = fractions.Fraction(bP_c, bP_p)
answer.sup = fractions.Fraction(bS_c, bS_p)
print(answer)
[14411518807585587/144115188075855872,14411518807585589/144115188075855872]

今回は

  • 引数が正の場合しか考慮していない
  • float.as_intger_ratio()で得られる分数が、最近点丸めで近似されたaと厳密に同じかどうか考慮していない
  • 分数演算の有効桁数に由来する計算誤差を考慮していない(分母、分子ともに巨大な整数になりがちである)
  • その後の分数の演算方法を考慮していない(作ったは良いが、使い方を考えていない)

などの問題がある
最後の分数の演算に関しては以下のようにするアイデアがあるが、検証はしていない

  • 分数 * 分数: 分数演算をしたあと、丸め制御で適切な浮動小数点数にする
  • 分数 * 浮動小数点数(下位ビットがない): 分子 * 小数を計算後、丸め制御で適切な浮動小数点数にする
  • 分数 * 浮動小数点数(下位ビットがある): 小数を分数に変換したあと、分数 * 分数に帰着させる