探索を逆方向から攻める

整数の平方根をとったとき，小数点以下 1 ～ 6 桁が全て同じ数字になるような最小の整数を求めよ。

求める数 x の平方根の整数部分を a 小数部分を b とする。つまりその整数は (a + b) ^ 2 = a ^ 2 + 2 * a * b + b ^2 であるとする。
小数部分を小数点以下 6 桁まで考えると上の式で求められる数値は，整数 x より小さな実数になる。差がきわめて近い物を探索する。

b = 0.111111 * 1:9 # b は 0.111111 ～ 0.999999 の 9 種類
b2 = 2 * b
b.sq = b ^ 2
n = 100000 # 探索範囲上限
for (a in 1:n) {
  x = a * b2 + b.sq # 求める整数から a^2 を引いたものの近似値
  if (any(abs(round(x)-x) < 1e-6)) { # round(x) が実際の求めるべき整数から a^2 を引いたもの
    index = which.min(abs(round(x)-x))
    cat("(", a, "+", b[index], ") ^ 2 = ", (a+b[index])^2, "\n")
  }
}

( 55557 + 0.111111 ) ^ 2 =  3086592595
( 83333 + 0.666666 ) ^ 2 =  6944500000
>
> 55557.111111 ^ 2
[1] 3086592595

> sqrt(3086592595)
[1] 55557.111111 確かに小数点以下 6 桁が同じ数字

求める整数は，3086592595

計算所要時間は n = 100000 まで探索して（2 つの解を求めて） 0.3 秒

素数の分布

素数の存在範囲（上限）n と階級幅 w を与えて，素数の度数分布を求める。（たいしておもしろくない）

> func = function(n, w) {
+   x = 2:n
+   f = table(ceiling(x[1-outer(x, x)] / w))
+   m = nchar(as.character(n))
+   label = as.integer(names(f)) * w
+   label = sprintf("%0*i-%0*i:", m, label-w+1, m, label)
+   for (i in seq_along(f)) {
+     cat(label[i], rep("*", f[i]), "\n", sep="")
+   }
+ }

> func(30, 5)
01-05:***
06-10:*
11-15:**
16-20:**
21-25:*
26-30:*

> func(100, 15)
001-015:******
016-030:****
031-045:****
046-060:***
061-075:****
076-090:***
091-105:*

三次元の合成数？

辺の長さが整数で、体積がnの直方体を考える。
縦、横、高さのいずれも長さが1以外の直方体を作ることができる整数nが
1≦n≦300の範囲にいくつあるかを求めよ。

いわれてみれば確かに，「三次元の合成数」といえなくもないが。
要するに，約数が 3 個以上あるのが該当する整数だ。

factorize = function(n) {
  f = NULL
  repeat {
    prime = TRUE
    if (n == 1) return(f)
    for (i in 2:max(sqrt(n), 2)) {
      if (n %% i == 0) {
        f = c(f, i)
        n = n %/% i
        prime = FALSE
        break
      }
    }
    if (prime) {
      return(c(f, n))
    }
  }
}
count = 0
for (i in 2:300) {
  k = factorize(i)
  if (length(k) >= 3) {
    cat(i, ": ", k, "\n")
    count = count+1
  }
}
cat("Ans. =", count)


> cat("Ans. =", count)
Ans. = 143

解禁日: 2015/03/16 10:00

逆ポーランド記法の電卓

いろいろ拡張した。

stack = NULL
push = function(x) {
  c(x, stack) ->> stack
}
pop = function() {
  x = stack[1]
  stack[-1] ->> stack
  return(x)
}
debug = TRUE # スタックの状態を表示
calc = function(eq) {
  cat(eq, "\n")
  func = c("sqrt", "log", "log10", "exp", "sin", "cos", "tan", "abs")
  eq = unlist(strsplit(eq, " "))
  for (i in eq) {
    if (grepl("^[0-9]+$", i)) {
      push(as.numeric(i))
    } else if (grepl("[-+*/]", i)) { # 四則演算
      eval(parse(text=sprintf("push(pop() %s pop())", i)))
    } else if ("^" == i) { # べき乗
      x = pop()
      eval(parse(text=sprintf("push(pop() %s x)", i)))
    } else if (any(func == i)) { # 1 引数の関数
      eval(parse(text=sprintf("push(%s(pop()))", func[which(i == func)])))
    }
    if (debug) {
      rs = rev(stack) # 「右」がスタックの「上」ということで
      cat("[", sprintf("%f", rs[1]), sprintf(", %f", rs[-1]), "]\n", sep="")
    }
  }
  cat(pop(), "\n")
  if (debug) {
    cat("[", stack, "]", sep="")
  }
}
calc("3 4 + 5 *")
calc("2 sqrt 1 + 2 sqrt 1 - *")
calc("2 3 ^")
calc("10 log")
calc("100 log10")
calc("2 exp")
calc("111111 2 ^")

解禁日：2015/04/13 10:00

制約条件のある部分和問題

以下に示す 30 個の整数を 10 個ずつに分け，それぞれの数値の和が 3 つとも等しくなるような分割のしかたを求めよ。

