裏 RjpWiki

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

整数問題 3 題(その2)

2015年03月30日 | ブログラミング

久しぶりに AWK(GAWK)で書いてみた。

やはり R で書くとすっきり書ける。

$ cat prog.awk
#!/usr/local/bin/gawk -f
function F(term,    a, b, c, n, count, d) {
  a = 3
  b = 0
  c = 2
  count = 1
  for (n = 3; ; n++) {
    d = a+b
    if (d % n == 0 && ++count == term) return d
    a = b
    b = c
    c = d
  }
}

function divisor(n,     i) {
  if (n % 2 == 0) return 2
  else if (n % 3 == 0) return 3
  maxitr = int(sqrt(n))
  for (i = 5; ; i += 4) {
    if (i > maxitr) return n
    else if (n % i == 0) return i
    i = i + 2
    if (i > maxitr) return n
    else if (n % i == 0) return i
  }
}

function G(n,     d) {
  for (;;) {
    d = divisor(n)
    n = n/d
    if (n == 1) return d
  }
}
 
function H(n,     i, tbl, j, prime, sum) {
  for (i = 1; i <= n; i++) {
    tbl[i] = i
  }
  tbl[1] = 0
  for (i = 2; i <= int(sqrt(n)); i++) {
    if (tbl[i]) {
      for (j = 2; j <= int(n / i); j++) {
        tbl[j*i] = 0
      }
    }
  }
  for (i = 1; i <= n; i++) {
    if (tbl[i] > 0) {
      prime[++count] = tbl[i]
    }
  }
  for (i = 1; i <= count; i++) {
    if (prime[i] <= count) {
      sum += prime[prime[i]]
    }
  }
  return sum
}

BEGIN {
  P = F(30)
  Q = G(P)
  R = H(Q)
  printf "P = %i, Q = %i, R = %i\n", P, Q, R
}

$ time ./prog.awk
P = 63088012960224, Q = 120739, R = 72834276
0.132u 0.005s 0:00.13 100.0%    0+0k 0+0io 0pf+0w

コメント

模範解答としてはどうかなあ

2015年03月30日 | ブログラミング

素因数分解なんだけど

> 対象の数を、まずは2でこれ以上割り切れなくなるまで割り尽くします。
> 次にその商を3で割り尽くします。
> さらにその商を4,5,6,・・・の順で割りつくしていけば...

4 や 6... で割ることを試みるのは,無駄だと知らない訳か?

5 以上の素数は 6n±1 (n = 1, 2, ...)で表すことができる。しかし,6n±1 が必ず素数であるとは限らない(25 や 49)。

さらには,模範解答には書いていないのだけど,調べる約数はどこまでにするか。
模範解答氏は,「対象の数-1」まで調べてしまいそう(藁)。

コメント

整数問題 3 題

2015年03月27日 | ブログラミング

Q1. 数列 F(n) を以下の漸化式で定義する。
F(0) = 3, F(1) = 0, F(2) = 2,
F(n) = F(n-2) + F(n-3) (n>2のとき)
このような数列において,n > 1 かつ F(n) の値が n が割り切れる場合を考える。
たとえば,n = 2, 3, 5, 7, 11, 13, 17, 19 のとき,F(n) の値はそれぞれ n で割り切れる。
19 は,このような性質を持つ 8 番目に小さな n で,この n に対する F(n) の値は 209 である。
このような性質を持つ 30 番目に小さな n に対する F(n) の値 P を求めよ。

Q2. 整数 n に対し,n の素因数のうちで最大のものを G(n) とする。たとえば,G(123456) = 643である。
G(P) の値 Q を求めよ。

Q3. 2, 3, 5, 7 はそれぞれ 1 番目,2 番目,3 番目,4 番目の素数である。k 番目の素数 p に対し,k もまた素数である場合に,p を「素数番目の素数」と呼ぶことにする。たとえば,3 と 5 は素数番目の素数であるが、2 と 7 はそうではない。
整数 n に対し,n 以下の素数番目の素数の和を H(n) と定義する。たとえば,H(100) = 317, H(1000) = 15489 である。
H(Q) の値 R を求めよ。

