0

最近、数年前に書いたコードでこれを見つけました。これは、適切な分母を決定し、元の実数と有理数の差が十分に小さいかどうかを確認することにより、(許容範囲内で) 実数を合理化するために使用されました。

明確にするために編集:実際には、すべての実際の値を変換したくありません。たとえば、分母の最大値として 14 を選択すると、7/15 に等しい実際の値はそのままになります。ここで書いたアルゴリズムの外部変数であるため、それは明確ではありません。

分母を取得するアルゴリズムは次のとおりです (疑似コード)。

denominator(x)
   frac = fractional part of x
   recip = 1/frac
   if (frac < tol)
      return 1
   else
      return recip * denominator(recip)
   end
end

連分数に基づいているようですが、再確認すると間違っていることがわかりました。(最終的には無限大を吐き出すだけなので、私にとってはうまくいきました。これは私が外部で処理しましたが、多くの場合、非常に遅くなります。) tol の値は、終了の場合または終了する数値の場合を除いて、実際には何もしません。近い。実数 - 有理数変換の許容範囲とは関係ないと思います。

高速であるだけでなく、理論的には失敗しないと確信している反復バージョンに置き換えました(最初はd = 1で、小数部は正を返すため、recipは常に> = 1です):

denom_iter(x d)
    return d if d > maxd 
    frac = fractional part of x
    recip = 1/frac
    if (frac = 0)
        return d
    else
        return denom_iter(recip d*recip)
    end
end

特定の許容範囲で可能なすべての値を確実に変換する maxd を選択する方法があるかどうか知りたいです。私は 1/tol を想定していますが、何かを見逃したくありません。また、このアプローチで分母のサイズを実際に制限する方法があるかどうかも疑問に思っています。これにより、一部の分母を maxd よりも大きくすることができます。

4

2 に答える 2

0

これに似た問題は、約 20 年頃に始まる近似セクションで解決されます。Bill Gosper の 連分数算術ドキュメントの 28 ページ。(参照:追記ファイル。1984行からのテキスト バージョンも参照してください。) 一般的な考え方は、2 つの分数が異なるまで、下限と上限の範囲制限数の連分数近似を計算し、値を選択することです。これら 2 つの近似値の範囲内です。これは、Gosper の用語を使用して、最も単純な分数を与えることが保証されています。

以下の python コード (プログラム「simpleden」) は、同様のプロセスを実装しています。(おそらく Gosper の提案した実装ほど良くはありませんが、メソッドが生成する結果の種類を確認するには十分です。) 実行される作業量は、Euclid のアルゴリズムの場合と同様です nビットであるため、プログラムはかなり高速です。テスト ケースの例 (プログラムの出力) は、コード自体の後に示されていますsimpleratio(vlo, vhi)vhi が vlo より小さい場合、ここに示す関数は -1 を返すことに注意してください。

#!/usr/bin/env python    
def simpleratio(vlo, vhi):
    rlo, rhi, eps = vlo, vhi, 0.0000001
    if vhi < vlo: return -1
    num  = denp = 1
    nump = den  = 0
    while 1:
        klo, khi = int(rlo), int(rhi)
        if klo != khi or rlo-klo < eps or rhi-khi < eps:
            tlo = denp + klo * den
            thi = denp + khi * den
            if tlo < thi:
                return tlo + (rlo-klo > eps)*den
            elif thi < tlo:
                return thi + (rhi-khi > eps)*den
            else:
                return tlo            
        nump, num = num, nump + klo * num
        denp, den = den, denp + klo * den
        rlo, rhi = 1/(rlo-klo), 1/(rhi-khi)

def test(vlo, vhi):
    den = simpleratio(vlo, vhi);
    fden = float(den)
    ilo, ihi = int(vlo*den), int(vhi*den)
    rlo, rhi = ilo/fden, ihi/fden;
    izok = 'ok' if rlo <= vlo <= rhi <= vhi else 'wrong'
    print '{:4d}/{:4d} = {:0.8f}  vlo:{:0.8f}   {:4d}/{:4d} = {:0.8f}  vhi:{:0.8f}  {}'.format(ilo,den,rlo,vlo, ihi,den,rhi,vhi, izok)

test (0.685, 0.695)
test (0.685, 0.7)
test (0.685, 0.71)
test (0.685, 0.75)
test (0.685, 0.76)
test (0.75,  0.76)
test (2.173, 2.177)
test (2.373, 2.377)
test (3.484, 3.487)
test (4.0,   4.87)
test (4.0,   8.0)
test (5.5,   5.6)
test (5.5,   6.5)
test (7.5,   7.3)
test (7.5,   7.5)
test (8.534537, 8.534538)
test (9.343221, 9.343222)

プログラムからの出力:

