Last Digits

別のネタを先に書こうと思っていたのだが、
リクエストがあったので、こっちを先に書くことにします。
問題はこちら

http://acm.pku.cn/JudgeOnline/showproblem?problem_id=2720

問題の要求はシンプルで、
整数入力 b i n (1<=b,i<=100, 1<=n<=7) が与えられたとき、
b^b^b^ ... ^b (つまりbのi回乗) の下n桁を求めよというもの。
Haskell的に書くと

f b 0 = 1
f b i = b ^ f b (i-1)

g b i n = f b i `mod` 10^n

なるgを求めよというものだ。
当然ながらこんなコードをそのまま書いていては
計算がいつまでたっても終わらない
…というかそもそも桁数が宇宙の粒子数をゆうに越えてしまうような気がする。
そこで、intに収まる範囲で下位桁だけ計算してやるということが必要になるが
これがなかなかの曲者である。

そこで、ここではとりあえずだまされたと思って
次のコードを見て頂きたい。

f _ 0 _ = 1
f b i mask = powi b (f b (i-1) mask) mask -- mask は 10^n

powi a b mask :: Int -> Int -> Int -> Int -- a^b `mod` mask を計算する

ご存知の通り累乗の計算というのはlognのオーダで計算できるので、
powiはbがどんな数であろうとほぼ一瞬で計算できる。
そのため、f b n mask も短時間で計算することが出来る。
このfは元のfを各計算ごとに下n桁の計算のみをするというものであり、
元と比べてもかなり単純な処理のままである。
これで答えが計算できれば素敵でしょう？

しかし、残念ながらというか当然ながらというか
この関数では正しい答えを計算できない。
というのもf b i mask の値をmaskで剰余をとってしまっているので、
要するに、a^b (mod n) と、a^(b+10^n) (mod 10^n)が等しいとは限らないので、
結果的に答えが合うとも限らないのである。
実際にf 10 3 1000 を計算してみると、

f 10 3 1000
= powi 10 (f 10 2 1000) 1000
= powi 10 (powi 10 (f 10 1 1000) 1000) 1000
= powi 10 (powi 10 (powi 10 (f 10 0 1000) 1000) 1000) 1000
= powi 10 (powi 10 (powi 10 1 1000) 1000) 1000
= powi 10 (powi 10 10 1000) 1000
= powi 10 0 1000 -- powi 10 10 1000 が 0 になる
= 1

これは0になるべきなのでおかしくなる。
しかし、ここであきらめてはいけない。
実は次に示すような性質がこの式にはあるのだ。

定理1：累乗の下位桁の周期に関する定理

∃a, ∀b>=a, ∀n>=2, ∀x, x^b ≡ x^(b+10^n) (mod 10^n)

定理2：累乗の下位桁の周期の下限に関する定理(1)

∀b, ∀n>=2, ∀x, x^(b+10^n) ≡ x^(b+2*10^n) (mod 10^n)

定理3：累乗の下位桁の周期の下限に関する定理(2)

∀b,x^(b-5)>=10^n, ∀n>=2, ∀x, x^b ≡ x^(b+10^n) (mod 10^n)

定理1は任意のxの累乗の下位桁nに対して、10^nの周期が存在するという
かなり大胆な定理…というか予想だ。
周期の最大値は明らかに10^nなので、
任意の周期はその約数になるということか。
n>=2となっているのは、n=1には簡単な反例があり、成り立たないからである。
定理2は定理1が正しいなら即座に導くことが出来る。
定理1におけるaは10^n以下として良いということがいえる。
定理3は定理2より強力な周期の下限の存在(というかこれも予想…)定理である。
"-5"というマジックナンバーがいかにも怪しいが、
正確な下限は良く分からないので、これで勘弁。
少なくともそんなに大きくは無いはずだということ。
定理1のaがかなり小さくても成り立つということがいえる。

証明であるが、ちょっと数学的には無理だった…ので
四色定理よろしく計算機にさせてみることにした。
といっても完全に証明するのは無理なので、
とりあえず今回必要な範囲に限って証明することにした。
x^b>=10^n,x^(b-1)<10^n なるx,b,nに対して
x^(b+5) `mod` 10^n = x^(b+5+10^n) `mod` 10^n であれば
定理1〜3はすべて証明されるので(これは…自明ですよね)
そのようなコードを作成すれば良い。

import Control.Monad

main = do
  res <- mapM proofProcess $ [(a,b) | a <- [2..100000], b <- [2..7]]
  let ok = and res
  putStrLn $ "the theorem is " ++ if ok then "TRUTH" else "FALSE"

proofProcess (x,n) = do
  let res = proof x n
  when (not res) $ putStrLn $ "x = " ++ show x ++ " , n = " ++ show n ++ " : NG"
  return res

proof x n = powi x b mask == powi x (b+mask) mask where
  mask = 10^n
  b = searchb 0
  searchb b | x^b >= mask = b+5
            | otherwise = searchb (b+1)

powi _ 0 _ = 1
powi a b mask =
  let t = powi a (b `div` 2) mask
      u = (t*t) `mod` mask
      v = if b `mod` 2 == 1 then (u*a) `mod` mask else u
  in v

とりあえずx<=100の範囲で証明が出来ればいいのだが、
せっかくなので100000まで証明してみた。
実際に実行して確かめていただきたい。
というわけで、上の定理はx<=100000,n<=7という条件を付けていただくと
正しいものだということができる。

さて、この定理を用いると、f b n mask の定義は
ごく一部のパラメータ以外に対してはほとんどうまく動くことが分かるだろう。
動かないものに対して補正をしてやるなりすれば正しいプログラムが出来る。
それを書くのは大した手間では無いと思うので、ここでは実装の詳細までは触れない。

• • • • -

☆あらかじめお断り☆
私には整数論の素養などほとんど無いので、
上の定理は一般には全くのでたらめであるかもしれません。
もし数学的に証明できた/誤りが証明できた、あるいは
こういう定理はとっくの昔に発見されているなどありましたら
情報をお寄せいただければ幸いです。