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



In [1]:

    
import pint as pn
from pint import roundfloat as rf
from pint import roundmode as rdm
import fractions



In [2]:

    
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

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



In [3]:

    
a = 0.9735930826234084
# 0.9735930826234084 例外値



In [4]:

    
itv_a = pn.interval(a)



In [5]:

    
format(itv_a.inf, '.17g')









    Out[5]:





'0.97359308262340838'



In [6]:

    
format(itv_a.sup, '.17g')









    Out[6]:





'0.97359308262340838'

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



In [7]:

    
itv_a.inf = rf.pred(a)
itv_a.sup = rf.succ(a)



In [8]:

    
format(itv_a.inf, '.17g')









    Out[8]:





'0.97359308262340827'



In [9]:

    
format(itv_a.sup,'.17g')









    Out[9]:





'0.97359308262340849'

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

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



In [10]:

    
aH, aL = rf.split(a)

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



In [11]:

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









    



False

以下はelseの続き

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



In [12]:

    
aS = rf.succ(a)
aP = rf.pred(a)



In [13]:

    
aS_c, aS_p = aS.as_integer_ratio()
aP_c, aP_p = aP.as_integer_ratio()



In [14]:

    
a_c, a_p = a.as_integer_ratio()

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



In [15]:

    
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は以下の性質を満たす

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

0.9735930826234084 例外値



In [16]:

    
bP_c / bP_p == a == bS_c / bS_p









    Out[16]:





False

これは常に成り立つ



In [17]:

    
bP_c / bP_p <= a <= bS_c / bS_p









    Out[17]:





True



In [18]:

    
bS_c / bS_p









    Out[18]:





0.9735930826234085



In [19]:

    
a









    Out[19]:





0.9735930826234084



In [20]:

    
bP_c / bP_p









    Out[20]:





0.9735930826234083

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

a <= bS



In [21]:

    
a_c * bS_p < bS_c * a_p









    Out[21]:





True

bP <= a



In [22]:

    
bP_c * a_p < a_c * bP_p









    Out[22]:





True

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

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



In [23]:

    
rf.pred(a) < bP_c / bP_p









    Out[23]:





False



In [24]:

    
bS_c / bS_p < rf.succ(a)









    Out[24]:





False

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



In [25]:

    
answer = pn.interval(0.)



In [26]:

    
answer.inf = fractions.Fraction(bP_c, bP_p)
answer.sup = fractions.Fraction(bS_c, bS_p)



In [27]:

    
print(answer)









    



[17538693776453097/18014398509481984,17538693776453099/18014398509481984]

今回は

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

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

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

つまり、分数は生成はするが、なるべく浮動小数点数に戻すようにする。