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

Pell 方程式 x^2-my^2=±4 の最小解 (=最小の正の整数解) を求める Maple プログラムの例。

連分数展開による方法

m の平方根の連分数展開を利用する方法。

solvePellEqCFrac4 := proc(m)
  local x, y, r, m1, d, p0, q0, p1, q1, p2, q2, s, t, a;

  # m が平方数ならエラー (非自明解なし)
  if issqr(m) then
    error("%1 is a square number.", m);
  end if;

  if m = 5 then
    x := 1; y := 1;
  elif m = 13 then
    x := 3; y := 1;
  else
    if modp(m, 4) = 0 then
      m1 := m/4;
    else
      m1 := m;
    end if;

    d := floor(sqrt(m1));
    p0 := 1; q0 := 0;
    p1 := d; q1 := 1;
    s := 0; t := 1; a := d;

    while true do
      s := a * t - s;
      t := floor((m1 - s * s) / t);
      a := floor((s + d) / t);
      p2 := a * p1 + p0;
      q2 := a * q1 + q0;
      p0 := p1; q0 := q1;
      p1 := p2; q1 := q2;

      if t = 4 and modp(m, 8) = 5 then
        x := p0; y := q0;
        break;
      elif t = 1 then
        if modp(m, 4) = 0 then
          x := 2*p0; y := q0; break;
        else
          x := 2*p0; y := 2*q0;
        end if;
        break;
      end if;
    end do;
  end if;

  r := x*x - m*y*y;

  [m, x, y, r]
end proc:

実行例

printf("[m, x, y, r]\n");
for m from 1 to 100 do
  if not issqr(m) then
    printf("%a\n", solvePellEqCFrac4(m));
  end if
end do

実行結果:

[m, x, y, r]
[2, 2, 2, -4]
[3, 4, 2, 4]
[5, 1, 1, -4]
[6, 10, 4, 4]
[7, 16, 6, 4]
[8, 2, 1, -4]
[10, 6, 2, -4]
[11, 20, 6, 4]
[12, 4, 1, 4]
[13, 3, 1, -4]
[14, 30, 8, 4]
[15, 8, 2, 4]
[17, 8, 2, -4]
[18, 34, 8, 4]
[19, 340, 78, 4]
[20, 4, 1, -4]
[21, 5, 1, 4]
[22, 394, 84, 4]
[23, 48, 10, 4]
[24, 10, 2, 4]
[26, 10, 2, -4]
[27, 52, 10, 4]
[28, 16, 3, 4]
[29, 5, 1, -4]
[30, 22, 4, 4]
[31, 3040, 546, 4]
[32, 6, 1, 4]
[33, 46, 8, 4]
[34, 70, 12, 4]
[35, 12, 2, 4]
[37, 12, 2, -4]
[38, 74, 12, 4]
[39, 50, 8, 4]
[40, 6, 1, -4]
[41, 64, 10, -4]
[42, 26, 4, 4]
[43, 6964, 1062, 4]
[44, 20, 3, 4]
[45, 7, 1, 4]
[46, 48670, 7176, 4]
[47, 96, 14, 4]
[48, 14, 2, 4]
[50, 14, 2, -4]
[51, 100, 14, 4]
[52, 36, 5, -4]
[53, 7, 1, -4]
[54, 970, 132, 4]
[55, 178, 24, 4]
[56, 30, 4, 4]
[57, 302, 40, 4]
[58, 198, 26, -4]
[59, 1060, 138, 4]
[60, 8, 1, 4]
[61, 39, 5, -4]
[62, 126, 16, 4]
[63, 16, 2, 4]
[65, 16, 2, -4]
[66, 130, 16, 4]
[67, 97684, 11934, 4]
[68, 8, 1, -4]
[69, 25, 3, 4]
[70, 502, 60, 4]
[71, 6960, 826, 4]
[72, 34, 4, 4]
[73, 2136, 250, -4]
[74, 86, 10, -4]
[75, 52, 6, 4]
[76, 340, 39, 4]
[77, 9, 1, 4]
[78, 106, 12, 4]
[79, 160, 18, 4]
[80, 18, 2, 4]
[82, 18, 2, -4]
[83, 164, 18, 4]
[84, 110, 12, 4]
[85, 9, 1, -4]
[86, 20810, 2244, 4]
[87, 56, 6, 4]
[88, 394, 42, 4]
[89, 1000, 106, -4]
[90, 38, 4, 4]
[91, 3148, 330, 4]
[92, 48, 5, 4]
[93, 29, 3, 4]
[94, 4286590, 442128, 4]
[95, 78, 8, 4]
[96, 10, 1, 4]
[97, 11208, 1138, -4]
[98, 198, 20, 4]
[99, 20, 2, 4]

素朴な方法

y に値を入れ、m*y^2±4 が平方数になるものを探します。

solvePellEqNaive4 := proc(m)
  local x, y, r, xm2, xp2;

  # m が平方数ならエラー (非自明解なし)
  if issqr(m) then
    error("%1 is a square number.", m);
  end if;

  y := 0;
  while true do
    y := y + 1;
    xm2 := m * y * y - 4;
    if issqr(xm2) then
      x := isqrt(xm2); break;
    end if;
    xp2 := xm2 + 8;  # m * y * y + 4
    if issqr(xp2) then
      x := isqrt(xp2); break;
    end if;
  end;

  r := x*x - m*y*y;

  [m, x, y, r]
end proc:

実行例

printf("[m, x, y, r]\n");
for m from 1 to 100 do
  if not issqr(m) then
    printf("%a\n", solvePellEqNaive4(m));
  end if
end do

実行結果は、連分数展開による方法と同じになります。

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

プロフィール

よしいず

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

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

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

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

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

リンク
月別アーカイブ
カテゴリ
最新記事
検索フォーム
RSSリンクの表示