余りが１になる割り算から元の数を求める方法

1000003 × B ÷ 1000000 = q 余り 1 となる最小のBは何か？

これを素朴に求めると、B＝1から開始して初めて余りが 1になるまでBをインクリメントするような事になるが、ずっと効率的なやり方があったのでメモっておく。

以下、euler 134 のネタバレあり

==== ==== ==== ====

project euler に、与えられた素数 p1, p2 について、末尾が p1 で、なおかつ p2 で割りきれる数 S をもの凄く大量に計算する問題があった。

合同式で書くとこうなる

S ≡ 0 (mod p2)
S ≡ p1 (mod m) (m はp1を越える最小の10の正数乗)

実は合同式自体あまり詳しくないので、"連立 合同"でググってみると、連立合同式の解法を紹介しているサイトがあった。

これをよく読んで、上の問題に当てはめてみると、

  S ≡ p1*p2*B (mod p2*m)
  ただし、B は p2*B ≡ 1 (mod m) を満たす整数)

を解けば良い事がわかる。

試しにやってみると、

p1 = 19, p2 = 23
23*B ≡ 1 (mod 100) より、B = 87
(19 × 23 × 87) ÷ (23 × 100) = 16 余り 1219
答え： 1219

となり、一応計算できてることが分かる。

ただし困ったのが、B を求めるやり方で、素朴にやると以下のようになるけど（Haskell）

ghci> find (\a ->1==mod (23*a) 100) [1..]
Just 87

これだと冒頭であげた1000000近辺の数だとかなり時間がかかる。

そこで更に調べてみると、「逆元（inverse）」なるものを求める方法があって、どうやらこれが使えそう。

ググりまわって Wikipedia でこの擬似コードを見つけた

remainder[1] := f(x)
remainder[2] := a(x)
auxiliary[1] := 0
auxiliary[2] := 1
i := 2
while remainder[i] > 1
    i := i + 1
    remainder[i] := remainder(remainder[i-2] / remainder[i-1])
    quotient[i] := quotient(remainder[i-2] / remainder[i-1])
    auxiliary[i] := -quotient[i] * auxiliary[i-1] + auxiliary[i-2]
inverse := auxiliary[i]

これを、以下のように解釈して、表計算ソフトに打ち込んでみる。

n	r	a	q
1	m	0	-
2	p2	1	-
・・・
i	MOD(ri-2, ri-1)	= - qi * ai-1 + ai-2	= DIV( ri-2, ri-1)

で、r が1 になったときの a が求める数字らしい。

試しに 1000003 でやってみると以下のようになる。

1000000	0
1000003	1
1000000	0	0
3	1	1
1	-333333	333333

-333333 が負数になっているので、mod (-333333,1000000)としてみると 666667 が得られる。 1000003 × 666667 ＝ 666669000001で、確かに 1000000 で割ると余りが 1 になる事が分かる。

もの凄く計算量が減ってる。ははは、これは面白い。

上の問題に戻って、p1 = 999983, p2 = 1000003 としてみると、mod (999983 * 1000003 * 666667) (1000003000000) = 666662999983 となり、1000003 で割りきれて末尾が 999983 になる数が確かに得られた。

上手く行きそうなので、表計算ソフトに書いた事を以下のような Haskell コードに置き換えてみる。

inverse a m = mod (f m 0 0 a 1 0) m where
    f 1  a  _  _  _  _  = a
    f r1 a1 q1 r2 a2 q2 = f r a q r1 a1 q1
      where r = mod r2 r1
            q = div r2 r1
            a = -q*a1 + a2

euler 134 に適用してみると、1秒で正答が得られた。（ちなみに「余り１になるB」をインクリメントしながら計算すると、4時間くらいかかった。）