
2ヶ月もこのブログを更新していない.忙しかった.問題のほうはさすがにここに来て解くペースがだいぶ落ちてきた.1日1問でも厳しい日が続いている.でも諦めたくない.
さて,上の計算を見てみよう.2の3乗の4乗の5乗の下一桁はなんだ,との問いかけだけど,答えられる? ちなみに,数の大きさは対数を取って見積もったら,489桁であることが分かる.
べき乗のべき乗のべき乗の...で表した巨大整数の剰余計算は今回のテーマ.この種の一般化した問題はどういうわけか,正面から答える書物は少ない.誰でも思いつく問題なのに,計算の仕方を言及する書物はあまりみかけない.だからプログラミングコンテスト問題として出されたのかもしれない.
問題の正確な仕様は ここ を参照すること.つまり,与えられた整数 a1, a2, a3, … , aN および 法 m に対し,次の式を計算すること.
a1^a2^a3^…^aN mod m
なお,各変数の取る範囲は 1 <= ai <= 1000, 1 <= N <= 10, 2 <= m <= 10000.
では,まず重要な定理を紹介しよう.
<離散対数定理 (Discrete logarithm theorem)>
gとmとが互いに素であれば,
gx ≡ gy (mod m)
が成り立つための必要十分条件は x ≡ y (mod φ(m)) である.なお,φ(m) はオイラー関数 (mと素であり,しかも m より小さい自然数の数を表す.) </離散対数定理>
離散対数定理に従えば,g と m が互いに素という条件を満足さえすれば,べき乗の高さは1段下げることができる.例えば,
275 (mod 5)
を計算するときに,値を 2x と設定すれば,2と法5は互いに素なので,
2x ≡ 275 (mod 5)
x ≡ 75 (mod 4) ← φ(5)=4
というふうになり,2段あるべき乗のところが1段に下がり,電卓でもつぎのように値を計算することができる.つまり,
75 = 16807
4で割った時のあまりは 3. ということで,x ≡ 3.
さらに,x ≡ 3 をもとの式のところに戻すと, 275 ≡ 2x ≡ 23 ≡ 8 ≡ 3 (mod 5).つまり,計算結果は 3 である.
このように,再帰手法を使えば,べき乗の高さは10段であろうと,繰り返し処理し最終的にべき乗の高さが1段までに下がり,計算できる訳だ.
さて,べき乗の底 (つまり g )と法 m とが互いに素でないときに,どう計算すれば良いのか.整数論関係の書物はこのことについてあまり言及しない.めんどくさいから書かないか,一般論としての議論は難しいか,理由は分からない.しかし,プログラミング問題は容赦なく,われわれの嫌がる所を聞いてくる.人を困らせるのがプログラミング問題の目的だから?
法が合成数の場合,問題解決に良く役立つのは中国人剰余定理 (Chinese remainder theorem) というもの.ここでもこの定理を使ってみよう.
<中国人剰余定理 (Chinese remainder theorem)>
連立1次合同式
x ≡ b1 (mod m1)
x ≡ b2 (mod m2)
…
x ≡ bn (mod mn)
の場合、 m1, m2, … mn が互いに素であれば、
m = m1m2…mn
を法として、ただ一つの解がある。 また、
(m/mi) xi ≡ 1 (mod mi)
の解 xi を求めることができ,
x ≡ b1 (m/m1) x1 + b2 (m/m2) x2 + … + bn (m/mn) xn (mod m)
とすれば、この x はもとの合同式すべてを満足する解となる。</中国人剰余定理>
この中国人剰余定理に従い,問題の法 m をいかに分割するかをつぎのように考える.
法 m を 互いに素な m1 と m2 との積に分け,m1 を底 g とも互いに素であるようにする.つまり,
m = m1 m2
m1 ⇔ m2 互いに素.このことを整数論では (m1, m2) = 1 と表記されることもある.
m1 ⇔ g 互いに素,つまり, (m1, g) = 1
とする.それで,
gx (mod m1) については, g と m1 は互いに素なので,離散対数定理により計算する.
gx (mod m2) については,以下のように考える.
底 g と法 m2 との最大公約数 gcd(g, m2) を d (d > 1 であることは明らか) とすると,gx は dx で割り切れる.もし,dx は m2 で割り切れるなら,gx も m2 で割り切れることになる.つまり,
gx ≡ 0 (mod m2)
となる.では,dx は m2 で割り切れるだろうか.
m2の大きさは本問題の仕様から分かるように,5000以下であり,素因数 d を含むならば,d の数は最大でも 12 である.なぜなら,
m2 < m/2 = 5000 < 8192 = 213 <= d13
という不等式が成立するからである.一方,gx のほうは, x の値がほとんどの場合 12 以上であろう.勿論,このことはプログラムの中できちんと確認しなければいけないが,x の値が 12 以上と分かれば,gx は 法 m2 で割り切れることになる.
これらの議論をつぎのようにまとめる.
剰余計算 gx (mod m) に対し,
1) x の値が12以下であれば,直接,剰余計算を行い,終了.
2) g と m とが互いに素であるかを確かめる,素であれば,離散対数定理により計算し終了.
3) 互いに素でなければ,g と m との共有素因子をすべて m から取り除き,残った値を m1 とする.
4) gx (mod m1) を離散対数定理により計算し,値を b1 とする.
5) 中国人剰余定理により,gx ≡ b1(m/m1)k (mod m) で最終結果を算出する.ただし, k は (m/m1)k ≡ 1 (mod m1) より得る.
上の計算は勿論べき乗の各高さに繰り返し適応すべきである.
では,冒頭の計算を例としてここで取り上げる.
2345 (mod 10) では, 2 と 10 が互いに素でないので,法 10 を 5 と 2 に分割する.
⇒ 2345 (mod 5) を離散対数定理により計算する.つまり, 345 (mod 4) を計算する.
⇒ 345 (mod 4) では 3 と 法 4 は互いに素であるので,続いて,また離散対数定理を適応する.つまり,45 (mod 2) を計算する.なお,法 2 はオイラー関数 φ(4)=2 とのことから.
⇒ 45 (mod 2) では,5 が12より小さいので,直接剰余計算を行い,答は勿論 0 である. 2進数表現では,偶数の下1桁は常に 0 なので.
⇒ 値を戻し,345≡30≡1 (mod 4) となり,
⇒ 2345≡21≡2 (mod 5) となる.
さらに,中国人剰余定理により,もとの計算 2345≡2*(10/5)*3≡12≡2 (mod 10) となる.
つまり,下一桁は 2 である.なお,3 は (10/5)k≡1 (mod 5) の解 k = 3 による.
以下は本問題を解くC言語プログラム.きちんと整理してなくて,見苦しいところもあろうが,整数論の多くの知識が詰みこまれていて,勉強にはなるだろう.
C Source Solution for UVa P10692 (Huge Mod) Online Judge
少し説明すると,
剰余計算関数 int mod(int a, int p, int m) は,ap (mod m) を直接高速計算する.
最大公約数計算関数 int gcd(int a, int b) は整数 a と b との最大公約数を算出する.
素因数分解関数 int prime_factor(int n) は整数 n (n > 1) のすべての異なる素因数を算出し,外部配列 factor[10] に格納する.異なる素因数の数は関数のリターン値に入っている.なお,高速化するため,素因数 2 と 3 は関数の中では,特別に処理している.
オイラー関数の計算 int phi(int n) (n > 1) は,上の素因数分解関数を利用して,オイラー関数 φ(n) の値を算出する.
合同式を解く関数 int congruence(int a, int m) は,a*x ≡1 (mod m) を解き, x の値をリターンする.
本問題の最重要関数 int huge_mod(int a0, int *a, int size, int m) は a0^a[0]^a[1]^a[2]^…^a[size-1] (mod m) を計算する関数.高速化するための工夫箇所もあるが,上述したアルゴリズムが含まれている.
メイン関数 int main(void) では入出力を司る部分.工夫として,1以降の入力をすべて打ち切るところにある.1のべき乗はつねに1だからである.
プログラミングコンテストの会場では,これらのソースリスト(整数論の基礎知識に関する部分)が手元になければ,数時間で本問題を解くことは大変困難であろう.
ちなみに,本プログラムの計算時間は 0.002秒とUVa判定プログラムが 報告 している.ゼロ秒には到達できていないが,遅い計算法ではなかろう.