裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

最寄りの素数

2016年02月26日 | ブログラミング

締め切りが 2月26日(金)10:00 AM なので,
2月26日(金)10:01 AM に公開されるように,自動投稿

それにしても,締め切りまで長いわ~~

最寄りの素数
【概要】
下図のように、1以上のすべての整数をならべます。
入力で指定された数に最も近い素数を求めてください。

【入力と出力】
入力は
14
のような形で、標準入力から来ます。
普通に10進数です。素数が来ることはありません。
出力は、入力の最寄りの素数です。
最寄りの素数が複数ある場合は、小さい順にコンマ区切りで出力してください。
入力が14のときの出力は
7,13,23
になります。


ブルートフォースで何とかなる。

f = function(n=1000000, range=7) {
   k = floor(sqrt(n)-1e-10)+1
   mx <- (k+range)^2
   tbl <- 1:mx
   tbl[1] <- 0
   for (i in 2:floor(sqrt(mx))) {
      if (tbl[i]) {
         mx2 <- mx %/% i
         tbl[2:mx2*i] <- 0
      }
   }
   mid = k^2-k+1
   nx = ifelse(n > mid, k^2-n+1, k)
   ny = ifelse(n > mid, k, n-(k^2-2*k+2)+1)
   Min = .Machine$integer.max
   M = D2 = NULL
   for (k2 in k+(-range:range)) {
      mid = k2^2-k2+1
      for (m in max(1, k2^2-2*k2+2):max(1, k2^2)) {
         if (tbl[m]) {
            mx = ifelse(m > mid, k2^2-m+1, k2)
            my = ifelse(m > mid, k2, m-(k2^2-2*k2+2)+1)
            d2 = (nx-mx)^2+(ny-my)^2
            if (Min >= d2) {
               Min = d2
               M = c(M, m)
               D2 = c(D2, d2)
            }
         }   
      }
   }
   cat(M[D2 == Min], sep=",")
}
con = file("stdin", "r"); f(as.integer(readLines(con, 1)))

コメント

合計が 0 になるのは何通り?

2016年02月18日 | ブログラミング

当初の締め切りが 2/18 10:00AM なので,その 1 分後に投稿されるように予約する

自然数 n に対して、次の等式を考えます。
    □1□2□3□4…□n = 0
四角の部分には、プラス(+)またはマイナス(-)の記号が入ります。
等式を成立させるような記号の入れ方を考えましょう。
例えば、n = 7 のとき、次の 8 通りの入れ方が可能です。

    +1+2-3+4-5-6+7 = 0
    +1+2-3-4+5+6-7 = 0
    +1-2+3+4-5+6-7 = 0
    +1-2-3-4-5+6+7 = 0
    -1+2+3+4+5-6-7 = 0
    -1+2-3-4+5-6+7 = 0
    -1-2+3+4-5-6+7 = 0
    -1-2+3-4+5+6-7 = 0

自然数 n に対し、等式を成立させる記号の入れ方の数を F(n) と定義します。例えば F(7) = 8 です。
同様に、F(3) = 2、F(8) = 14 となることが確かめられます。
標準入力から、自然数 n (1 ≦ n ≦ 50)が与えられます。
標準出力に F(n) の値を出力するプログラムを書いてください。

f = function(n, s) { # f(n, 0) で求まる
    if (n == 0) s == 0
    else if (abs(s) > n*(n+1)/2) 0
    else Recall(n-1, s-n) + Recall(n-1, s+n)
}

> f(3,0)
[1] 2
> f(7,0)
[1] 8
> f(8,0)
[1] 14

しかし,これでは,n = 50 なんてことになると,いつまで待っても答が返ってこない

そこで,以下のプログラムにする
瞬殺である

g = function(n) {
    d = 1
    for (m in 1:(n+1)) {
        len = length(d)
        i = ceiling(len/2)
        a = ifelse(i > len, 0, d[i])
        e = c(rep(0, 2*m), d)
        d = e+rev(e)
    }
    a
}

> options(digits=16)
> g(3)
[1] 2
> g(7)
[1] 8
> g(8)
[1] 14
> g(48)
[1] 1140994231458
> g(49)
[1] 0

gawk で --bignum を指定して

$ cat a.awk
function f(n) {
    len = 1
    d[len] = 1
    for (m = 1; n+1 >= m; m++) {
        i = int((len+1)/2)
        a = i > len ? 0 : d[i]
        for (j = 1; 2*m >= j; j++) {
            e[j] = 0
        }
        for (j = 1; len >= j; j++) {
            e[j+2*m] = d[j]
        }
        len += 2*m
        for (j = 1; len >= j; j++) {
            d[j] = e[j]+e[len+1-j]
        }
    }
    print a

}
BEGIN{
    f(48)
    f(100)
    f(101)
    f(102)
    f(103)
}

$ gawk --bignum --file=a.awk
1140994231458
1731024005948725016633786324
0
0
13252206944539668255309614628

