仕事の部屋

2015年11月8日日曜日

151108(2)

Ruby

Polite number(2)

politeness について、
オンライン整数列大辞典の

A069283（http://oeis.org/A069283/list）
と比較し、答え合わせしてみる。

require 'prime'

def politeness(n)
  return 0 if n < 2
  factor = Prime.prime_division(n)
  o_factor = factor.clone
  o_factor.shift if o_factor[0][0] == 2
  # 奇数の約数の個数 - 1
  o_factor.inject(1){|a, n| a * (n[1] + 1)} - 1
end

def A069283(n)
  (0..n).map{|i| politeness(i)}
end
ary = A069283(100)

# OEIS A069283のデータ
ary0 =
[0,0,0,1,0,1,1,1,0,2,1,1,1,1,1,3,0,1,2,1,1,3,1,1,
 1,2,1,3,1,1,3,1,0,3,1,3,2,1,1,3,1,1,3,1,1,5,1,1,1,
 2,2,3,1,1,3,3,1,3,1,1,3,1,1,5,0,3,3,1,1,3,3,1,2,1,
 1,5,1,3,3,1,1,4,1,1,3,3,1,3,1,1,5,3,1,3,1,3,1,1,2,
 5,2]
# 一致の確認
p ary == ary0

151108

Ruby

約数の出力

Polite number を考えていて、約数の出力について気になった。

haruyaさんのコード（http://d.hatena.ne.jp/haruya12/20150918/1442570781）

がシンプルで勉強になった。

require 'prime'

def d(m, k)
  return m unless @factor[k]
  (0..@factor[k][1]).map{|e| d(m * @factor[k][0] ** e, k + 1)}.flatten
end

n = 630
# 順番を逆にする方が後の結果が綺麗になる
@factor = Prime.prime_division(n).reverse
p @factor
p d(1, 0)
p d(1, n % 2 == 0 ? 1 : 0)
p d(1, n % 3 == 0 ? 2 : 0)
p d(1, n % 5 == 0 ? 3 : 0)
p d(1, n % 7 == 0 ? 4 : 0)

出力結果
[[7, 1], [5, 1], [3, 2], [2, 1]]
[1, 2, 3, 6, 9, 18, 5, 10, 15, 30, 45, 90, 7, 14, 21, 42, 63, 126, 35, 70, 105, 210, 315, 630]
[1, 2, 3, 6, 9, 18, 5, 10, 15, 30, 45, 90]
[1, 2, 3, 6, 9, 18]
[1, 2]
1

2015年11月7日土曜日

151107(3)

Ruby

Polite number(1)

2の冪を除けばよく、

f(n) = n + log(n + log(n, 2), 2).to_i

（https://en.wikipedia.org/wiki/Polite_number）

と表される。
オンライン整数列大辞典の
A138591（http://oeis.org/A138591/list）
と比較し、答え合わせしてみる。

include Math

def f(n)
  n + log(n + log(n, 2), 2).to_i
end

def A138591(n)
  (1..n).map{|i| f(i)}
end
ary = A138591(71)

# OEIS A138591のデータ
ary0 =
[1,3,5,6,7,9,10,11,12,13,14,15,17,18,19,20,21,22,
 23,24,25,26,27,28,29,30,31,33,34,35,36,37,38,39,
 40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,
 56,57,58,59,60,61,62,63,65,66,67,68,69,70,71,72,
 73,74,75,76,77]
# 一致の確認
p ary == ary0

151107(2)

Ruby

(-1)^k (n / k) の和

オンライン整数列大辞典の
A059851（http://oeis.org/A059851/list）
と比較し、答え合わせしてみる。

def s(n)
  (1..n).inject(0){|s, k| s - (n / k) * (-1) ** k}
end

def A059851(n)
  (0..n).map{|i| s(i)}
end
ary = A059851(74)

# OEIS A059851のデータ
ary0 =
[0,1,1,3,2,4,4,6,4,7,7,9,7,9,9,13,10,12,12,14,12,
 16,16,18,14,17,17,21,19,21,21,23,19,23,23,27,24,
 26,26,30,26,28,28,30,28,34,34,36,30,33,33,37,35,
 37,37,41,37,41,41,43,39,41,41,47,42,46,46,48,46,
 50,50,52,46,48,48]
# 一致の確認
p ary == ary0

151107

Ruby

n / k の和

オンライン整数列大辞典の
A006218（http://oeis.org/A006218/list）
と比較し、答え合わせしてみる。

def s(n)
  (1..n).inject(0){|s, k| s + n / k}
end

def A006218(n)
  (0..n).map{|i| s(i)}
end
ary = A006218(59)

# OEIS A006218のデータ
ary0 =
[0,1,3,5,8,10,14,16,20,23,27,29,35,37,41,45,50,52,
 58,60,66,70,74,76,84,87,91,95,101,103,111,113,119,
 123,127,131,140,142,146,150,158,160,168,170,176,
 182,186,188,198,201,207,211,217,219,227,231,239,
 243,247,249]
# 一致の確認
p ary == ary0

2015年11月4日水曜日

151104

Ruby

Kolakoski sequence(2)

短くしてみた。

def A000002(n)
  ary = [1, 2, 2]
  i = 2
  while ary.size < n
    ary += [1 + i % 2] * ary[i]
    i += 1
  end
  ary[0..n - 1]
end
ary = A000002(108)

# OEIS A000002のデータ
ary0 =
[1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,2,1,2,1,1,2,
 1,2,2,1,1,2,1,1,2,1,2,2,1,2,2,1,1,2,1,2,2,1,2,1,1,
 2,1,1,2,2,1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,1,2,
 2,1,2,1,1,2,2,1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,
 2,1,2,1,1,2,1,2,2]
# 一致の確認
p ary == ary0

2015年11月3日火曜日

151103(2)

Ruby

Kolakoski sequence(1)

オンライン整数列大辞典の
A000002（http://oeis.org/A000002/list）
と比較し、答え合わせしてみる。

def A000002(n)
  ary = [1, 2]
  b_ary = [2]
  # 次の文字
  str = 1
  ary += b_ary
  while ary.size < n
    f_ary, b_ary = b_ary, []
    (0..f_ary.size - 1).each{|i|
      b_ary += [str] * f_ary[i]
      # 次の文字
      str = 1 + str % 2
    }
    ary += b_ary
  end
  ary[0..n - 1]
end
ary = A000002(108)

# OEIS A000002のデータ
ary0 =
[1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,2,1,2,1,1,2,
 1,2,2,1,1,2,1,1,2,1,2,2,1,2,2,1,1,2,1,2,2,1,2,1,1,
 2,1,1,2,2,1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,1,2,
 2,1,2,1,1,2,2,1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,
 2,1,2,1,1,2,1,2,2]
# 一致の確認
p ary == ary0

151103

Ruby

Gauss circle problem(4)

f(10 ** i) を調べてみた。
iが9でも計算が合ったが、速度がかなり遅くなった。
（実行時間は7時間弱ほどかかる。）

require 'bigdecimal/math'
include BigMath

def f(r)
  q = r * r
  m = (BigDecimal.new('2').sqrt(12) * r / 2).to_i
  # 正方形の内部またはx軸,y軸上の格子点
  s = (2 * m + 1) * (2 * m + 1) + (r - m) * 4
  # 上記以外
  (m + 1..r).each{|x| s += 8 * (BigDecimal.new((q - x * x).to_s).sqrt(12)).to_i}
  s
end

(0..9).each{|i| p f(10 ** i)}

出力結果
5
317
31417
3141549
314159053
31415925457
3141592649625
314159265350589
31415926535867961
3141592653589764829

2015年11月2日月曜日

151102

Ruby

Gauss circle problem(3)

f(10 ** i) を調べてみた。
計算は速いが、iが9で合わなくなった。
Math.sqrtの計算の精度の問題だろう。

def f(r)
  q = r * r
  m = (Math.sqrt(2) * r / 2).to_i
  # 正方形の内部またはx軸,y軸上の格子点
  s = (2 * m + 1) * (2 * m + 1) + (r - m) * 4
  # 上記以外
  (m + 1..r).each{|x| s += 8 * Math.sqrt(q - x * x).to_i}
  s
end

# f(10 ** 9)の計算結果は間違い
(0..9).each{|i| p f(10 ** i)}

出力結果
5
317
31417
3141549
314159053
31415925457
3141592649625
314159265350589
31415926535867961
3141592653589765277

2015年11月1日日曜日

151101

Ruby

Gauss circle problem(2)

Gauss circle problem(1)のコードはnを大きくすると遅くなるので、
別の方法で求めるようにしてみた。
オンライン整数列大辞典の
A000328（http://oeis.org/A000328/list）
と比較し、答え合わせしてみる。

def f(r)
  q = r * r
  m = (Math.sqrt(2) * r / 2).to_i
  # 正方形の内部またはx軸,y軸上の格子点
  s = (2 * m + 1) * (2 * m + 1) + (r - m) * 4
  # 上記以外
  (m + 1..r).each{|x| s += 8 * Math.sqrt(q - x * x).to_i}
  s
end

def A000328(n)
  (0..n).map{|i| f(i)}
end
ary = A000328(46)

# OEIS A000328のデータ
ary0 =
[1,5,13,29,49,81,113,149,197,253,317,377,441,529,
 613,709,797,901,1009,1129,1257,1373,1517,1653,
 1793,1961,2121,2289,2453,2629,2821,3001,3209,3409,
 3625,3853,4053,4293,4513,4777,5025,5261,5525,5789,
 6077,6361,6625]
# 一致の確認
p ary == ary0

登録: 投稿 (Atom)