Help us understand the problem. What is going on with this article?

正則連分数を用いたモジュラ逆数の計算

最近,拡張ユークリッド互除法を用いたモジュラ逆数の計算を学んで,正則連分数と関連付けられそうだなぁ,と思ったので書きました.実装はHaskellです.

問題

素数 $p$ と,$2$ 以上 $p$ 未満の整数 $m$ が与えられる.$mx\equiv1\mod p$ を満たす整数 $x$ を求めよ.ただし $0\lt x\lt p$ とする.

このような $x$ は存在して一意に定まりますが,証明は述べません.$x$ は法 $p$ に関する $m$ のモジュラ逆数と呼ばれます.
条件である $mx\equiv1\mod p$ は,「ある非負整数 $q$ が存在して $mx-pq=1$ 」と言い換えることができます.$q$ が非負整数であることは,$a,x,p$が正整数であることからわかります.以下では,この形で条件を扱います.

理論

正則連分数

有理数1 $r$ の正則連分数による表現は次のようなものです.
$$r = a_0 + \frac1{\displaystyle a_1+\frac1{\displaystyle a_2+\frac1{\displaystyle \cdots a_{n-1}+\frac1{a_n}}}}$$
このとき,$r=[a_0;a_1,a_2,\cdots,a_n]$と書きます.「$\cdots$」に無茶がある気がするので,再帰的な定義も記します.
$$[a_0;] = a_0$$$$[a_0;a_1,a_2,\cdots,a_n] = a_0 + \frac1{[a_1;a_2,a_3,\cdots,a_n]}$$
さて,次のように数列$(q_n),(x_n)$を定めます.
$$q_0=0\quad q_1=x_0 = 1\quad x_1 = a_0$$$$q_n=a_{n-1}q_{n-1}+q_{n-2}\quad x_n=a_{n-1}x_{n-1}+x_{n-2}$$すると,$n$ についての帰納法によって
$$[a_0;a_1,a_2,\cdots,a_{n-1}]=\frac{x_n}{q_n}$$となることがわかります.証明は省略します.2

正則連分数とモジュラ逆数の関係

結論としては,上記問題設定のもとで,
$$\frac pm = [a_0;a_1,a_2,\cdots,a_n]$$のとき,
$$\frac {x^\prime}{q^\prime} = [a_0;a_1,a_2,\cdots,a_{n-1}]$$とおけば,$x^\prime,q^\prime$ はそれぞれ $p,m$ より小さく,3
$$mx^\prime - pq^\prime = (-1)^n$$が成り立ちます.
$n$ が偶数の場合には $x=x^\prime,\ q=q^\prime$ です.$n$ が奇数であれば,両辺に$-1$をかけて,$x=p-x^\prime,\ q=m-q^\prime$ とすればよいです.よって問題が解けました.4

証明

$n$ についての帰納法によって証明します.$m\ge2$ で $p$ が素数であることから,$n\ge1$は容易に確認できます.5

$n=1$ の場合は,
$$\frac pm = [a_0; a_1] = a_0 + \frac1{a_1} = \frac{p-1}m + \frac1m$$となります.6 すなわち,
$$\frac{x^\prime}{q^\prime} = [a_0;] = a_0 = \frac{p-1}m$$より,$x^\prime = p-1,\ q^\prime = m$です.このとき,
$$mx^\prime - pq^\prime = m(p-1) - pm = -1 = (-1)^1$$となっています.

$n=k-1\ge1$ で成り立つとします.また,
$$[a_0;a_1,a_2,\cdots,a_{k-1}]=\frac{x^{\prime\prime}}{q^{\prime\prime}}$$とおきます.すなわち,
$$q^\prime x^{\prime\prime}-x^\prime q^{\prime\prime} = (-1)^{k-1}$$を仮定します.このとき,$n=k$ では
$$\frac pm=\frac{a_{k-1}x^\prime+x^{\prime\prime}}{a_{k-1}q^\prime+q^{\prime\prime}}$$から,