> ./simpleden 
   8/  13 = 0.61538462  vlo:0.68500000      9/  13 = 0.69230769  vhi:0.69500000  ok
   6/  10 = 0.60000000  vlo:0.68500000      7/  10 = 0.70000000  vhi:0.70000000  ok
   6/  10 = 0.60000000  vlo:0.68500000      7/  10 = 0.70000000  vhi:0.71000000  ok
   2/   4 = 0.50000000  vlo:0.68500000      3/   4 = 0.75000000  vhi:0.75000000  ok
   2/   4 = 0.50000000  vlo:0.68500000      3/   4 = 0.75000000  vhi:0.76000000  ok
   3/   4 = 0.75000000  vlo:0.75000000      3/   4 = 0.75000000  vhi:0.76000000  ok
  36/  17 = 2.11764706  vlo:2.17300000     37/  17 = 2.17647059  vhi:2.17700000  ok
  18/   8 = 2.25000000  vlo:2.37300000     19/   8 = 2.37500000  vhi:2.37700000  ok
 114/  33 = 3.45454545  vlo:3.48400000    115/  33 = 3.48484848  vhi:3.48700000  ok
   4/   1 = 4.00000000  vlo:4.00000000      4/   1 = 4.00000000  vhi:4.87000000  ok
   4/   1 = 4.00000000  vlo:4.00000000      8/   1 = 8.00000000  vhi:8.00000000  ok
  11/   2 = 5.50000000  vlo:5.50000000     11/   2 = 5.50000000  vhi:5.60000000  ok
   5/   1 = 5.00000000  vlo:5.50000000      6/   1 = 6.00000000  vhi:6.50000000  ok
  -7/  -1 = 7.00000000  vlo:7.50000000     -7/  -1 = 7.00000000  vhi:7.30000000  wrong
  15/   2 = 7.50000000  vlo:7.50000000     15/   2 = 7.50000000  vhi:7.50000000  ok
8030/ 941 = 8.53347503  vlo:8.53453700   8031/ 941 = 8.53453773  vhi:8.53453800  ok
24880/2663 = 9.34284641  vlo:9.34322100   24881/2663 = 9.34322193  vhi:9.34322200  ok

範囲内の最も単純な分数ではなく、分母のサイズの上限を考慮して最適な近似値を求める場合は、def test(vlo, vhi)前方のすべてのコードを置き換える次のようなコードを検討してください。

def smallden(target, maxden):
    global pas
    pas = 0
    tol = 1/float(maxden)**2
    while 1:
        den = simpleratio(target-tol, target+tol);
        if den <= maxden: return den
        tol *= 2
        pas += 1

# Test driver for smallden(target, maxden) routine
import random
totalpass, trials, passes = 0, 20, [0 for i in range(20)]
print 'Maxden  Num Den     Num/Den       Target        Error    Passes'
for i in range(trials):
    target = random.random()
    maxden = 10 + round(10000*random.random())
    den = smallden(target, maxden)
    num = int(round(target*den))
    got = float(num)/den
    print '{:4d}  {:4d}/{:4d} = {:10.8f}  = {:10.8f} + {:12.9f}  {:2}'.format(
        int(maxden), num, den, got, target, got - target, pas)
    totalpass += pas
    passes[pas-1] += 1
print 'Average pass count: {:0.3}\nPass histo: {}'.format(
    float(totalpass)/trials, passes)

プロダクション コードでは、pas(etc.) へのすべての参照を削除します。つまり、パス カウント コードを削除します。

ルーチンsmalldenには、許容される分母の目標値と最大値が与えられます。maxden分母の可能な選択が与えられた場合、 のオーダーの許容誤差を1/maxden²達成できると仮定するのは合理的です。次の典型的な出力に示されているパス カウント (ここでtarget、 とmaxdenは乱数によって設定されています) は、半分以上の時間でそのような許容値にすぐに達したことを示していますが、他のケースでは 2 倍、4 倍、または 8 倍の許容値が使用されたため、への余分な呼び出しsimpleratio。10000 回のテスト実行からの出力の最後の 2 行が、20 回のテスト実行の完全な出力の後に表示されていることに注意してください。

Maxden  Num Den     Num/Den       Target        Error    Passes
1198    32/ 509 = 0.06286837  = 0.06286798 +  0.000000392   1
2136   115/ 427 = 0.26932084  = 0.26932103 + -0.000000185   1
4257   839/2670 = 0.31423221  = 0.31423223 + -0.000000025   1
2680   449/ 509 = 0.88212181  = 0.88212132 +  0.000000486   3
2935   440/1853 = 0.23745278  = 0.23745287 + -0.000000095   1
6128   347/1285 = 0.27003891  = 0.27003899 + -0.000000077   3
8041  1780/4243 = 0.41951449  = 0.41951447 +  0.000000020   2
7637  3926/7127 = 0.55086292  = 0.55086293 + -0.000000010   1
3422    27/ 469 = 0.05756930  = 0.05756918 +  0.000000113   2
1616   168/1507 = 0.11147976  = 0.11147982 + -0.000000061   1
 260    62/ 123 = 0.50406504  = 0.50406378 +  0.000001264   1
3775    52/3327 = 0.01562970  = 0.01562750 +  0.000002195   6
 233     6/  13 = 0.46153846  = 0.46172772 + -0.000189254   5
3650  3151/3514 = 0.89669892  = 0.89669890 +  0.000000020   1
9307  2943/7528 = 0.39094049  = 0.39094048 +  0.000000013   2
 962   206/ 225 = 0.91555556  = 0.91555496 +  0.000000594   1
2080   564/1975 = 0.28556962  = 0.28556943 +  0.000000190   1
6505  1971/2347 = 0.83979548  = 0.83979551 + -0.000000022   1
1944   472/ 833 = 0.56662665  = 0.56662696 + -0.000000305   2
3244   291/1447 = 0.20110574  = 0.20110579 + -0.000000051   1
Average pass count: 1.85
Pass histo: [12, 4, 2, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

10000 回のテスト実行からの出力の最後の 2 行:

Average pass count: 1.77
Pass histo: [56659, 25227, 10020, 4146, 2072, 931, 497, 233, 125, 39, 33, 17, 1, 0, 0, 0, 0, 0, 0, 0]
于 2013-08-14T23:11:45.860 に答える