ペル方程式の問題です。
って言ってますけど、わからないので、情報を調べていてこの方程式がペル方程式という物だということを知ったんですけどね。
解く方法はわかったんですが、残念なことに、どうしてこうすれば答が出るのか理解できてません。
そういうの多いんですけど、ちょっと悔しいね。
もとのペル方程式は、
x^2 - Dy^2 = ±1の形をしています。
これの最小解を求めるには、二次体(Q√D)がなんたらかんたらで、結局は連分数√Dの連分数展開の結果を使ってごにょごにょであるとのこと。詳しくはペル方程式でググってください。
やりかただけ書くと、まずは√Dの連分数展開を求めます。平方根の連分数展開は、小数点以下が必ずループするので、こんな形になります。
{a0: a1, a2, ... an, ,a1, ...}で、この、n-1番目までの結果を分数にしたものの分子と分母がそれぞれxとyの最小解になっているんだそうな。
連分数展開は前の問題で作ったのがあるけど、ちょっと改変が必要。
あと、結果は場合によっては、x^2 - Dy^2 = -1 の解が先に出る。
ペル方程式の解は無数にあって、n番目のもの xn,ynはここでもとめた最小解x0,y0と
xn - √Dyn = (x0 - √Dy0)^n
という関係があって、しかも、-1 の次は +1 であるとのことなので、2番目を求めればよい。
さて、よくわかってませんが、これを実装しました。
式の中では
√3 + 5のような形が基本データになるので、これを、
----------
4
√X + Mとして、基本形としています。
----------
N
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ;; Problem 66 : 2011/11/29 | |
| ;; "Elapsed time: 67.348225 msecs" | |
| ;; treat all num as | |
| ;; (sqrt X)+M | |
| ;; ----------- | |
| ;; N | |
| (defn int-floor [x] | |
| (Math/round (Math/floor x))) | |
| (defn inverse-rationalize [X M N] | |
| (let [new-N (/ (- X (* M M)) N) | |
| new-M (- M)] | |
| [new-M new-N])) | |
| (defn make-proper [X M N] | |
| (let [nume (+ (Math/sqrt X) M) | |
| a (int-floor (/ nume N)) | |
| new-M (- M (* a N))] | |
| [a new-M])) | |
| (defn get-cf-list-loop | |
| "Get continued fraction parms of Sqrt X" | |
| ([X] (get-cf-list-loop X 0 1 nil)) | |
| ([X M N top] | |
| (let [[a prop-M] (make-proper X M N) | |
| [new-M new-N] (inverse-rationalize X prop-M N)] | |
| (lazy-seq | |
| (cons a | |
| (if (= [new-M new-N] top) | |
| nil | |
| (get-cf-list-loop X | |
| new-M | |
| new-N | |
| (if (empty? top) | |
| [new-M new-N] | |
| top)))))))) | |
| (defn calc-cf-partial-fraction [cf-list] | |
| (reduce #(+ (/ 1 %1) %2) | |
| (drop 1 (reverse cf-list)))) | |
| (defn check-norm [X Y D] | |
| (= 1 (- (* X X) (* Y Y D)))) | |
| (defn my-numerator [n] | |
| (if (integer? n) | |
| n | |
| (numerator (rationalize n)))) | |
| (defn my-denominator [n] | |
| (if (integer? n) | |
| 1 | |
| (denominator (rationalize n)))) | |
| (defn solv-pell-eq [n] | |
| (let [fraction (calc-cf-partial-fraction (get-cf-list-loop n)) | |
| num (my-numerator fraction) | |
| den (my-denominator fraction)] | |
| (if (check-norm num den n) | |
| [num den n] | |
| [(+ (* num num) (* den den n)) (* 2 den num) n]))) | |
| (defn square? [n] | |
| (let [rt (Math/round (Math/sqrt n))] | |
| (= n (* rt rt)))) | |
| (defn pe-66 [limit] | |
| (reduce #(if (> (first %1) (first %2)) %1 %2) | |
| (map solv-pell-eq | |
| (filter (complement square?) (range 10 (inc limit)))))) | |
- int-floor
引数nを超えない最大の整数を返します。
- inverse-rationalize
基本形の分数の逆数を求めます。ただの逆にするのではなく、分母を有理化します。有理化すると形は基本形になり、平方根の中身Xは変らないので、更新したMとNの組を返します。
- make-proper
基本形の分数を帯分数化します。括り出した整数と新しいMを返します。
- get-cf-list-loop
引数Xの連分数展開の結果を返します。ただし、ループの終端までで返ります。
- calc-cf-partial-fraction
連分数展開のリストをもらって、分数に戻します。
- check-norm
x^2-Dy^2が1になっているかどうかを返します。
名前のつけかたを間違えていますね、norm-is-1? だったな。
- my-numerator
- my-denominator
Clojureは分数表記を扱えて、分子と分母を求める関数がある(
numerator,denominator)のですが、引数に整数を渡すと例外を出してしまうので、整数でも動くようにしたもの。マルチメソッドにすべきだったかな。
- solv-pell-eq
問題を解くメイン関数。
Dの値を引数に取って、√Dを連分数展開。
ループの直前までの分数の値を計算。
分母と分子をペル方程式にあてはめて、1ならそれ、-1なら2番目を返す。返るのは、 x,y,D。
- square?
平方数かどうか
- pe-66
1からlimitまでの平方数以外の数について、ペル方程式が1になる解を求めて、xの値が最大になるものを返します。
0 コメント:
コメントを投稿