2016-01-31

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

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

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と厳密に同じかどうか考慮していない
分数演算の有効桁数に由来する計算誤差を考慮していない(分母、分子ともに巨大な整数になりがちである)
その後の分数の演算方法を考慮していない(作ったは良いが、使い方を考えていない)

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

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

2016-01-05

pythonの区間演算ライブラリ進捗2016-01-05

実は大体できた(唐突)。
今回はvmathクラスの紹介

github.com

import pint as pn

from pint import vmath

a = 2.

itv_a = pn.interval(a)

print(itv_a)

[2.0,2.0]

vmath.sqrt(a)

1.4142135623730951

floatを代入するとfloatが返ってくる

vmath.sqrt(itv_a)

[1.414213562373095,1.4142135623730951]

intervalを代入するとintervalが返ってくる
じゃあ、少し変なことをしてみる

itvmath = pn.vmath()

itvmath.arg_type = pn.interval

このitvmathのmath関数を使う時は代入する前にinterval型に変換するように設定

itvmath.sqrt(a)

[1.414213562373095,1.4142135623730951]

返ってくるのはinterval型になる。
itvmath.sqrt(pn.interval(a))と同じ

itvmath.sqrt(itv_a)

[1.414213562373095,1.4142135623730951]

itvmath2 = pn.vmath()

itvmath2.return_type = pn.interval

itvmath2.sqrt(a)

[1.4142135623730951,1.4142135623730951]

return_typeはreturnする直前に変換を行う
return return_type(answer)となるそのため、変換するタイミングが違うitvmath.sqrt(2.)とitvmath2.sqrt(2.)は結果が少し違う
このようにmath関数の振る舞いを自由に変更することができるのがvmathクラスの特徴です。
今はinterval型とビルトインされている型にしか対応してませんが、そのうち増やす予定

2015-12-21

メモ:as_interger_ratio()では正確に元の分数には戻せない

as_interger_ratio()では正確に元の分数には戻せない

a = 0.1 を区間型に代入する際にうまく包含できないので、色々実験してみた。

何が問題なのか?

import pint as pn

x=pn.interval(0.1)

と、2進数小数で表現不可能な数字なので下限と上限は

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

'0.10000000000000001'

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

'0.10000000000000001'

となり、0.1はうまく包含できない

解決策になりそうだったもの

pythonのfloat型にはas_integer_ratio()という関数があり、

a = 0.1

x,y = a.as_integer_ratio()

print(x)

3602879701896397

print(y)

36028797018963968

x/y

0.1

となり、0.1の分母と分数を取得できる
これを用いれば、さきほどのinterval型への代入問題を解決できそうだが

sup = pn.roundfloat.rddiv(x, y, pn.roundmode.up)

inf = pn.roundfloat.rddiv(x, y, pn.roundmode.down)

format(sup,'.17g')

'0.10000000000000001'

format(inf,'.17g')

'0.10000000000000001'

となり、0.1を包含するように計算できない。
なぜならば a=0.1 とした時点で 0.1 から、0.1よりずれた小数に近似されるため、as_interger_ratio で取得できる分母と分子は 0.1よりずれた小数に合うように取得されるためである。

誤解回避

0.1...となる分数を取得することは可能だが
厳密に0.1になるような分数はas_integer_ratioでは取得できない
(このように厳密に計算するとき以外は大した問題にならないと思われる)

2015-12-21

メモ:twosumアルゴリズムの研究

twosumアルゴリズムの研究

from pint import roundfloat as rf

from pint import roundmode as rdm

a=0.1

b=0.2

誤差が発生するパターン

rf.rdadd(a, b, rdm.nearest)

0.30000000000000004

pythonの区間演算ライブラリ進捗2015-11-16

成果

github.com
ソースコードが長くなりそうだったので分離した。
本当は一つのファイルで行きたかったけど、管理や読み込みなどを考えたら分離せざるを得なくなった。

github.com

list型を継承したinterval型も入る便利な配列を作る。
に該当するmcmatrixクラス
標準のmathを継承し、interval型対応と精度保証された値を出力するmathクラスを作る。の該当する予定のvmathクラスを作成

反省

独自の配列は結局、list型を継承して作ることにしたので、パフォーマンスは期待できない。
今後改良していく予定はもちろんあるが、numpyのndarrayと連携しようとすると、死にそう。

今後の方針

とりあえず、mcmatrixクラスに内積(dot)や転置、zerosやonesやらの実装をしないとまともに使えないので、さっさと実装する。

2015-10-25

pythonの区間演算ライブラリ進捗2015-10-25

成果

PEP8に対応した。
Githubにコードを公開した。

github.com

反省

PEP8対応のため、手動で途中まで直したが、autopep8コマンドで十分だった。(勉強にはなったけど)

今後の方針

list型を継承したinterval型も入る便利な配列を作る。
標準のmathを継承し、interval型対応と精度保証された値を出力するmathクラスを作る。

2015-10-13

pythonの区間演算ライブラリ進捗2015-10-13

成果

numpyのndarray型に対して、numpyのndarrayの加減乗っぽく計算する、仮想丸め方向付き計算クラスroundarrayを実装

反省

ndarrayの除算っぽい計算が実装できなかった(数値計算にとって必要なのかもよくわからない……)
double型の行列に対してtwosum,twoproductで誤差を拾おうとしているが、たぶん無理そう(そのうち別記事でまとめる予定)。

今後の方針

最初に掲げていた「interval型で簡単にplotしたい」は怪しくなってきた。
今後はnumpy依存部分と、標準ライブラリのみの部分を分けて開発をしていく。

numpy依存部分

引き続き「interval型で簡単にplotしたい」を実現するためにmatplotlibやらndarrayと戦っていく。
(でも、効率の良いinterval型の配列は中心値形式を利用しなければできそうにないので、時間がかかりそう)
しかし、使えるものは何でも使っていく予定

標準ライブラリのみ部分

「標準ライブラリでどこまで精度保証できるか?」を追求していく。
しかし、アルゴリズム以外の部分での高速化(BLASを使うetc...)はあまり求めないでいく。
ここまで来たら、「pythonのみで実装する」も方針に入れていこうと思う。
少し、コード周りを整理したら、githubで公開する予定。

ミンミンの日記

o108minminの近況や感想など

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

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

pythonの区間演算ライブラリ進捗2016-01-05

メモ:as_interger_ratio()では正確に元の分数には戻せない

as_interger_ratio()では正確に元の分数には戻せない

何が問題なのか?

解決策になりそうだったもの

誤解回避

メモ:twosumアルゴリズムの研究

twosumアルゴリズムの研究

pythonの区間演算ライブラリ進捗2015-11-16

成果

反省

今後の方針

pythonの区間演算ライブラリ進捗2015-10-25

成果

反省

今後の方針

pythonの区間演算ライブラリ進捗2015-10-13

成果

反省

今後の方針

numpy依存部分

標準ライブラリのみ部分