system.time({
  F = function(term) {
    a = 3
    b = 0
    c = 2
    n = 2
    count = 1
    repeat {
      n = n + 1
      d = a + b
      if (d %% n == 0 && (count = count + 1) == term) return(d)
      a = b
      b = c
      c = d
    }
  }
  P = F(30)

  G = function(n) {
    div = function(n) {
      if (n%%2 == 0) return(2)
      else if (n%%3 == 0) return(3)
      maxitr = floor(sqrt(n))
      i = 1
      repeat {
        i = i + 4
        if (i > maxitr) return(n)
        else if (n%%i == 0) return(i)
        i = i + 2
        if (i > maxitr) return(n)
        else if (n%%i == 0) return(i)
      }
    }
    result = NULL
    repeat {
      d = div(n)
      result = append(result, d)
      n = n/d
      if (n == 1) break
    }
    return(max(result))
  }
  Q = G(P)

  H = function(n) {
    tbl = 1:n
    tbl[1] = 0
    for (i in 2:floor(sqrt(n))) {
      if (tbl[i]) tbl[2:(n%/%i) * i] = 0
    }
    prime = tbl[tbl > 0]
    sum(prime[prime[!is.na(prime[prime])]])
  }
  R = H(Q)

  cat(sprintf("P = %.f, Q = %.f, R = %.f\n", P, Q, R))
})

実行例

P = 63088012960224, Q = 120739, R = 72834276
   ユーザ   システム       経過  
     0.117      0.000      0.126

コメント

分数の計算(和)

2015年03月24日 | ブログラミング

文字列 "a/b" の形式で表される 2 つの分数の和を,同じ形式で表示する
分母が 1 の場合は,分子だけを表示する

func = function(s1, s2) {
  euclid = function(m, n) { # ユークリッドの互除法
    while ((temp = n%%m) != 0) {
      n = m
      m = temp
    }
    return(m)
  }
  a = as.integer(unlist(strsplit(s1, "/")))
  b = as.integer(unlist(strsplit(s2, "/")))
  denominator = a[2] * b[2]
  numerator = a[1] * b[2] + b[1] * a[2]
  gcm = euclid(denominator, numerator)
  denominator = denominator / gcm
  numerator = numerator / gcm
  if (denominator == 1)
    cat(numerator)
  else
    cat(numerator, denominator, sep = "/")
  cat("\n")
}
func("5/6", "1/10") # 14/15
func("1/3", "2/3")  # 1
func("1/3", "2/7")  # 13/21
func("2/8", "3/5")  # 17/20
func("3/10", "1/6") # 7/15
func("3/4", "5/8")  # 11/8
func("2/5", "2/3")  # 16/15

コメント

探索を逆方向から攻める(2)

2015年03月23日 | ブログラミング

探索を逆方向から攻める」は不完全であった。条件を満たす最小の整数を漏らしていた。

b = 0.111111 * 1:9 # b は 0.111111 ~ 0.999999 の 9 種類
b2 = 2 * b
b.sq = b^2
n = 1e+05 # 探索範囲上限
for (a in 1:n) {
  x = a * b2 + b.sq # 求める整数から a^2 を引いたものの近似値
  if (any(abs(round(x) - x) < 1e-05)) { # round(x) が実際の求めるべき整数から a^2 を引いたもの(チェックを緩めた)
    index = which.min(abs(round(x) - x))
    y = sqrt(round((a + b[index])^2)) # 実際に平方根をとって,小数点以下に 6 個の数字が並ぶかをチェックする
    z = unlist(strsplit(as.character(y), ""))
    if (any(grepl("\\.", z))) {
      i = which(z == ".")
      u = z[(i + 1):(i + 6)] # 小数点以下 6 桁を取り出す
      if (all(u[1] == u)) { # 全部同じか
        options(digits=15)
        cat("sqrt(", round((a + b[index])^2), ") =", y, "\n")
        break
      }
    }
  }
}

実行結果

sqrt( 48273160 ) = 6947.88888800044

実行時間は,0.052 秒であった。

コメント

sample 関数

2015年03月19日 | ブログラミング

1 ~ n までの整数乱数を発生する sample 関数を使うとき,第1引数はベクトルを指定しない方が速い。

system.time({
    set.seed(1)
    a = sample(1:1e8, 1e5, replace=TRUE)
})