\begin{align}
mx^\prime-pq^\prime
&= (a_{k-1}q^\prime+q^{\prime\prime})x^\prime - (a_{k-1}x^\prime+x^{\prime\prime})q^\prime \\
&= -(q^\prime x^{\prime\prime}-x^\prime q^{\prime\prime}) = (-1)^k
\end{align}

となります.よって,すべての正整数 $n$ に対して $mx^\prime - pq^\prime = (-1)^n$ が成り立ちます. (証明終わり)

実装

Haskellで実装します.まずは正則連分数への展開から.

import Data.Ratio

type CFrac = [Integer]
toCFrac :: Rational -> CFrac
toCFrac a
  | r == 0 = [q]
  | otherwise = q : toCFrac (d % r)
  where
    (n, d) = numDenom a
    (q, r) = divMod n d

numDenom :: Rational -> (Integer, Integer)
numDenom a = (numerator a, denominator a)

イメージとしては,(あくまでイメージですが)
$$\mathrm{toCFrac}\left(\frac nd\right) = q +
\frac1{\displaystyle\mathrm{toCFrac}\left(\frac dr\right)}$$のような感じです.こんどは正則連分数を有理数に戻します.

import Prelude hiding (toRational)

toRational :: CFrac -> Rational
toRational [a] = a % 1
toRational (a : as) = (a % 1) + recip (toRational as)

これはそのまま
$$\mathrm{toRational}([a_0;a_1,a_2,\cdots,a_n])=a_0+\frac1{\mathrm{toRational}([a_1;a_2,\cdots,a_n])}$$と考えていいですね.7
これらを使って,モジュラ逆数を計算します.

modInv :: Integer -> Integer -> Integer
modInv p m = x where
  (x, q) = if even n then (x', q') else (p - x', m - q')
  (x', q') = numDenom r
  r = toRational $ init cf
  cf = toCFrac $ p % m
  n = length cf - 1

理論で使った文字そのままで,特に面白いところもないですね.

ちょっとだけ面白いことをします.実装から連分数を削除します.

modInv' p m 
  | minus = p - numerator frac
  | otherwise = numerator frac
  where
    (frac, minus) = f p m 
    f p m
      | r == 1 = (a % 1, True)
      | otherwise = ((a % 1) + recip next, not minus)
      where
        (a, r) = divMod p m
        (next, minus) = f m r

ついでにRationalも消してしまいましょう.

modInv'' p m
  | minus = p - x
  | otherwise = x
  where
    (x, q, minus) = f p m
    f p m
      | r == 1 = (a * 1 + 0, 1, True)
      | otherwise = (a * x + q, x, not minus)
      where
        (a, r) = divMod p m
        (x, q, minus) = f m r

ユークリッドの互除法で最後の状態$(1,0)$のひとつ手前$(?,1)$を捨てて逆順にたどっていることになります.


  1. 無理数でも表現可能ですが,$a_0,a_1,\cdots$ の列が無限に長くなります.ここでは有理数のみ扱います. 

  2. $[a_1;a_2,a_3,\cdots,a_n]$ に対して同様に数列 $(q_n^\prime),(x_n^\prime)$ を定め,$q_n=x_{n-1}^\prime,\ x_n=a_0x_{n-1}^\prime+q_{n-1}^\prime$ から証明します. 

  3. 証明していませんが,連分数の性質などから明らかです. 

  4. 分数と連分数の相互変換が計算できることを前提にしています.まぁ,そりゃそう. 

  5. $m$ と $p$ が互いに素なので $\frac pm$ は整数ではありません.よって $n\ge1$ です 

  6. かなり無茶はありますが,正則連分数による表現の定義を厳密にしてやると確認できます.ここでは深く追究しません. 

  7. $\mathrm{toRational}(x)=x$ になってしまっていますが,これはtoRational表現の変換だからです. 

cunitac
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした