8 バイト実数では，普通に階乗を計算しようとすると 171! は計算できない。
> prod(1:170)
[1] 7.257415615308e+306
> prod(1:171)
[1] Inf
> factorial(171)
[1] Inf
警告メッセージ:
In factorial(171) : 'gammafn' 中の値が範囲を超えています
フィボナッチ数列のときと同じく，階乗の先頭一桁の分布がベンフォードの法則に従うかどうかを見たいとき，170 個の数字の分布では説得力がない。
そこで，
問題: 特別なライブラリや関数などを使わずに，10000000! の先頭の数字を求めよ。
プログラム例は，この記事の「コメント」として投稿しておく。
さて，そのようなプログラムを若干修正して 1! ～ 10000000! の先頭の数字 1e7 個の分布を調べる。
度数 3011187 1761416 1249153 969498 789733 669636 578730 512327 458320
パーセント 0.301119 0.176142 0.124915 0.096950 0.078973 0.066964 0.057873 0.051233 0.045832
理論値 0.301030 0.176091 0.124939 0.096910 0.079181 0.066947 0.057992 0.051153 0.045757
フィボナッチ数列の場合は理論値にごくごく近かったが，階乗の場合は少しズレがある。なぜかな？
特別なライブラリや関数を用いずに，8 の 10000000 乗（8 ^ 1e7）の先頭 4 桁と末尾 4 桁を求めよ。
解答例はコメントを参照
フィボナッチ数列は急速に大きくなり，普通に計算するとオーバーフローしてしまう。
8 バイト実数での有効数字 16 桁なので，79 項目は 14472334024676220 と正確に表示できるが，それ以降は上位桁しかわからない。1476 項目は 1.306989e+308 となるが，1477 項目は計算されない（Inf になる）。
ベンフォードの法則を実例で示すために，フィボナッチ数列の各項の先頭の数字を求めて度数分布を求めるとき，1476 項までの結果では面白くない（その先どうなるか分からないじゃないかという意見もあるかも知れないし）。
多倍長演算によれば計算はできる。R だと，gmp ライブラリには fibnum 関数がある。しかし，それを使わなくても，目的は達成できる。
問題: 特別なライブラリや関数を使わずに，フィボナッチ数列の 1e7 項目の先頭の数字を求めよ。もちろん，計算時間は長くても数分以内とする。なお，答えは 1 である。
プログラム例は，この記事の「コメント」として投稿しておく。
さて，そのようなプログラムを使って，「フィボナッチ数列の 1e7 項までの先頭の数字の度数分布を調べる」と以下のようになる。
度数 3010300 1760918 1249382 969101 791817 669467 579913 511531 457571
パーセント 0.301030 0.176092 0.124938 0.096910 0.079182 0.066947 0.057991 0.051153 0.045757
理論値 0.301030 0.176091 0.124939 0.096910 0.079181 0.066947 0.057992 0.051153 0.045757
「難しいパズルを解きたい」と思うのは，「パズルを解くのが面白いから」というのが動機だろう
パズルを解く参加者に制限を加えたり，正解が開示されなかったり，期限（？）が過ぎたら問題自体が参照できなくなったり，最初提示された締め切り日がズルズルズルズルと引き延ばされたり，挙げ句の果てに急に終わってしまう，というような，パズルを解くということと無関係（正反対）な制約があるのでは，なんだかなあと思う。
それらの制約が，「優秀なコンピュータサイエンティストの発掘」，「求職者来たれ!」ということなんだろうけど，その問題自体がその目的を満たすにはほど遠い，低レベル（？）なものでであるとしたら，なにをかいわんや（アハ!）。そんなレベルの求職者にどの程度の待遇を用意するのか？？
今後も，どしどし，解答するぞ!
それが，低レベルなものなら，反面教師になるだろう。
同じ解答を得られる，もっと高精度・高速なプログラムを得た人は，本来の窓口から投稿すべし。
主催者は，ちょっと運営方針を見直した方がよいんじゃないか？？
「サルベジオン社で宇宙船のデータを救え」とのことで...
問題 1. キーは昇順になっているので，二分法で探索する。
求める値は "V406435859539156181269150751031"
library(gmp)
url = "http://salvageon.textfile.org/?db=1&index=%s"
key = as.bigz("208050656559285601386927895421059705239114932023754")
begin = as.bigz("0")
end = as.bigz("1267650600228229401496703205375")
for (i in 0:110) {
cat(i, as.character(end-begin), "\n")
med = (begin+end) %/% 2
m = scan(sprintf(url, med), ""); Sys.sleep(1)
k = as.bigz(substr(m[3], 2, nchar(m[3])))
if (k > key) {
end = med
} else if (key > k) {
begin = med
} else {
print(m)
break
}
}
問題 2. 奇数のキーは後ろ半分にまとまって，順序正しく並んでいる。
求める値は "V1101943557675920722238136981003"
これは，プログラムを組まずに求めた。
まず，レコードの個数が 2^100 なのに気づけば，ちょうど真ん中が何になっているか，気になる。真ん中の次，その次と見てみれば，ははあ～んと気づいた。
2023636070998557444542586045 は 初項 = 1，項差 = 2 の数列の何番目か？
> (as.bigz("2023636070998557444542586045")+1)/2
[1] 1011818035499278722271293023
対応するコードは
> as.bigz("633825300114114700748351602688")+(as.bigz("2023636070998557444542586045")+1)%/%2 - 1
Big Integer ('bigz') :
[1] 634837118149613979470622895710
コードは 634837118149613979470622895710
scan("http://salvageon.textfile.org/?db=2&index=634837118149613979470622895710", "")
問題 2 で，キーが偶数だったら，ちょっと面倒だったはず。
1 円，2 円，10 円，52 円，82 円，280 円の切手を合計が 2007 円になるように買う方法を求めよ。ただし，買う枚数を最小にせよ。
素直に考えて，280 円切手を 7 枚，10 円切手を 4 枚，2 円切手を 3 枚，1 円切手を 1 枚，つまり，280*7 + 10*4 + 2*3 + 1*1 =2007 ということで，15 枚が取りあえずの候補。
280 円切手だけは 2007 %/% 280 = 7 枚までしか買えないがそのほかは，まあ他の制約はあるけど，15 枚以下でないと候補にはならない。そこで，効率なんかは考えずにプログラムをして，以下のプログラムを得る。
fee = 2007
for (i280 in 0:(fee %/% 280)) {
for(i82 in 0:15) {
for (i52 in 0:15) {
for (i10 in 0:15) {
for (i2 in 0:15) {
i1 = fee-(i280*280+i82*82+i52*52+i10*10+i2*2)
if (i1 < 0) break
if (sum(c(i1, i2, i10, i52, i82, i280)) <= 15) {
cat(i1, i2, i10, i52, i82, i280, "\n")
}
}
}
}
}
}
どれくらい計算時間が掛かるか分からないが（いや，たいしたことないことは分かる），実行して見ると，以下の結果が表示される。
1 0 1 2 6 5
1 3 0 3 2 6
1 3 4 0 0 7
つまり
1 円切手 1 枚，2 円切手 0 枚，10 円切手 1 枚，52 円切手 2 枚，82 円切手 6 枚，280 円切手 5 枚。計 15 枚で 2007円。
1 円切手 1 枚，2 円切手 3 枚，10 円切手 0 枚，52 円切手 3 枚，82 円切手 2 枚，280 円切手 6 枚。計 15 枚で 2007円。
1 円切手 1 枚，2 円切手 3 枚，10 円切手 4 枚，52 円切手 0 枚，82 円切手 0 枚，280 円切手 7 枚。計 15 枚で 2007円。
のいずれでもよい。
2 桁の整数のうち，2 進数表記と BCD 表記に含まれる 1 の個数が同じになる整数はいくつあるか。
馬鹿馬鹿しい問題だけど...
# 2 進数表記に含まれる 1 の個数
bin = function(n) {
count = 0
repeat {
count = count + n %% 2
n = n %/% 2
if (n == 0) return(count)
}
}
# 0～9 の BCD 表記に含まれる 1 の個数
count1 = numeric(10)
for (d in 0:9) count1[d+1] = bin(d)
# 2 桁の整数の 2 進数表記と BCD 表記に含まれる 1 の個数が同じかどうか
same = logical(99)
for (n in 11:99) {
same[n] = bin(n) == count1[n %/% 10+1] + count1[n %% 10+1]
}
sum(same) # 同じになる数の個数
(1:99)[same] # そのような数のリスト
実行結果
> sum(same) # 同じになる数の個数
[1] 20
> (1:99)[same] # そのような数のリスト
[1] 12 13 18 19 24 25 26 27 38 39 48 49 52 53 70 71 78 79 98 99
エネルギー最短経路の問題と同じ。ただ，対称行列であるところが違うが。
効率がよいか悪いか分からない新しい言語（？）を覚えるまでもない。
route = matrix(c(
0, 200, 200, 390, 160, 220,
200, 0, 160, 220, 220, 310,
200, 160, 0, 310, 220, 310,
390, 220, 310, 0, 470, 550,
160, 220, 220, 470, 0, 220,
220, 310, 310, 550, 220, 0), ncol=6, byrow=TRUE)
station = c("東京", "新宿", "渋谷", "三鷹", "錦糸町", "北千住")
dimnames(route) = list(station, station)
# 料金表
route
東京 新宿 渋谷 三鷹 錦糸町 北千住
東京 0 200 200 390 160 220
新宿 200 0 160 220 220 310
渋谷 200 160 0 310 220 310
三鷹 390 220 310 0 470 550
錦糸町 160 220 220 470 0 220
北千住 220 310 310 550 220 0
library(e1071)
# 東京から始まり東京で終わる場合
d = permutations(length(station))
d = d[d[,1] == 1,]
fee = apply(d, 1, function(from) sum(route[cbind(from, c(from[-1], from[1]))]))
p = which.min(fee)
fee[p]
apply(d[fee == fee[p],], 1, function(n) station[c(n, 1)])
# 東京から始まり錦糸町で終わる場合
d = permutations(length(station))
d = d[d[,1] == 1 & d[,6] == 5,]
fee = apply(d, 1, function(from) sum(route[cbind(from[-6], from[-1])]))
p = which.min(fee)
fee[p]
apply(d[fee == fee[p],], 1, function(n) station[n])
径路は列方向に見る
東京から始まり東京で終わる場合
> fee[p] # 料金
[1] 1390
> apply(d[fee == fee[p],], 1, function(n) station[c(n, 1)])
[,1] [,2] [,3] [,4]
[1,] "東京" "東京" "東京" "東京"
[2,] "渋谷" "新宿" "北千住" "北千住"
[3,] "三鷹" "三鷹" "錦糸町" "錦糸町"
[4,] "新宿" "渋谷" "渋谷" "新宿"
[5,] "錦糸町" "錦糸町" "三鷹" "三鷹"
[6,] "北千住" "北千住" "新宿" "渋谷"
[7,] "東京" "東京" "東京" "東京"
東京から始まり錦糸町で終わる場合
> fee[p] # 料金
[1] 1260
> apply(d[fee == fee[p],], 1, function(n) station[n])
[,1] [,2]
[1,] "東京" "東京"
[2,] "渋谷" "新宿"
[3,] "三鷹" "三鷹"
[4,] "新宿" "渋谷"
[5,] "北千住" "北千住"
[6,] "錦糸町" "錦糸町"
1～2014 までの数を 2 進数で表したときに，連続する 2 数の 1 の個数が同じであるようなペアは何組存在するか
連続する数のうち，(1, 2)，(5, 6), (9, 10), ...
のように, (4(i-1)+1, 4(i-1)+2)，i = 1, 2, ... の 2 数の 1 の個数が等しくなる。
よって，n までには ifelse((n-1)%%4, (n-1)%/%4+1, (n-1)%/%4) 組存在する
n = 2014 ならば，
> n = 2014
> ifelse((n-1)%%4, (n-1)%/%4+1, (n-1)%/%4)
[1] 504
504 組存在する
愚直なプログラムを書いて確かめるなら，
> func = function(n) {
+ count = 0
+ repeat {
+ count = count + n %% 2
+ n = n %/% 2
+ if (n == 0) return(count)
+ }
+ }
> sum(diff(sapply(1:n, func)) == 0)
[1] 504
20 分で答よということなので，できるだけ無精をすることにすると...
library(gmp)
p = 1
repeat {
p = nextprime(p)
if (p > 100) break
print(as.integer(p))
}
「パスカルの三角形」は「右上の数と左上の数の和」を配置するが，「和」ではなく「排他的論理和」を使う場合を考える。
2014番目の「0」が出力されるのは何段目になるか？
規則性もあるだろうが，たかが 2014 番目なので，多めに見積もって，馬鹿正直なプログラムを書く。
それでも，解を得るのに 0.06 程度しかかからないのだから，これ以上の策を弄するのは無駄というものだ。
system.time({
m = 100
a = matrix(0, m, m+1)
a[1, 2] = 1
for (i in 2:m) {
for (j in 2:(i+1)) {
a[i, j] = xor(a[i-1, j-1], a[i-1, j])
}
}
b = 1:m - rowSums(a)
p = which.min(2014 > cumsum(b))
print(p) # 解
print(sum(b[1:(p-1)])) # 念のための確認
print(sum(b[1:p])) # 念のための確認
})
[1] 75 # 75 段目
[1] 1980 # 74 段目の終わりまでに 1980 個
[1] 2047 # 75 段目の終わりまでに 2047 個
所要時間
ユーザ システム 経過
0.058 0.003 0.061
「どの桁にも 0, 3, 6, 9 が表れない」という条件を満たす正の整数を小さい順に並べて以下のような数列を作る。
1, 2, 4, 5, 7, 8, 11, 12, 14, 15, 17, 18, 21, 22, ...
この数列の 30 番目の項は 58 である。では，この数列の 10^7 番目の項を求めよ。
n 番目の項を手計算で求める。
0, 3, 6, 9 を含まない数だけを考えるということは，
要するに 1, 2, 4, 5, 7, 8 という記号で表される数を，0, 1, 2, 3, 4, 5 という 5 つの数字で表される 6 進数を考えるということ
出現する数字 6 進数字
1 0
2 1
4 2
5 3
7 4
8 5
6 進数字を使って表示すると
1 桁のものは 6 個 0，1, 2, 3, 4, 5
2 桁のものは 6^2 個 ただし，00, 01, ..., 55 であることに注意
3 桁のものは 6^3 個 ただし，000, 001, ..., 555 であることに注意
4 桁のものは 6^4 個 以下同様
5 桁のものは 6^5 個
6 桁のものは 6^6 個
7 桁のものは 6^7 個
8 桁のものは 6^8 個
ここまでで，累計 2015538 個の数（項）である（10 進数で）
目的の 10^7 番目は，9 桁の 6 進数 000000000 を 1 番目として，10^7-2015538 = 7984462 番目の数，数としては 007984461（10 進数）
7984462 の 6 進表示（9 桁）は 443045034
これを逆変換すると 775178157 になる。これが 10^7 番目の項
ということで，プログラム化
func = function(n) { # n 番目の数を求める
k = sum(n > cumsum(6^(1:8))) # k+1 桁の 6 進数の中にある
m = sum(6^(1:k)) # k 桁までの 6 進数の個数の総数
i = n - m - 1 # 求める数の k+1 桁の 6 進表記
num = numeric(k+1)
for (j in 1:(k+1)) {
num[j] = i %% 6
i = i %/% 6
if (i == 0) break
}
as.numeric(paste(c(1, 2, 4, 5, 7, 8)[rev(num)+1], collapse=""))
}
func(10^7) # 775178155
func(10^6) # 44244455
func(10^5) # 1858755
func(10^4) # 115155
func(10^3) # 5455
func(10^2) # 255
func(3333) # 24244
時間が掛かるが，真っ正直に計算して，正しいことを確かめるプログラム
func0 = function(n) {
count = 0
for (i in 1:1e7) {
a = unlist(strsplit(as.character(i), ""))
if (all(! (a %in% c(0, 3, 6, 9)))) {
# print(i)
count = count+1
if (count == n) break
}
}
print(i)
}
func0(3333) # 24244
func0(10^2) # 255
func0(10^3) # 5455
func0(10^4) # 115155
func0(10^5) # 1858755
library(gmp)
func = function(n) {
a = factorize(n) # 素因数分解
u = unique(a) # ユニークな素因数
b = numeric(length(u)) # 同じ素因数の個数
for (i in 1:length(a)) {
suf = which(as.character(a[i]) == u)
b[suf] = b[suf]+1
}
return(prod(b+1)) # 「各素因数の個数+1」の積
}
func(20)
func(472)
func(1073741824)
func("3572947927495273") # 大きな数は文字列で
library(e1071)
d = permutations(7)
a = d[,1] ^((d[,2]/10+d[,3]/10)^(-d[,4]/10)-d[,5]^(-d[,6]-d[,7]/10))
p = which.min(abs(a-exp(1)))
d[a==a[p],]
a[p]
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 2 1 3 4 5 7 6
[2,] 2 3 1 4 5 7 6
> プログラムを書いて、19670808という10進数で表されている数を2進数に変換して登場する最後の数字を調べてくれ！
偶数なんだからさぁ，0 に決まっているじゃん。
って，「登場する最後の数字」って，末尾の数字のことなんだろ？
PVアクセスランキング にほんブログ村
