Jacobi記号を計算するRubyプログラム

Jacobi(ヤコビ)記号を計算するRubyプログラムの紹介。

アルゴリズム

Legendre記号やJacobi記号の計算は、相互法則を利用することにより、第一補充法則、第二補充法則、および(a/n)においてa=1の場合(あるいは、平方数の場合)に帰着します。実際、以下の手順(アルゴリズム)によって計算できます。

  1. (a/n)のaが負ならば、-1を分離して第一補充法則で計算する。
  2. (a/n)においてa = 1 (あるいは、aが平方数)ならば、終了。そうでなければ、以下を実行する。
  3. (a/n)においてa > nならば、aをnで割った最小正の余りに置き換え、ステップ2に戻る。そうでなければ、次へ。
  4. (a/n)のaが偶数ならば、2の因子を分離して第二補充法則で計算し、ステップ2に戻る。そうでなければ、次へ。
  5. (a/n)のaが奇数ならば、平方剰余の相互法則で(a/n)の上下を交換し、ステップ2に戻る。

a, nはどんどん小さくなり、最後にはa = 1になって(あるいは、aが平方数になって)、計算が終了します。

なお、aが負の場合をステップ3で対応すれば、ステップ1は省略可能です。

相互法則を利用してJacobi記号を計算するRubyプログラム

#!ruby -Ks

# Jacobi記号の計算
def jacobi_symbol(a, n)

  # Jacobi記号の値
  val = 1

  # 未定義の引数への対応
  if n < 3 or n % 2 == 0
    return  #=> nil
  end

  # a, nが互いに素でないときは0を返す
  if gcd(a, n) != 1
    return 0
  end

  # 第一補充法則を適用:
  if a < 0
    val *= (-1) ** ( (n-1)/2 )
    # val *= (n % 4 == 3) ? -1 : 1
    a = -a
  end

  while not_square?(a) do
    if a > n
      a %= n
    elsif a % 2 == 0
      a /= 4 while a % 4 == 0  # (4/n)=1
      if a % 2 == 0
        # 第二補充法則を適用
        val *= (-1) ** ( (n*n-1)/8 )
        # val *= (n % 8 == 3 or n % 8 == 5) ? -1 : 1
        a /= 2
      end
    else
      # 相互法則を適用
      val *= (-1) ** ( (a-1)*(n-1)/4 )
      # val *= (a % 4 == 3 and n % 4 == 3) ? -1 : 1
      n, a = a, n; 
    end
  end

  # (a/n)におけるaが平方数の場合(whileループ脱出)
  val
end

# 平方数でないかどうかの判定
def not_square?(a)
  (a > 0 && Math.sqrt(a).round ** 2 == a) ? false : true
end

# 最大公約数
def gcd(a, b)
  a = -a if a < 0
  b = -b if b < 0
  while b > 0
    a, b = b, a % b
  end
  a
end

# ======================================

if $0 == __FILE__ 

puts jacobi_symbol(439, 1201)  #=> -1

end

途中経過を表示

相互法則をどこで適用したのかといった、計算の途中経過を表示するように修正を加えたものが以下のプログラムです。

#!ruby -Ks

# Jacobi記号の計算
def jacobi_symbol(a, n)

  # Jacobi記号の値
  val = 1

  # 計算前の状態の表示
  print_progress("計算前:", val, a, n)

  # 未定義の引数への対応
  if n < 3 or n % 2 == 0
    print_progress("未定義です。", "nil")
    return  #=> nil
  end

  # a, nが互いに素でないときは0を返す
  if gcd(a, n) != 1
    print_progress("(a/n)においてgcd(a, n)>1の場合: ", 0)
    return 0
  end

  # 第一補充法則を適用:
  if a < 0
    val *= (-1) ** ( (n-1)/2 )
    # val *= (n % 4 == 3) ? -1 : 1
    a = -a
    print_progress("第一補充法則を適用:", val, a, n)
  end

  while not_square?(a) do
    if a > n
      a %= n
      print_progress("(a/n)のaをnで割った余りに置き換え:", val, a, n)
    elsif a % 2 == 0
      a /= 4 while a % 4 == 0  # (4/n)=1
      if a % 2 == 0
        # 第二補充法則を適用
        val *= (-1) ** ( (n*n-1)/8 )
        # val *= (n % 8 == 3 or n % 8 == 5) ? -1 : 1
        a /= 2
        print_progress("第二補充法則を適用:", val, a, n)
      else
        print_progress("(4/n)=1を適用:", val, a, n)
      end
    else
      # 相互法則を適用
      val *= (-1) ** ( (a-1)*(n-1)/4 )
      # val *= (a % 4 == 3 and n % 4 == 3) ? -1 : 1
      n, a = a, n;
      print_progress("相互法則を適用:", val, a, n)
    end
  end

  # (a/n)におけるaが平方数の場合(whileループ脱出)
  puts("(a/n)におけるaが平方数の場合:", val)
  val
end

# 平方数でないかどうかの判定
def not_square?(a)
  (a > 0 && Math.sqrt(a).round ** 2 == a) ? false : true
end

# 最大公約数
def gcd(a, b)
  a = -a if a < 0
  b = -b if b < 0
  while b > 0
    a, b = b, a % b
  end
  a
end

# 途中経過の表示
def print_progress(msg, val, a=nil, n=nil)
  puts msg
  if a && n
    puts (val < 0) ? "-(#{a}/#{n})" : "(#{a}/#{n})";
  else
    puts val
  end
end

# ======================================

if $0 == __FILE__ 

jacobi_symbol(439, 1201)

end

実行結果は、次の通りです。

計算前:
(439/1201)
相互法則を適用:
(1201/439)
(a/n)のaをnで割った余りに置き換え:
(323/439)
相互法則を適用:
-(439/323)
(a/n)のaをnで割った余りに置き換え:
-(116/323)
(4/n)=1を適用:
-(29/323)
相互法則を適用:
-(323/29)
(a/n)のaをnで割った余りに置き換え:
-(4/29)
(a/n)におけるaが平方数の場合:
-1

次回につづく。

【theme : 算数・数学の学習
【genre : 学校・教育

プロフィール

よしいず

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

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

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

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

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

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