« スイートピー | トップページ | 本屋散歩の帰りがけ »

2014年12月 2日 (火)

e の近似分数を求める(少々力まかせな方法)

http://yamada-kuebiko.cocolog-nifty.com/blog/2014/11/post-e62b.html
2014年11月28日 (金)「連分数から近似分数を作る」
ここの終わりの方でこんなことを書きました↓

★ところで
「πとeの話」YEO・エイドリアン[著]、久保儀明+蓮見亮[訳]、青土社
この本のp.89にこんな記述があります。

eに小数点以下4桁まで正確な数値を与える、ある程度の調和美を備えている簡潔な分数は878/323である。

あれ?
連分数から作った近似分数の中に「878/323」はないですね。
そうなのか。連分数に依らずに、「878/323」を生み出すようなアルゴリズムを考えなくっちゃいけませんね。
あらまあ。いいけどさ。

  193 / 71  = 2.71830985915493
  878 / 323 = 2.7182662538699690402476780185759
1264 / 465 = 2.71827956989247
こんな風に、878/323が登場してくれるとよかったのですが。

★では、878/323が出てくるような近似分数づくりのプログラムを考えます。
αを無理数とし、m、nを整数とします。
近似分数というのは
α≒n/m
のことですね。
そこで、
mα≒n
としてαに整数mをかけ、その答えが整数に「近かったら」n/mを近似分数とすればよいわけです。
どの程度の差なら「近い」と考えるかは、自分で決めればいい。

e=2.71828171828…
me=n

m=1 n=2.71828171828…
m=2 n=5.43656343656
m=3 n=8.15484515484
m=4 n=10.87312687312
・・・
m=7 n=19.02797202796
・・・
メンドクサ

m=7の場合、大甘に見てn=19とし
n/m=19/7=2.714285714285714
これで小数点以下2桁目まで一致する近似分数だ、といっても構いません。
判断は人任せ。
(ナント)連分数からの近似分数での第4次近似として現れる近似分数でもあります。

というわけで、プログラムの姿は見えてきました。
ただ注意しないといけないのは、上の例でのm=2と4の場合。
10/4は10と4に公約数があるので既約ではない。公約数2で割ればm=2に帰着します。
プログラムを書くときに、n/mは既約である、という条件を付けなければいけませんね。
そのために、2数の最大公約数を求める必要があります。
こんな方針で10進BASICで書いたのが下のプログラムです。
「gcd」という関数はとても便利な関数なので、独立させてあちこちで使い回したい。
で「プログラム本体の外に書くよ」という宣言をしておくわけです。
そうしないと、インタープリタが、「gcd」なんて言われても知らないよ、となってしまいます。
gcdは公約数があればその値を、ない場合「1」を返す関数で、ユークリッドの互除法という方法を使って再帰的に簡潔に書きました。なお、gcdに与える2数を大きさ順にしておくべきかなどということは考える必要はありません。どう入れても大丈夫です。{gcdについては、稿を改めて書きます。}

!----------------------------------
DECLARE EXTERNAL FUNCTION gcd

LET a = 2.7182818284590452353602874713527
!LET a = 3.1415926535897932384626433832795
LET MAX = 500
LET D =0.0001

LET count=0
FOR m = 1 TO MAX
   LET n = INT((a*m)+0.5)       !四捨五入
   IF (ABS(a-n/m) <D) AND (gcd(n,m)=1) THEN
      PRINT n;"/";
      PRINT m;"=";
      PRINT n/m
      LET count =count+1
   END IF
NEXT m

PRINT
PRINT count

END
!----------------------------------
EXTERNAL FUNCTION gcd(a, b)
    IF b=0 THEN
        LET gcd = a
    ELSE
        LET gcd = gcd(b, MOD(a,b))
    END IF
END FUNCTION
!----------------------------------

結果ですが
LET MAX = 500{mは500まで}
LET D =0.0001
という条件で走らせると

  193 /  71 = 2.71830985915493
  492 / 181 = 2.7182320441989
  666 / 245 = 2.71836734693878
  685 / 252 = 2.71825396825397
  791 / 291 = 2.71821305841924
  859 / 316 = 2.71835443037975
  878 / 323 = 2.71826625386997
1052 / 387 = 2.71834625322997
1071 / 394 = 2.71827411167513
1090 / 401 = 2.71820448877805
1139 / 419 = 2.71837708830549
1177 / 433 = 2.71824480369515
1245 / 458 = 2.71834061135371
1264 / 465 = 2.71827956989247
1283 / 472 = 2.71822033898305

15
こうなって、15の近似分数が見つかりました。
メデタイ!ですね。878/323 があります。
878 / 323 = 2.71826625386997
小数点以下4桁まで正しい。当たり前ですね、D =0.0001 で判断していますから。

LET MAX = 10000
LET D =0.000001
でやると65個近似分数が得られまして

2721 / 1001 = 2.71828171828172
3985 / 1466 = 2.71828103683492
4178 / 1537 = 2.71828236824984
5635 / 2073 = 2.71828268210323
6706 / 2467 = 2.71828131333604
・・・
26525 / 9758 = 2.71828243492519
26631 / 9797 = 2.71828110646116
26718 / 9829 = 2.71828263302472
26911 / 9900 = 2.71828282828283
27017 / 9939 = 2.71828151725526
こうでした。
道具はできた、後はもう遊ぶだけです。
どうぞ。

« スイートピー | トップページ | 本屋散歩の帰りがけ »

理科おじさん」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

« スイートピー | トップページ | 本屋散歩の帰りがけ »

2021年5月
            1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30 31          
サイト内検索
ココログ最強検索 by 暴想
無料ブログはココログ