exp(pi * sqrt(n))が整数に近くなるnを探す in diff DouKaku?

1以上200未満の整数nのうち、 exp(pi * sqrt(n))がほとんど整数であるようなnを求めるコードを書いてください。なお、expは底がeである指数関数 - Wikipedia、 piは円周率、sqrtは平方根です。また「ほとんど整数である」とは整数からプラスマイナス0.0001の範囲にあることとします。

Pythonで34行のスクリプトを書いて得られた出力の例が下のようになります。

37 199148647.999978
58 24591257752.000000
67 147197952743.999999
163 262537412640768744.000000

この問題は光成さんに教えて頂いた e^{π*sqrt{163}}≒26253741640768744 が元になっています。ご協力ありがとうございました。

あああ、すみません。組み終わった後に ".00..." という表記を自分好みの "0.00.." という表記に書き換えたときにエンバグしちゃったみたいで、37 が表示されなくなってました。プラス評価もらったのに申し訳ない。

--- a.orig	Sat Sep 15 15:12:33 2007
+++ a.py	Sat Sep 15 15:12:50 2007
@@ -16,7 +16,7 @@
     for n in xrange(1, 200):
         x = decimal_exp(pi * Decimal(n).sqrt()) # result
         i = x.quantize(Decimal("1.0"), rounding=ROUND_HALF_UP) # nearest integer
-        if (x - i).quantize(Decimal("0.00001"), rounding=ROUND_DOWN) == 0:
+        if (x - i).quantize(Decimal("0.0001"), rounding=ROUND_DOWN) == 0:
             print n, x
 
 if __name__ == '__main__':

すみません、さっきのコードだと(pi * sqrt(n))^k / k! がたまたま整数に近かったときに
誤動作する可能性がありました(今回の解の探索範囲では問題ないのですが)。

--- eps.scm.orig        2007-09-16 00:31:43.000000000 +0900
+++ eps.scm     2007-09-16 00:32:29.000000000 +0900
@@ -39,9 +39,11 @@
       (define (odd-term i pn pd fract)
         (let* ((new-pn (* pn pin))
                (new-pd (* pd pid i))
-               (af (fract/ (* new-pn q1n) (* new-pd q1d)))
+               (an (* new-pn q1n))
+               (ad (* new-pd q1d))
+               (af (fract/ an ad))
                (new-fract (+ fract af)))
-          (if (< 0.0 af epsilon)
+          (if (and (< 0.0 af epsilon) (< an ad))
               new-fract
               (even-term (+ i 1) new-pn new-pd new-fract))))
       (define (even-term i pn pd fract)
@@ -49,7 +51,7 @@
                (new-pd (* pd pid i))
                (af (fract/ new-pn new-pd))
                (new-fract (+ fract af)))
-          (if (< 0.0 af epsilon)
+          (if (and (< 0.0 af epsilon) (< new-pn new-pd))
               new-fract
               (odd-term (+ i 1) new-pn new-pd new-fract))))
       (odd-term 1 1 1 0))))
@@ -60,3 +62,4 @@
                  (let1 x (fract-exp-pi* n)
                    (<= (abs (- x (round x))) 0.0001)))
                (iota 199 1)))

このバグはPython2.5.2で修正されました。というわけでパッチです。

--- a.py    Thu Feb 28 02:32:31 2008
+++ a.py    Thu Feb 28 02:32:46 2008
@@ -1,10 +1,11 @@
 from decimal import *
-import itertools
+import itertools, sys
 
-def divmod(x, y): # workaround for python2.5 decimal mod bug.
-    div = x // y
-    mod = x - div * y
-    return div, mod
+if sys.version_info[:3] <= (2, 5, 1):
+    def divmod(x, y): # workaround for python2.5 decimal mod bug.
+        div = x // y
+        mod = x - div * y
+        return div, mod
 
 def decimal_exp(x):
     div, mod = divmod(x, Decimal("2.3025850929940456840179914546844"))

DouKaku?

exp(pi * sqrt(n))が整数に近くなるnを探す

Posted feedbacks - diff

Index

Feed

Other

Link

Pathtraq