exp(pi * sqrt(n))が整数に近くなるnを探す
Posted feedbacks - Scheme
Gaucheだと任意精度の小数が扱えないので基本的に有理数で計算しています。 ただしまともにやるととんでもない時間がかかるので、以下のように計算をはしょっています。 (1) 問題は「整数に近いかどうか？」というものなので、exp(pi*sqrt(n))を級数展開した ときの各項では小数部分のみ計算している。これで、各項の最終的な精度はdoubleで十分になる。 (2) 分母、分子の桁数が大きくなるのを防ぐため、偶数項と奇数項で計算を分けている (偶数項ではnの平方根はいらなくなるから)。 とはいえ、これでも結構時間がかかっていて、MacPro 3GHzで1500秒程度でした。
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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 | (use srfi-1)
(define (calc-pi k)
(do ((n 1 (+ n 1))
(a 1 (* a (- (* 4 n n) (* 2 n))))
(b 2 (* b 16))
(c 1 (* c n))
(p 0 (+ p (/ a (* b c c (- (* 2 n) 1))))))
((<= k n) (* 6 p))))
(define (%sqrt* a n)
(let ((root (sqrt a)))
(if (= (* root root) a)
(inexact->exact root)
(do ((i 1 (+ i 1))
(xn a (/ (+ xn (/ a xn)) 2)))
((<= n i) xn)))))
(define sqrt* (cut %sqrt* <> 10))
(define pi* (calc-pi 50))
(define (fract/ n d)
(define (str->dbl str)
(string->number
(string-append "0." (substring str 0 (min 15 (string-length str))))))
(let ((dstr (number->string d))
(nstr (number->string (modulo n d))))
(let ((d* (str->dbl dstr))
(n* (str->dbl nstr)))
(/. n* d* (expt 10 (- (string-length dstr) (string-length nstr)))))))
(define (%fract-exp-pi* n epsilon)
(let ((q1 (sqrt* n)))
(let ((pin (numerator pi*))
(pid (denominator pi*))
(q1n (numerator q1))
(q1d (denominator q1)))
(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)))
(new-fract (+ fract af)))
(if (< 0.0 af epsilon)
new-fract
(even-term (+ i 1) new-pn new-pd new-fract))))
(define (even-term i pn pd fract)
(let* ((new-pn (* pn pin n))
(new-pd (* pd pid i))
(af (fract/ new-pn new-pd))
(new-fract (+ fract af)))
(if (< 0.0 af epsilon)
new-fract
(odd-term (+ i 1) new-pn new-pd new-fract))))
(odd-term 1 1 1 0))))
(define fract-exp-pi* (cut %fract-exp-pi* <> 0.00001))
(print (filter (lambda (n)
(let1 x (fract-exp-pi* n)
(<= (abs (- x (round x))) 0.0001)))
(iota 199 1)))
|
herumi
#3416()
Rating0/2=0.00
Pythonで34行のスクリプトを書いて得られた出力の例が下のようになります。
この問題は光成さんに教えて頂いた e^{π*sqrt{163}}≒26253741640768744 が元になっています。ご協力ありがとうございました。[ reply ]