モジュラ逆數[逆数]や冪剩餘[剰余]を求める (Ruby)

ほんたうはマクロ無しのExcelとGoogle Spreadsheetで求めたいのだが、とりあへずRubyでやる。

追記 2014-12-04
Google Spreadsheetでやった。Google Apps Scriptを使はないのを前提としてゐる。
cf. RSA模型 https://docs.google.com/spreadsheets/d/1DwIFrzkfVf3xQkRqQR6po33aQELyttZuM2-pFGYzwUE/edit?usp=sharing

冪剩餘[剰余]

まず冪剩餘からやる。ひとつ、再歸[再帰]でひとつづつ掛けてゆく。

RubyVM::InstructionSequence.compile_option = {
  tailcall_optimization: true,
  trace_instruction:     false,
}

eval <<RUBY
def mod_exp b, e, m, accm = 1
  return accm if b == 0 || e == 0
  mod_exp b, e - 1, m, (accm * b) % m
end
RUBY

n, m = 6753, 111112
p mod_exp(n, m, m) # => 45768

そのままやると呼び出しstackが溢れてしまふので、末尾呼び出し最適化を有効にしてやる。RubyでClojure風loop-recur (末尾再帰) http://c4se.hatenablog.com/entry/2013/09/05/184937 でも調べたが、compilerオプションで有効にできる。それを実行時に指定し、evalしてやる。

ふたつ、モンゴメリ乗算をつかふ。

def mod_exp b, e, m, accm = b
  return accm if b == 0 || e == 0
  accm <<= 1
  accm *= accm if e % 2 == 1
  mod_exp b, e >> 1, m, accm % m
end

n, m = 6753, 111112
p mod_exp(n, m, m) # => 45768

かうすれば再歸が圧倒的に浅くなるので、まずstack overflowにはならない。とうぜん圧倒的に速い。モンゴメリ乗算はO(log(n)) の String.prototype.repeat (JavaScript) http://c4se.hatenablog.com/entry/2014/03/02/003653 でもやった。

モジュラ逆數[逆数]

まず拡張したユークリッドの互除法を行なふ。その前にユークリッドの互除法で最大公約数を求めるのははかうやる。

def gcd n, m
  return m if n % m == 0
  gcd m, n % m
end

n, m = 6753, 111112
p gcd(n, m) # => 1

簡単だ。では拡張したユークリッドの互除法で、x * n + y * m = gcd(n, m)となるxとyを求める。

def e_gcd n, m, x = 0, last_x = 1, y = 1, last_y = 0
  return [n, last_x, last_y] if m == 0
  q = n / m
  e_gcd m, n % m, last_x - q * x, x, last_y - q * y, y
end

n, m = 6753, 111112
d, x, y = e_gcd n, m
puts "#{x} * n + #{y} * m = #{x * n + y * m} = #{d}"
  # => 7289 * n + -443 * m = 1 = 1

これを今後つかふ。

モジュラ逆數を算出する。或るmod m (mの剩餘[剰余]類) に於いて、n * k = 1 mod mとなるkを求める。

def e_gcd n, m, x = 0, last_x = 1, y = 1, last_y = 0
  return [n, last_x, last_y] if m == 0
  q = n / m
  e_gcd m, n % m, last_x - q * x, x, last_y - q * y, y
end

def mod_inv n, m
  d, x, y = e_gcd n, m
  raise "#{n} & #{m} aren't coprime." if d != 1
  x % m
end

n, m = 6753, 111112
k = mod_inv n, m
puts "#{n} * #{k} ≡ #{(k * n) % m} mod #{m}"
  # => 6753 * 7289 ≡ 1 mod 111112