system.time({
    set.seed(1)
    b = sample(1e8, 1e5, replace=TRUE)
})

all(a == b)

実行結果

> system.time({
+     set.seed(1)
+     a = sample(1:1e8, 1e5, replace=TRUE)
+ })
   ユーザ   システム       経過  
     0.111      0.092      0.201
>
> system.time({
+     set.seed(1)
+     b = sample(1e8, 1e5, replace=TRUE)
+ })
   ユーザ   システム       経過  
     0.001      0.000      0.001
>
> all(a == b)
[1] TRUE

なぜそうなのかは,sample のソースを見ればわかる。

コメント

英語で数を読む

2015年03月16日 | ブログラミング

符号付き 32 ビット整数を英語で読めとのお達しで...

o1 = c("One", "Two", "Three", "Four", "Five", "Six", "Seven", "Eight", "Nine", "Ten", "Eleven",
  "Twelve", "Thirteen", "Fourteen", "Fifteen", "Sixteen", "Seventeen", "Eighteen", "Nineteen")
o2 = c("", "Twenty", "Thirty", "Forty", "Fifty", "Sixty", "Seventy", "Eighty", "Ninety")
o3 = "Hundred"
o4 = c("Thousand", "Million", "Billion")
count = function(n) {
  cat("n =", n, "\n")
  ans = NULL
  if (0 > n) {
    ans = "Negative"
    n = -n
  }
  s = unlist(strsplit(as.character(abs(n)), ""))
  ns = length(s)
  if (12 > ns) s = c(rep(0, 12 - ns), s)
  s = matrix(as.numeric(s), ncol = 3, byrow = TRUE)
  if (all(s == 0)) ans = "Zero"
  for (i in 1:4) {
    if (any(s[i, ] != 0)) {
      if (s[i, 1] > 0) ans = c(ans, o1[s[i, 1]], o3)
      if (s[i, 2] >= 2) {
        ans = c(ans, o2[s[i, 2]])
        if (s[i, 3] != 0) ans = c(ans, o1[s[i, 3]])
      }
      if (1 >= s[i, 2] && (s[i, 2] + s[i, 3] != 0)) ans = c(ans, o1[s[i, 2] * 10 + s[i, 3]])
      if (3 >= i) ans = c(ans, o4[4 - i])
    }
  }
  paste(ans, collapse = " ")
}
# テストデータ
q = c(7, 123, 4567, 89012, 0, -34, -5678901, 1111111111)
ans = c("Seven", "One Hundred Twenty Three", "Four Thousand Five Hundred Sixty Seven", "Eighty Nine Thousand Twelve",
  "Zero", "Negative Thirty Four", "Negative Five Million Six Hundred Seventy Eight Thousand Nine Hundred One",
  "One Billion One Hundred Eleven Million One Hundred Eleven Thousand One Hundred Eleven")
for (i in seq_along(q)) {
  a = count(q[i])
  cat("a: ", a, "\n")
  cat("b: ", ans[i], "\n")
  if (a != ans[i]) cat("Wrong!!\n\n")
}

結果は一応 OK

n = 7
a:  Seven
b:  Seven

n = 123
a:  One Hundred Twenty Three
b:  One Hundred Twenty Three

n = 4567
a:  Four Thousand Five Hundred Sixty Seven
b:  Four Thousand Five Hundred Sixty Seven

n = 89012
a:  Eighty Nine Thousand Twelve
b:  Eighty Nine Thousand Twelve

n = 0
a:  Zero
b:  Zero

n = -34
a:  Negative Thirty Four
b:  Negative Thirty Four

n = -5678901
a:  Negative Five Million Six Hundred Seventy Eight Thousand Nine Hundred One
b:  Negative Five Million Six Hundred Seventy Eight Thousand Nine Hundred One

n = 1111111111
a:  One Billion One Hundred Eleven Million One Hundred Eleven Thousand One Hundred Eleven
b:  One Billion One Hundred Eleven Million One Hundred Eleven Thousand One Hundred Eleven

コメント

カプレカ数(その2)

2015年03月16日 | ブログラミング

カプレカ数」での,二番目の定義によるものを出力するプログラム。

かなり時間のかかるプログラムであり,5 番目のカプレカ数を出力するまでには 35 秒かかる。

