元祖ワシ的日記

一つ前のエントリ


法Nでの逆数を求めるその２

はつまるところ、

\[ax + by = 0 \pmod N\]

を考え$a, b$をどんどん小さくしてゆき、上記等式を満たすように$x, y$を変えていくアルゴリズムと言えるでしょう。ここで説明する手法は両辺を2で割りより早く$a, b$を小さくするという手法です。ところで一つ前のエントリの補題により$a, b$は互いに素です。よって$a, b$が共に偶数ということはありません。しかし例えば$a, y$が共に偶数なら

\[(a/2)x + b(y/2) = 0 \pmod N\]

と係数をより小さくすることができます。ここでちょっと工夫をします。$\pmod N$で考えるのなら上記の式は

\[ax + b(y+N) = 0 \pmod N\]

と$y$に$N$を足しても成立します。$N$が奇数ならば$y+N$は偶数になりますので、つまるところ$y$が偶数だろうが奇数だろうが、$a$が偶数ならば$a/2$を考える事が出来ます。$b$についても同様で$x$が偶数ならば$x/2$とすればよく、奇数ならば$(x+N)/2$を考えればよいのです。 この作業の結果$a, b$は共に奇数になるはずです。今借りに$a>b$とすれば $a' = a-b, b' = b$とさらに係数を小さくすることができます。$a-b$は偶数ですから必ず2で割れるのでさらに係数を小さくしていくことができます。この作業によりいつかはどちらかの係数は1になります。そのときの$x$が求める答えです。

このアルゴリズムは$N$が奇数のときにしか使えませんが、シフトと加算減算だけで実行できるのでとても高速です。

サンプルコード

#!/usr/bin/env ruby



def inv(t, n)
  a = n
  b = t
  x = 1
  y = 0



  while b != 1
    shift = 0
    while (b & (1 << shift)) == 0
      if((x & 1) == 1)
        x = x + n
      end
      x /= 2
      shift += 1
    end
    b = b >> shift
    


    shift = 0
    while (a & (1 << shift)) == 0
      if((y & 1) == 1)
        y = y + n
      end
      y /= 2
      shift += 1
    end
    a = a >> shift
    


    if b > a
      x = x + y
      b = b - a
    else
      y = x + y
      a = a - b
    end
  end
  x
end



n = 3801307799
(1..n).each do |t|
  s = inv(t, n)
  puts("#{t} * #{s} = #{t * s % n} (mod #{n})")
end