Problem 188

正整数aのbによる"超冪",もしくは"テトレーション"をa↑↑bまたはbaと表現し,これを次のように再帰的に定義する.

a↑↑1 = a
a↑↑(k+1) = a(a↑↑k)

よって例えば3↑↑2=33=27,3↑↑3=327=7625597484987,そして3↑↑4はおおよそ103.6383346400240996*10^12という値になる.

1777↑↑1855の下8桁を求めよ.

解説

テトレーション(tetration)は自らの冪を指定された回数だけ反復する演算で,通常の冪を3^3=3↑3のように上向き矢印で表記すると(初期の一部のコンピュータでは上向き矢印を冪乗演算子に使っていたらしい),

  • 3↑↑3 = 3↑3↑3
  • 3↑↑4 = 3↑3↑3↑3

といった具合に展開できる.注意しなければならないのは,展開された冪は右から計算されるということ.つまり,3↑3↑3は(3↑3)↑3 = 27^3ではなく,3↑(3↑3) = 3^27 = 7625597484987だということ(参考:テトレーション - Wikipedia).
つまり,問題の1777↑↑1855について例えば

t <- 1777
for(i in 1:1855) t <- t^1777 mod 1e8

のようなことをしても,これは左から計算していることになってしまうので正しい答は得られない.しかし,値を保持しつつ右側から計算していくとあっという間に値が爆発して,多倍長整数を使ってもまずメモリが足りなくなる.
数論における基本的な定理にオイラーの定理というものがある.これは,nを正整数,aとnを互いに素な正の整数としたとき,

  • aφ(n) ≡ 1 (mod n)

が成立するというもの.ここでφ(n)はオイラーのφ関数(オイラーのφ関数 - Wikipedia)と呼ばれるもので,n以下の自然数のうちnと互いに素な自然数の個数を返す(e.g. φ(6) = 2).
1777と1e8は互いに素なので,

  • 1777^(φ(1e8)) ≡ 1 (mod 1e8)

が成り立つ.オイラーのφ関数は数nの素因数分解がn = Πp^kと与えられているとき(pは素因数,kはその指数),φ(n) = nΠ(1-1/p)として計算できる.
1e8の場合,素因数分解は1e8 = 10^8 = 2^8 * 5^8というように簡単に求められ,φ(1e8) = 1e8 * (1 - 1/2) * (1 - 1/5) = 4e7となる.
よって,展開したテトレーションの右側から冪を計算していく場合,今回は4e7を法とする剰余だけを保持していけばいいことになる.

コード

library(gmp)
## 1e8 = 2^8 + 5^8
## i.e. phi(1e8) = 2^8*(1-1/2) * 5^8*(1-1/5)
phi <- 1e8 * (1-1/2) * (1-1/5)

tetratemod8 <- function(b, t){
  v <- as.bigz(1)
  modulus(v) <- phi
  while(t > 0){
    if(t == 1) modulus(v) <- 1e8
    v <- pow.bigz(b, v) 
    t <- t - 1
  }
  return(v)
}
tetratemod8(1777, 1855)

55ms.
modulus()はbigzクラスのオブジェクトに法を設定する関数で,法を設定するとそれ以降値は剰余になる.次の例を見れば動作は分かるだろう.

> a <- as.bigz(10)
> a
[1] "10"
> modulus(a) <- 7
> a
[1] "(3 %% 7)"
> a + 4
[1] "(0 %% 7)"