N = 631764
for (i in 0:N) {
  str = unlist(strsplit(as.character(i), ""))
  big = as.numeric(paste(sort(str, decreasing = TRUE), collapse=""))
  small = as.numeric(paste(sort(str), collapse=""))
  if (big - small == i) print(i)
}

[1] 0
[1] 495
[1] 6174
[1] 549945
[1] 631764
   ユーザ   システム       経過  
    34.915      0.202     34.773

コメント

カプレカ数

2015年03月16日 | ブログラミング

0 から 1000 までのカプレカ数を求めよとのことだ...

カプレカ数とは馴染みのないものであるが,定義は 2 通りあるそうだ。

1. 2 乗して前の部分と後ろの部分に分けて(偶数桁 2n 桁である場合は先頭 n 桁と末尾 n 桁に分け,奇数桁 2n + 1 桁である場合は先頭 n 桁と末尾 n + 1 桁に分けて)和を取ったとき,元の値に等しくなるもの。

2. 桁を並べ替えて最大にしたものから最小にしたものの差を取ったとき,元の値に等しくなるもの。

今回は,前者のことをいっているようだ。

Wikipedia で調べると
カプレカ数
1, 9, 45, 55, 99, 297, 703, 999, 2223, 2728, 4879, 4950, 5050, 5292, …(オンライン整数列大辞典の数列 A006886 http://oeis.org/A006886)が簡単に見つかる。
プログラムを書いて,0 ~ 5300 の範囲を探索すると,得られる出力と違うところがある。

a1 = 0:5300
a2 = a1^2
old = options(scipen=10)
a3 = strsplit(as.character(a2), "")
options(old)
a4 = sapply(a3, function(x) {
            y = unlist(x)
            n = length(y)
            if (n == 1) {
              as.numeric(paste(y, collapse=""))
            } else {
              m = n %/% 2
              as.numeric(paste(y[1:m], collapse="")) + as.numeric(paste(y[(m+1):n], collapse=""))
            }
          })
a1[a1 == a4]

> a1[a1 == a4]
 [1]    0    1    9   45   55   99  297  703  999 2223 2728 4950 5050

0 は明白な カプレカ数だと思うがそれはさておき,プログラムによる出力には 4879 はない。4879^2 = 2380 と 4641 なので,2380 + 4641 = 7021 になるので,カプレカ数ではないはず。2380 の最後の 0 を無視するというルールならば,238 + 4641 = 4879 になるのでカプレカ数ということになるのだが。同様に,5292^2 = 2800 と 5264 なので,28 + 5264 = 5292 ではあるが?と根拠を求めてネットサーフィンすると,

Karekar numbers
http://www.numbersaplenty.com/set/Kaprekar_number/

に,

Note that the second part can start with zero:  5292^2 = {28} {005264}  and 28+5264=5292.

と書いてあった。後ろ半分の数字の前に 0 があってもよいということは,前の部分の後ろの 0 は無視してもよいということである(無視しない場合も考えよということ)。このルールも加えてプログラムを修正する。

プログラムに for を使うことになったがやむを得ない。

N = 5300
a1 = 0:N
a2 = a1^2
old = options(scipen=10)
a3 = strsplit(as.character(a2), "")
options(old)
a4 = a5 = integer(N+1)
for (i in seq_along(a1)) {
  y = unlist(a3[[i]])
  n = length(y)
  if (n == 1) {
    a4[i] = a5[i] = as.numeric(paste(y, collapse=""))
  } else {
    m = n %/% 2
    a4[i] = as.numeric(paste(y[1:m], collapse="")) + # 前半分の後ろの 0 を無視しない
            as.numeric(paste(y[(m+1):n], collapse=""))
    a5[i] = as.numeric(paste(gsub("0*$", "", y[1:m]), collapse="")) + # 前半分の後ろの 0 を無視する
            as.numeric(paste(y[(m+1):n], collapse=""))
  }
}
a1[a1 == a4 | a1 == a5] # 前半分の後ろの 0 を無視しない場合と無視する場合のどちらか一方でも可ということで

> a1[a1 == a4 | a1 == a5]
 [1]    0    1    9   45   55   99  297  703  999 2223 2728 4879 4950 5050 5292

Wikipedia などの例示と一致した。

10 万までで 5 秒弱かかって,以下の結果を得る。

> print(a1[a1 == a4 | a1 == a5])})
 [1]     0     1     9    45    55    99   297   703   999  2223  2728  4879  4950  5050  5292  7272
[17]  7777  9999 17344 22222 38962 77778 82656 95121 99999

コメント

完全数を探せ

2015年03月11日 | ブログラミング

「10000 以下の最大の完全数を答えよ」ということなのだが,完全数は小さい方から順に,6, 28, 496, 8128, 33550336, 8589869056, 137438691328, 2305843008139952128, 2658455991569831744654692615953842176, 191561942608236107294793378084303638130997321548169216, ... ということは,ちょっと調べるとすぐ分かる。
「それをいっちゃ~,おしめえよ!」なので,以下のプログラム

n = 10000
repeat {
  m = 1:n # その数までの自然数
  m = m[n %% m == 0] # 割りきれるなら約数(その数自身も含む)
  if (sum(m) == 2*n) { # 「その数の約数の和が,その数の2倍」が完全数の定義
    print(n) # 見つかったらプリントしてループから脱出
    break
  }
  n = n-2 # 完全数は,(今のところ)偶数であるという事実を利用
}

0.4 秒ほどで,8128 が表示される。
(ideone では 1.1 秒ほどかかる)

以下のようにすれば,所要時間は半分となる。

n = 10000
repeat {
  m = 1:(n/2) # その数の半分(約数 2 の相棒)までの自然数(28 なら,約数 2 と 14 の関係)
  m = m[n %% m == 0] # 割りきれるなら約数(その数自身を除く)
  if (sum(m) == n) { # 「その数の約数(その数自身を除く)の和が,その数に等しい」が完全数の定義
    print(n) # 見つかったらプリントしてループから脱出
    break
  }
  n = n-2 # 完全数は,(今のところ)偶数であるという事実を利用
}

コメント

特定の性質を持つ整数

2015年03月09日 | ブログラミング

数字の 0 から 9 までを一度ずつ使ってできる数で,
0 を除く全ての一桁の数で割り切れ,
この数に含まれる任意の隣り合う二桁でも割り切れる,
そんな数を探せ。

 library(e1071)
 d = cbind(permutations(9), 0) # 2 と 5 を因数に含むので,末尾が 0 の数だ
 w = 10^(9:0) # 10 進数の重み
 x = rowSums(t(t(d)*w)) # 10 進数として扱う
 f1 = x %% (7*8*9) == 0 # 2 ~ 9 のすべてを因数として持つか?
 d = d[f1,]
 x = x[f1]
 f12 = x %% (d[,1] * 10 + d[,2]) # 隣り合う 2 桁の数値を因数として持つか
 f23 = x %% (d[,2] * 10 + d[,3])
 f34 = x %% (d[,3] * 10 + d[,4])
 f45 = x %% (d[,4] * 10 + d[,5])
 f56 = x %% (d[,5] * 10 + d[,6])
 f67 = x %% (d[,6] * 10 + d[,7])
 f78 = x %% (d[,7] * 10 + d[,8])
 f89 = x %% (d[,8] * 10 + d[,9])
 f9x = x %% (d[,9] * 10)
 f2 = f12+f23+f34+f45+f56+f67+f78+f89+f9x == 0
 x[f2]

0.2 秒(!!)ほどで見つかる。3912657840

上のプログラムは冗長ともいえるが,下のように書き直したプログラムは,短いけどお行儀は悪い(良いのか?)

 library(e1071)
 d = permutations(9) # 2 と 5 を因数に含むので,末尾が 0 の数だ
 w = 10^(9:1) # 10 進数の重み
 x = colSums(t(d)*w) # 10 進数として扱う
 f1 = x %% (7*8*9) == 0 # 2 ~ 9 のすべてを因数として持つか?
 d = d[f1,]
 x = x[f1]
 x[rowSums(x %% (d*10 + cbind(d[,-1], 0))) == 0]


コメント

プログラム書法(2)

2015年03月06日 | ブログラミング

前の記事(プログラム書法)で,似たようなことを行う 2 つの関数を 1 つにまとめたが,今回はさらに,もう 2 つの関数を抱き合わせにする。