コメント

成長と分裂

2016年02月18日 | ブログラミング

締め切りが 2/18 10:00 AM ということなので,1 分後に投稿されるように予約

バイバイマンのサイズを整数値で表します。
バイバイマンは 1 日毎に体の大きさが倍増します。
一番最初の大きさを「1」とし、1 日経つと 2 倍の「2」、さらに 1 日後は「4」というように、1 日毎に 2 倍になります。
また、一定の大きさをこえると分裂します。
バイバイマンは「1」→「2」→「4」→「8」→「16」と、大きさが 10 を超えたところで分裂します。
分裂後のサイズは、「16」なら「1」と「6」のように 10 の位の数と 1 の位の数になります。
分裂したバイバイマンは、再び大きくなります。
このようなルールで増えるバイバイマンの数を、1 日ごとに調べて、100 日目までを出力してください。
1 日目から 20 日目までは,1, 1, 1, 1, 2, 3, 3, 3, 5, 8, 9, 9, 13, 21, 26, 27, 35, 55, 73, 80 のようになります。

ブルートフォースでやると,とんでもない時間がかかる。
問題をよく調べて見ると,サイズが 1, 2, 4, 8, 6 のものがいくつあるかを数えておいて,次の段階でどのサイズがいくつになるかを整理するとよいことがわかる。
ということで,以下のような簡単なプログラムを得ることができる。
a[1]~a[5] は今の時点でサイズが 1, 2, 4, 8, 6 のものがいくつあるかを表す。次の段階ではそれぞれが 2 倍の大きさの場所に数えられる。ただし,前の段階でサイズ 8 のものは,次の段階では サイズ 1, 6 のものとして,前の段階でサイズ 6 のものは,次の段階では サイズ 1, 2 のものとして数えられる。

a = b = integer(5)
a[1] = 1
for (i in 1:100) {
    cat(sum(a), "\n", sep="")
    b[1] = a[4] + a[5]
    b[2] = a[1] + a[5]
    b[3:5] = a[2:4]
    a = b
}

コメント

共通の友人数

2016年02月17日 | ブログラミング

締め切りはずいぶん先で(鬼が笑うが)2月17日(水)AM10:00
なので,この解答は 2月17日(水)AM10:01 の予約投稿とする。


Facebookなどで表示される「共通の友達」。
あなたの手元には、誰と誰が友達であるかを記したリストがあります。
このリストを元に、「共通の友達」の最大数を求めてください。

例えば、次のような友人関係の場合、共通の友達の最大数は3人になります。
1行目がこのリストに登場する友人番号の最大値、
2行目以降に友人関係がある場合に、そのペアを空白区切りでセットされています。
5
1 2
1 3
1 4
3 4
2 3
2 5
4 5

難易度が「星 4 つのうちの星 3 つ」となっているが,そんなに難しいか?

個人ごとに友達リストを作って,2 人ずつの組合せで共通の友人数を数え,その最大値を出力する(出題文まんまじゃないか)

func = function(fn) {
   con = file(fn, "r")
   n = as.integer(readLines(con, 1))
   f = vector("list", n)
   pair = readLines(con)
   for (i in seq_along(pair)) {
      p = as.integer(unlist(strsplit(pair[i], " ")))
      p1 = p[1]
      p2 = p[2]
      f[[p1]] = append(f[[p1]], p2)
      f[[p2]] = append(f[[p2]], p1)
   }
   n.max = 0
   for (i in 1:(n-1)) {
      for (j in (i+1):n) {
         n.max = max(n.max, length(intersect(f[[i]], f[[j]])))
      }
   }
   cat(n.max)
}
func("test-data1.R")
func("test-data2.R")
func("test-data3.R")

コメント

R の SVG ロゴを base で描いてみる

2016年02月14日 | ブログラミング

オリジナルの奴は,ハロにグラデーションがあるけど...

Quartz (PDF)

r_wkt_gist_file <- "https://gist.githubusercontent.com/hrbrmstr/07d0ccf14c2ff109f55a/raw/db274a39b8f024468f8550d7aeaabb83c576f7ef/rlogo.wkt"
if (!file.exists("rlogo.wkt")) download.file(r_wkt_gist_file, "rlogo.wkt")
s = unlist(strsplit(readLines("rlogo.wkt", warn=FALSE), "[()]"))
s = lapply(s[c(4, 6, 10, 12)], function(s) as.data.frame(matrix(as.numeric(unlist(strsplit(s, "[ ,]"))), ncol=2, byrow=TRUE)))
rng = Reduce(range, s)
old = par(mai=rep(0,4))
plot(c(0, rng[2]), c(rng[1], 0), type="n", bty="n", asp=1, xaxt="n", yaxt="n", xlab="", ylab="")
junk = mapply(function(s, col) polygon(s, col=col, border=NA), s, c("#b8babf", "#ffffff", "#1e63b5", "#ffffff"))
par(old)

コメント