スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

Pell方程式の最小解を求めるプログラム

Pell方程式 x^2-my^2=±1 の最小解を求めるプログラムの例。

素朴な方法

y に値を入れ、m*y^2±1 が平方数になるものを探す。Rubyによるプログラム。

m = 5

y = 0
loop do
  y += 1
  xm2 = m * y * y - 1
  xm = Math.sqrt(xm2).floor
  if xm * xm == xm2 then 
    puts "m = #{m}, x = #{xm}, y = #{y}"
    puts "x^2 - m * y^2 = -1"
    break
  end
  xp2 = xm2 + 2  # = m * y * y + 1
  xp = Math.sqrt(xp2).floor
  if xp * xp == xp2 then  
    puts "m = #{m}, x = #{xp}, y = #{y}"
    puts "x^2 - m * y^2 = 1"
    break
  end
end

連分数展開による方法(1)

m の平方根の連分数展開を利用する方法。途中で無理数の計算が出てくる。Mapleによるプログラム。

solvePellEq := proc(m) 
  local alpha, a, p0, p1, p2, q0, q1, q2, r;

  alpha := sqrt(m); 
  a := floor(alpha); 
  p0 := 0; p1 := 1; 
  q0 := 1; q1 := 0; 
  p2 := a*p1+p0; q2 := a*q1+q0;
  r := p2*p2+(-1)*m*q2*q2;
  while r <> 1 and r <> -1 do 
    alpha := 1/(alpha-a);  # 無理数の計算
    a := floor(alpha); 
    p0 := p1; p1 := p2; 
    q0 := q1; q1 := q2; 
    p2 := a*p1+p0; 
    q2 := a*q1+q0;
    r := p2*p2+(-1)*m*q2*q2;
  end do; 
  p2, q2, r
end proc:

> solvePellEq(199);
                         16266196520, 1153080099, 1

連分数展開による方法(2)

mの平方根の連分数展開において、補助的な数列 (s_n), (t_n) を導入することで、無理数の計算が回避できる。Rubyによるプログラム。

m = 199

d = Math.sqrt(m).floor
p0 = 1; q0 = 0;
p1 = d; q1 = 1;
s = 0; t = 1; a = d
begin
  s = a * t - s
  t = (m - s * s) / t
  a = (s + d) / t
  p2 = a * p1 + p0
  q2 = a * q1 + q0
  p0 = p1; q0 = q1
  p1 = p2; q1 = q2
end while t > 1
r = p0 * p0 - m * q0 * q0
puts
puts "m = #{m}, x = #{p0}, y = #{q0}"
puts "x^2 - #{m} * y^2 = #{r}"

【theme : 数学
【genre : 学問・文化・芸術

プロフィール

よしいず

Author:よしいず
MATHEMATICS.PDFというウェブサイトを運営しています。

管理の都合上、トラックバックとコメントはオフにしてあります。ブログ経験者なら分かっていただけると思いますが、スパム(アダルトやその他の宣伝)ばかりなのが現実です。

リンクは自由です。当サイトの記事に対する間違いの指摘・意見・感想などを述べた記事からのリンクは歓迎です。ただし、ブログ記事アップ直後はミスが多く、頻繁に修正します。場合によっては削除する可能性もあります。その際、何も断りもなく修正・削除しますがご了承ください。内容を参考にする場合には投稿後一週間ほど様子を見てからにしてください(笑)。

記事の間違いを指摘するときは、その具体的箇所、理由(仕様に反するなど)・根拠(参考にした文献など)、代替案(同じ結果を得るための正しいやり方)も教えてください。そうしないと、(指摘される側および第三者はその時点では無知の状態なので、)どこが間違いなのか分かりませんし、本当に間違っているのかどうかが判断・検証できません。実際、間違いだと指摘されたことが結局は正しかったというケースもありますので。

このブログのタイトル一覧

リンク
月別アーカイブ
カテゴリ
最新記事
検索フォーム
RSSリンクの表示
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。