裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

「奇数和平方数」問題

2018年03月22日 | ブログラミング

「奇数和平方数」問題
http://riverplus.hatenablog.com/entry/2018/02/25/143351

自然数 n に対し、連続する n 以下の奇数の和が平方数(自然数の 2 乗で表される数)となるものを探します。
例えば n = 10 の場合、以下の 7 通りです。
項が 1 つだけのものもカウントしている点に注意して下さい。
  1 = 1^2

  1+3 = 2^2

  1+3+5 = 3^2

  9 = 3^2

  1+3+5+7 = 4^2

  7+9 = 4^2

  1+3+5+7+9 = 5^2

このような表し方の数を F(n)と定義します。例えば F(10) = 7 です。
同様に、F(20)=14, F(30)=23, F(1000)=1272 となることが確かめられます。
F(10^7) の値を求めて下さい。
=================================================

ピタゴラス数を数えることになる。

euclid = function(m, n)    {
  repeat {
    temp = n %% m
    if (temp == 0) {
      break
    }
    n = m
    m = temp
  }
  return(m)
}

f = function(N) {
  N = ceiling(N / 2)
  count = 0
  for (m in 2:floor(sqrt(N - 1))) {
    for (n in 1:m) {
      if (m %% 2 != n %% 2 && euclid(m, n) == 1) {
        c = m ^ 2 + n ^ 2
        if (c > N) {
          break
        }
        count = count + N %/% c
      }
    }
  }
  cat(N + count * 2)
}

f(10000000) # 27368012, 3.073s

=================================================

Python だと,若干速いが,アルゴリズムが同じなので本質的な差はない。

import scipy as sp
import math

def f(N):
  N = int(sp.ceil(N / 2))
  count = 0
  for m in range(1, int(sp.floor(sp.sqrt(N - 1)))):
    for n in range(m + 1, N):
      if m % 2 == n % 2:
        continue
      if math.gcd(m, n) == 1:
        c = m ** 2 + n ** 2
        if c > N:
          break
        count += N // c
  return N + 2 * count

import time
s = time.time()
print(f(10000000)) # 27368012, 1.421s
print(time.time() - s)

解説・別解
http://riverplus.hatenablog.com/entry/2018/03/19/220705

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

「メダルツアー」問題

2018年03月21日 | ブログラミング

「メダルツアー」問題
http://riverplus.hatenablog.com/entry/2018/02/25/143351

自然数 w, h に対し、横と縦の長さがそれぞれ w, h となる格子状の道を考えます。
さらに、自然数 m に対し、左から m 番目の縦の道の、各交差点の中間にメダルが 1 個ずつ置いてあります。
例えば下図は、w, h, m = (4, 3, 2) のときの様子です。


この図の左下の点からスタートして、右上のゴールまで最短距離で向かいます。
途中で通過したメダルはすべて拾います。
取り得る経路を等しい確率で選ぶとき、拾うメダルの個数の期待値を F(w, h, m) と定義します。
例えば F(3, 2, 2) = 0.5 です。
取り得る経路は下図の 10 通りで、そのうちメダルを 2 個拾うものが 1 通り、1 個拾うものが 3 通り、0 個拾うものが 6 通りです。
したがって、期待値は 2×(1/10)+1×(3/10)+0×(6/10) = 0.5 と計算できます。


同様に、F(1, 1, 1)=0.5, F(2, 1, 3)=1/3, F(4, 3, 2)=0.6 となることが確かめられます。
小数点以下が 7 桁になるよう四捨五入した F(10, 10, 3) の値を求めて下さい。

=================================================================

末尾に示す簡単なプログラムで w, h の小さい部分についての解を求めると結果は以下のようになる
この問題において,m は解に無関係であることはすぐに分かる
関係性は,ずばり,h/(w+1)


f = function(w, h,  m) {
  g = function(k) {
    x = y = integer(wh+2)
    x[1] = y[1] = 1
    p = integer(wh)
    p[k] = 1
    for (i in 1:wh) {
      if (p[i] == 1) {
        x[i+1] = x[i]+1
        y[i+1] = y[i]
      } else {
        x[i+1] = x[i]
        y[i+1] = y[i]+1
      }
    }
    count = 0
    for (i in 1:wh1) {
      if (x[i] == m && x[i+1] == m && y[i+1]-y[i] == 1) {
        count = count+1
      }
    }
    return(count)
  }
  wh = w+h
  wh1 = wh+1
  a = combn(wh, w, g)
  b = table(a)
  l = length(b)
  value = as.integer(names(b))
  e  = 0
  for (i in 1:l) {
    e = e+value[i]*b[i]
  }
  return(e/sum(b))
}

別解,解説

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ナンプレを解くプログラム

