def modinv(a, m) # compute a^-1 mod m if possible
raise "NO INVERSE - #{a} and #{m} not coprime" unless a.gcd(m) == 1
return m if m == 1
m0, inv, x0 = m, 1, 0
while a > 1
inv -= (a / m) * x0
a, m = m, a % m
inv, x0 = x0, inv
end
inv += m0 if inv < 0
inv