関数名を覚えるのが面倒ということは些細な理由である。関数のメンテナンスが「間違いなく」,「容易に」できるというメリットがある。また,短い関数は見通しがよいので,バグが入りにくい(か?)。

前回のプログラムで,wilcox.test, t.test の呼び出し時に paired=FALSE,t.test のときに var.equal=FALSE をわざわざ指定していたが,デフォルトなのでそれらを取り除くと,4 つの検定関数の呼び出し時の引数指定が全く同じなので,なにかうまい方法があるんじゃないかと思う人もいよう。そのうまい方法とは,以下のように,変数(func)に関数オブジェクト(wilcox.test など)を代入することである。

# 全ての変数に wilcox.test / t.test / kruskal.test / oneway.test を行う関数
Rep = function(data, group, method = c("wilcox", "t", "kruskal", "oneway"), sort = TRUE) {
  name = method = match.arg(method)
  if (method == "wilcox") {
    func = wilcox.test
  } else if (method == "t") {
    func = t.test
  } else if (method == "kruskal") {

    func = kruskal.test; name = "chi-squared"
  } else { # if (method == "oneway")
    func = oneway.test; name = "F"
  }
  old = options(warn = -1)
  ans = sapply(data, function(d) func(d ~ group)$statistic)
  options(old)
  ans = cbind(ans)
  rownames(ans) = colnames(data)
  if (sort) ans = cbind(ans[order(ans[,1], decreasing = TRUE), ])
  colnames(ans) = name
  ans
}

# 関数の実行
# library(kernlab)
# data(spam)
# Rep(spam[c(1:10, 16, 19)], spam$type)
# Rep(spam[c(1:10, 16, 19)], spam$type, sort=FALSE)
# Rep(spam[c(1:10, 16, 19)], spam$type, method = "t")
# Rep(spam[c(1:10, 16, 19)], spam$type, method = "t", sort=FALSE)

Rep(iris[1:4], iris$Species, method="k") # method は省略形でも指定できる
Rep(iris[1:4], iris$Species, method="k", sort=FALSE)
Rep(iris[1:4], iris[,5], method="o") # group はベクトルでなければならない
Rep(iris[1:4], iris[,5], method="o", sort=FALSE)

 

コメント

プログラム書法

2015年03月03日 | ブログラミング

まことにもっていらぬお節介の極みであるが,「二群間で差のある変数を特定する」で提示されている 2 つの関数をまとめて,データとグループ変数の指定を少し柔軟にし,結果をソートするかしないかの引数を追加したり...
ほとんど原型をとどめなくなちゃったけど,キレイにはなったかな。
「多群間で差のある変数を特定する例」でも,この関数を下敷きに使うこともできよう。

# 全ての変数にウィルコクソンの順位和検定または t 検定を行う関数
Rep <- function(data, group, method = c("wilcox", "t"), sort = TRUE) {
  method <- match.arg(method)
  old <- options(warn = -1)
  ans <- if (method == "wilcox")
    sapply(data, function(d) wilcox.test(d ~ group, paired = FALSE)$statistic)
  else
    sapply(data, function(d) t.test(d ~ group, paired = FALSE, var.equal = FALSE)$statistic)
  options(old)
  ans <- cbind(ans)
  rownames(ans) <- colnames(data)
  if (sort) ans <- cbind(ans[order(ans[,1], decreasing=TRUE), ])
  colnames(ans) <-  ifelse(method == "wilcox",  "W", "t")
  ans
}

# データの読み込み
library(kernlab)
data(spam)

# 関数の実行(データの指定)
Rep(spam[c(1:10, 16, 19)], spam$type)
Rep(spam[c(1:10, 16, 19)], spam$type, sort=FALSE)
Rep(spam[c(1:10, 16, 19)], spam$type, method = "t")
Rep(spam[c(1:10, 16, 19)], spam$type, method = "t", sort=FALSE)

コメント

探索を逆方向から攻める

2015年03月02日 | ブログラミング

整数の平方根をとったとき,小数点以下 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 秒

コメント

素数の分布

2015年03月02日 | ブログラミング

素数の存在範囲(上限)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:*

コメント