2018年03月07日 | ブログラミング

ナンプレを解くプログラム

締め切りが 2018/03/07 10:00 AM なので,その 1 分後に投稿されるように予約

【概要】
簡単なナンプレの問題を解くプログラムを書いてみましょう。

【問題】
次のようなナンプレがあります。
・2x2のマス(Box)が縦横ふたつずつ並んだ4x4のマスがある
・それぞれのマスには1~4の数字が入り、それらは行/列の4つのマスに同じ数字がくることはない(行/列の合計は10となる)
・2x2のマス(Box)にも同じ数字がくることはなく、合計は10となる
・最初に1~4の数字がランダムにひとつずつ設定されている(同じ行/列/Boxに複数の数字が設定されることは無い)

初期状態の例(_は数字の入っていない状態を表します)



このナンプレを解くプログラムを書いて、数字を埋めてください。

【入出力】
CSV形式のテストデータを標準入力で受け取る想定で、結果を標準出力するプログラムを作ってください。

上図のナンプレでは、以下のようなCSVファイルが入力されます。

    1.    1,0,0,0
    2.    0,0,3,0
    3.    0,0,0,4
    4.    0,2,0,0

これに対して、以下のようなデータを標準出力から出力してください。

    1.    1,3,4,2
    2.    2,4,3,1
    3.    3,1,2,4
    4.    4,2,1,3
【入出力サンプル】
Input
1,0,0,0
0,0,3,0
0,0,0,4
0,2,0,0

Output
1,3,4,2
2,4,3,1
3,1,2,4
4,2,1,3

================================

ブルートフォースで簡単に書く

permutations = function (n)
{
  z = matrix(1)
  for (i in 2:n) {
    x = cbind(z, i)
    a = c(1:i, 1:(i - 1))
    z = matrix(0, ncol = ncol(x), nrow = i * nrow(x))
    z[1:nrow(x), ] = x
    for (j in 2:i - 1) {
      z[j * nrow(x) + 1:nrow(x), ] = x[, a[1:i + j]]
    }
  }
  z
}

f = function(x) {
  g = function(a) {
    b = which(a != 0, arr.ind = TRUE)
    x = b[1] + 2 * b[2] - 2
    y = a[b]
    A = permutations(4)
    A = A[A[, x] == y, ]
    return(A)
  }
  x = matrix(x, ncol = 4, byrow = TRUE)
  A = g(x[1:2, 1:2])
  B = g(x[1:2, 3:4])
  C = g(x[3:4, 1:2])
  D = g(x[3:4, 3:4])
  for (i in 1:6) {
    for (j in 1:6) {
      for (k in 1:6) {
        for (m in 1:6) {
          # 横方向のチェック
          if (length(unique(c(A[i, 1], A[i, 3], B[j, 1], B[j, 3]))) != 4) next
          if (length(unique(c(A[i, 2], A[i, 4], B[j, 2], B[j, 4]))) != 4) next
          if (length(unique(c(C[k, 1], C[k, 3], D[m, 1], D[m, 3]))) != 4) next
          if (length(unique(c(C[k, 2], C[k, 4], D[m, 2], D[m, 4]))) != 4) next
          # 縦方向のチェック
          if (length(unique(c(A[i, 1], A[i, 2], C[k, 1], C[k, 2]))) != 4) next
          if (length(unique(c(A[i, 3], A[i, 4], C[k, 3], C[k, 4]))) != 4) next
          if (length(unique(c(B[j, 1], B[j, 2], D[m, 1], D[m, 2]))) != 4) next
          if (length(unique(c(B[j, 3], B[j, 4], D[m, 3], D[m, 4]))) != 4) next
          cat(sprintf("%d,%d,%d,%d\n", A[i, 1], A[i, 3], B[j, 1], B[j, 3]))
          cat(sprintf("%d,%d,%d,%d\n", A[i, 2], A[i, 4], B[j, 2], B[j, 4]))
          cat(sprintf("%d,%d,%d,%d\n", C[k, 1], C[k, 3], D[m, 1], D[m, 3]))
          cat(sprintf("%d,%d,%d,%d\n", C[k, 2], C[k, 4], D[m, 2], D[m, 4]))
        }
      }
    }
  }
}

# f(c(1, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 4, 0, 2, 0, 0))
x = scan(file("stdin", "r"), sep=",")
f(x)

入力
1,0,0,0
0,0,3,0
0,0,0,4
0,2,0,0

出力
1,3,4,2
2,4,3,1
3,1,2,4
4,2,1,3

入力
0,0,0,4
0,2,0,0
0,0,3,0
1,0,0,0

出力
3,1,2,4
4,2,1,3
2,4,3,1
1,3,4,2

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村