znorder.rb

require "prime"

# 参考: PARI/GPのソースコード
# http://pari.math.u-bordeaux.fr/
# http://pari.math.u-bordeaux.fr/cgi-bin/gitweb.cgi?p=pari.git;a=tree
# src/basemath/arith1.cのznorder()関数
# src/basemath/gen2.cのZ_pval()関数

# ord_m(a)を求める
# すなわちa^n ≡ 1 (mod m)となる最小のn≧1
def znorder(a, m)
  factors = Prime.prime_division(m)
  o = 1
  factors.each do |p, e|
    o = lcm(o, zp_order(a, p, e))
  end
  o
end

# ord_{p^e}(a)を求める
# すなわちa^n ≡ 1 (mod p^e)となる最小のn≧1
#
# 次の定理を利用している:
#  1 + v_p(2) <= v_p(a - 1) <= e ならば ord_{p^e}(a) = p^{e - v_p(a - 1)}
def zp_order(a, p, e)
  pe = p ** e
  if p == 2
    if e == 1
      return 1
    elsif e == 2
      return a % 4 == 1 ? 1 : 2
    end
    if a % 4 == 1
      op = 1
    else
      op = 2
      a = (a * a) % pe
    end
  else
    op = fp_order(a % p, p)
    if e == 1 then return op end
    a = powmod(a, op, pe)
  end
  if a == 1
    op
  else
    op * (p ** (e - pval(a - 1, p)))
  end
end

# ord_p(a)を求める
# すなわちa^n ≡ 1 (mod p)となる最小のn≧1
#
# a^{p-1} ≡ 1 (mod p)となる(フェルマーの小定理)のでp-1を素因数分解しp-1の約数から答えを見つける
def fp_order(a, p)
  o = p - 1
  factors = Prime.prime_division(o)
  factors.reverse_each do |l, e|
    t = o / (l ** e)
    y = powmod(a, t, p)
    if y == 1
      o = t
    else
      j = (1..e).find {
        y = powmod(y, l, p)
        y == 1
      }
      o = t * (l ** j)
    end
  end
  o
end

# (a^n) mod m
def powmod(a, n, m)
  if n == 0
    1
  elsif n.even?
    powmod((a*a) % m, n / 2, m)
  else
    (powmod((a*a) % m, n / 2, m) * a) % m
  end
end

# p進付値
def pval(x, p)
  (n, y) = lval(x, p)
  n
end

# x = q^i * y (yはqで割り切れない) と表示するときの(i, y)を返す
def lval(x, q)
  raise ArgumentError, "x == 0" if x == 0
  (z, r) = x.divmod(q)
  if r != 0
    [0, x]
  else
    (v, y) = lval(z, q*q)
    (zz, rr) = y.divmod(q)
    if rr != 0
      [2 * v + 1, y]
    else
      [2 * v + 2, zz]
    end
  end
end

def lcm(a, b)
  a * b / gcd(a, b)
end

def gcd(a, b)
  if b == 0
    a
  else
    gcd(b, a % b)
  end
end