138, 104, 139, 124, 157, 128, 134, 150, 157, 81,
159, 159, 98, 134, 147, 154, 113, 117, 138, 156,
116, 147, 81, 150, 95, 131, 156, 91, 105, 132

以下のプログラムで，簡単に求めることができるが，解は無数にある。

func = function(i, Sum) {
    if (i == n+1) {
        invisible(Sum == Result)
    }
    else if (func(i+1, Sum)) {
        invisible(TRUE)
    }
    else if (func(i+1, Sum+x[i])) {
        TRUE ->> w[i]
        invisible(TRUE)
    } else {
        invisible(FALSE)
    }
}

x0 = c(138, 104, 139, 124, 157, 128, 134, 150, 157, 81, 159, 159,
98, 134, 147, 154, 113, 117, 138, 156, 116, 147, 81, 150, 95,
131, 156, 91, 105, 132)
Result = sum(x0)/3

repeat {
    x = sample(x0)
    n = length(x0)
    w = logical(n)

    func(1,0)
    a1 = sort(x[w]) # 解1
    if (length(a1) != 10) next
    x = x[!w] # 残りの 20 個
    n = length(x)
    w = logical(n)

    func(1,0)
    a2 = sort(x[w]) # 解2
    a3 = sort(x[!w]) # 解3
    if (length(a1) == length(a2) && length(a1) == length(a3)) {
        cat(length(a1), sum(a1)); print(a1)
        cat(length(a2), sum(a2)); print(a2)
        cat(length(a3), sum(a3)); print(a3)
        break
    }
}

10 1297 [1]  98 104 116 124 128 132 138 139 159 159
10 1297 [1]  91  95 105 117 131 147 147 150 157 157
10 1297 [1]  81  81 113 134 134 138 150 154 156 156

とか

10 1297 [1]  81  98 113 116 132 134 154 156 156 157
10 1297 [1]  81 104 124 128 131 134 138 139 159 159
10 1297 [1]  91  95 105 117 138 147 147 150 150 157

とか

解禁日: 2015/03/16 10:00

1!+2!+3!+…+2015! の下三桁

1!+2!+3!+…+2015!の下三桁を求めてください。

a = numeric(2015)
f = 1
for (i in seq_along(a)) {
  a[i] = f = (f * i) %% 1000
}
sum(a) %% 1000

答え：313

2015/03/07 10:00 解禁

双子素数

差が 2 である 2 つの素数の組

x = 2:1000
x = x[1-outer(x, x)]
s = diff(x) == 2
cat(sub(", $", "", paste(sprintf("(%s), ", apply(cbind(x[s], (x[-1])[s]), 1, paste, collapse=", ")), collapse="")))

モールス符号（アルファベット）

カナ文字については既出

出題は，まずはモールス符号を平文化すればわかる。

m = c(
"01", "1000", "1010", "100", "0", "0010", "110", "0000", "00", "0111",
"101", "0100", "11", "10", "111", "0110", "1101", "010", "000", "1",
"001", "0001", "011", "1001", "1011", "1100",
"010101", "110011", "001100", "101011", "100001", "10010", "011010", "10110", "101101",
"01111", "00111", "00011", "00001", "00000",
"10000", "11000", "11100", "11110", "11111")
# カタカナ，数字
e = c(
letters, ".", ",", "?", "!", "-", "/", "@", "(", ")",
"1", "2", "3", "4", "5", "6", "7", "8", "9", "0")
# アルファベット，数字　→　モールス符号
em = function(str) {
    m[sapply(unlist(strsplit(str, "")), grep, k)]
}
# モールス符号　→　アルファベット，数字
me = function(str) {
    paste(k[sapply(str, function(s) grep(sprintf("^%s$", s), m))], collapse="")
}

me(c("0100", "1010", "11", "10010", "00111", "00000", "10000", "11100", "10010", "11000", "11000", "10000", "00000"))

答は 19940520 だな

三角形ができる 3 つの数の組み合わせ数

入力 n に対して 1≦a≦b≦c≦n（a, b, c, n はすべて整数）を満たす a, b, c の中で，これらを 3 辺とする三角形が成り立つ場合の組み合わせを求めるプログラムを書け。

n = 3 のときは，7 である。さて，n = 1000 のときの答えを求めよ。

R で素直に書くと，とても時間が掛かる。

> system.time({
+     a=0
+ for (i in 1:1000)
+ for (j in i:1000)
+ for (k in j:1000) {
+     a=a+(2*max(i,j,k) < i+j+k)
+ }
+ print(a)
+ })
[1] 83708750
   ユーザ   システム       経過  
   399.066      1.916    397.616

でも，GAWK で書くと，もっと時間が掛かる。

BEGIN {
for (i = 1; i <= 1000; i++)
for (j = i; j <= 1000; j++)
for (k = j; k <= 1000; k++) {
mx=i > j ? i : j
a+=2*(mx > k ? mx : k) < i+j+k
}
print a
}

このプログラムを収めたファイルを a.awk として，

> date;gawk -f a.awk;date
2015年 2月15日 日曜日 21時43分39秒 JST
83708750
2015年 2月15日 日曜日 21時52分13秒 JST

514 秒掛かった。

R も GAWK もインタープリタだけど，R の方が優れているのか？？

解禁日 2015/03/02 10:00