大きな数の先頭と末尾（3）

「大きな数の先頭と末尾（2）」を awk で書いてみると，「大きな数の先頭と末尾（2）」が R の特殊な機能を使っているのではないことが理解できるだろう。

# file name : func.awk
BEGIN {
    print fun(3863080011, 2613515386) # 254361
    print fun(21321331, 1234653876) # 81
    print fun(521340345, 1396720193) # 84765625
    print fun(1843165372, 2835135645) # 79967232
    print fun(417934669, 3961963772) # 11298161
    print fun(7564929, 1434531250) # 1
    print fun(3713420107, 616334320) # 57768001
    print fun(57564748, 243756997) # 3328
    print fun(2249407224, 3483719867) # 32397824
    print fun(444893257, 3572472608) # 25452801
}

function mul(m, n,   km, kn, am, an, i, j, carry, num, ans) { # 8 桁以下の，m, n の掛け算を行い，下 8 桁を返す
    km = m
    kn = n
    for (i = 1; i <= 8; i++) {
        am[i] = km % 10
        km = int(km / 10)
        an[i] = kn % 10
        kn = int(kn / 10)
    }
    for (i = 1; i <= 8; i++) {
        for (j = 1; j <= 8; j++) {
            ans[i+j-1] += am[i]*an[j]
        }
    }
    carry = 0 # 繰り上がり
    num = 0
    for (i = 1; i <= 8; i++) { # 繰り上がり処理
        ans[i] += carry
        carry = int(ans[i]  / 10)
        num += (ans[i] % 10) * 10^(i-1)
    }
    return num # 答の下 8 桁を返す
}

function fun(m, n,   ans) {
    m = m % 1e8
    ans = 1
    for (;;) {
        if (n % 2 == 1) {
            ans = mul(ans, m) # もとの n の ビットが立っているところを掛け算
        }
        n = int(n / 2) # もとの n のビットが立っているところを探すため半分ずつにしていく
        if (n == 0) return ans
        m = mul(m, m) # 被巾数は倍々にしていく
    }
    return ans
}

$ gawk -f func.awk
254361
81
84765625
79967232
11298161
1
57768001
3328
32397824
25452801

大きな数の先頭と末尾（2）

偶然，「大きな数の先頭と末尾」の後に類似問題が出たようだ。

「大きな数の先頭と末尾」では本当は先頭1桁だけが必要なのだが，4 桁求めるようにしていた。

更に，べき乗される数もそんなに大きいものは考えていなかった。

なので，仕様を満たすために，以下のプログラムを書くことになった。

問題：a, x を符号なし 32 bit 整数として，べき乗計算 a^x の下 8 桁のみを高速に計算するプログラム

mul = function(m, n) { # 8 桁以下の，m, n の掛け算を行い，下 8 桁を返す
  dec = function(k) { # 8 桁の整数を一桁ずつ配列に落とす
    a = numeric(8)
    for (i in 1:8) {
      a[i] = k %% 10
      k = k %/% 10
    }
    return(a)
  }
  am = dec(m)
  an = dec(n)
  ans = numeric(17) # 掛け算の答が入る
  for (i in seq_along(am)) {
    for (j in seq_along(an)) {
      ans[i+j-1] = ans[i+j-1]+am[i]*an[j]
    }
  }
  carry = 0 # 繰り上がり
  for (i in 1:8) { # 繰り上がり処理
    ans[i] = ans[i]+carry
    carry = ans[i]  %/% 10
    ans[i] = ans[i] %% 10
  }
  return(sum(ans[1:8] * 10^(0:7))) # 答の下 8 桁を返す
}
 
func = function(m, n) {
  m = m %% 1e8
  ans = 1
  repeat {
    if (n %% 2 == 1) {
      ans = mul(ans, m) # もとの n の ビットが立っているところを掛け算
    }
    n = n %/% 2 # もとの n のビットが立っているところを探すため半分ずつにしていく
    if (n == 0) return(ans)
    m = mul(m, m) # 被巾数は倍々にしていく
  }
  ans
}

func(3863080011, 2613515386) # 254361
func(21321331, 1234653876) # 81
func(521340345, 1396720193) # 84765625
func(1843165372, 2835135645) # 79967232
func(417934669, 3961963772) # 11298161
func(7564929, 1434531250) # 1
func(3713420107, 616334320) # 57768001
func(57564748, 243756997) # 3328
func(2249407224, 3483719867) # 32397824
func(444893257, 3572472608) # 25452801

以上全てを計算するのに，小生のポンコツ・マシンは 0.126 秒を要した。

 

一休さんのとんち話

一休さんのとんち話にもあったが，「一日目は1粒，二日目は倍の2粒，三日目は更にその倍の4粒というように増やして行くとき，一日に 8000万粒以上になるのは何日目か」というもの。

プログラムしてはいけない。ただ，そのためには log10(2)≒0.3010 を記憶していなければならないが（大学受験生なら知っているだろう）。

80000000 ≦ 2^(n-1)
3+7/log10(2) ≦ n-1
27.3 ≦ n

よって，28日目
なお，27日目までの累計も 80000000 を超える。
sum(2^0 + 2^1 + 2^2 + … + 2^(n-1)) = 2^n - 1