脳log[ProjectEuler: 2011-03-15～]

最終更新: 2012-02-29T17:34+0900

♪ [ProjectEuler] Problem 72

 Problem 72

前のバージョンはここ(20110312p01.02)。step2を最適化*すると倍くらい速くなって 1分半(ウチのPC基準で。C++だとコンマ1秒)。それに伴い素因数の組み合わせを求める必要がなくなったので、100万要素の配列の配列が 100万要素の Fixnum配列になってメモリ使用量が再び数MBレベルになった。

LIMIT = 1_000_000
pfs = Array.new(LIMIT+1, 1)
count = LIMIT*(LIMIT+1)/2 - 1 # 1) 2からLIMIT以下の分母 d につき d 通りの分子を予め計上する。
2.upto(LIMIT){|d| # 0) 2からLIMIT以下の分母 d に対して、
	print d,"\r"
	d.step(LIMIT, d){|_| pfs[_] *= -d } if pfs[d] == 1
	# 2) 分子が d で分母が d の倍数になるものを加減する。
	count += pfs[d]/pfs[d].abs * (LIMIT/d)*(LIMIT/d+1)/2 if pfs[d].abs == d
}
p count

一週間後にはこのコードが理解できなくなってること請けあい。

ググった>「Dreamshire | Project Euler Problem 72 Solution」。φ関数の値を合計するとかで、Rubyでも10秒未満。Problem 69に出てきた関数だけど、その問題はインチキしたから理解してないんよね。

LIMIT = 1000000
phi = (0..LIMIT).to_a
2.upto(LIMIT){|n|
	n.step(LIMIT, n){|m|
		phi[m] *= (n-1.0)/n
	} if phi[n] == n
}
p phi.inject(&:+).floor-1

繰り返しの構造は似てるのにこの実行時間の差はなんだ？ と思ったら print d,"\r" の有無だけだった。俺のも同等に速い。

最終更新: 2011-03-14T23:14+0900

♪ [ProjectEuler] Problem 73, 74

 Problem 73

昨日(20110312p01)の延長。Problem 72の解答と同じ部分より違う部分を見つける方が難しい。

LIMIT = 12000
pfs = Array.new(LIMIT+1){ [] }
count = 0
2.upto(LIMIT){|d| # 0) LIMIT以下のすべての分母 d に対して、
	print d,"\r"
	d.step(LIMIT, d){|_| pfs[_] << d } if pfs[d].empty?
	n_min, n_max = d/3, (d-1)/2 # (n_min, n_max]
	next unless n_min < n_max
	count += n_max - n_min # 1) とりあえず一通りの分子を計上し、
	# 2) 8分の6など通分可能なものを差し引きする。
	(1..(pfs[d].size)).each{|r|
		cms = pfs[d].combination(r).map{|pf| pf.inject(&:*) }
		count -= (-1)**(r%2+1) * cms.map{|cm| n_max/cm - n_min/cm }.inject(&:+)
	}
}
p count

 Problem 74

問題文のヒントを最大限に利用したが数分かかる。総当たりでなく組み合わせ単位でテストしてその順列を計上したらマシになるかも。順列を考えるときに先頭の桁に 0を置くようなミスを犯しそうだがね。

factorial = [1,1] # [0!,1!,...]
factorial.push(factorial.size*factorial.last) until 9 < factorial.size

chain_length = lambda{
	memo = {
		169 => 3, 363601 => 3, 1454 => 3,
		871 => 2, 45361 => 2,
		872 => 2, 45362 => 2
	}
	f = lambda{|start|
		return memo[start] if memo.has_key?(start)
		next_ = start.to_s.chars.map{|c| factorial[c[0]-?0] }.inject(&:+)
		return memo[start] = 1 + (start == next_ ? 0 : f.call(next_))
	}
}.call

p (1...1_000_000).inject(0){|sum,n| sum += 1 if chain_length.call(n) == 60; sum }

最終更新: 2011-03-15T00:15+0900

♪ [ProjectEuler][Ruby] Problem 71, 72

 Problem 71

「HCF(n,d)=1」には用語の説明があると思ったんだけどなかった。n/dの nと dは最大公約数が 1の既約分数だとすると意味がとおるので HCF=Highest Common Factorだと決めた(でっちあげ)。

分数とはなんぞやだとか切断だとか小難しく考えてしまったが(実際には考えられるほど知らない)、ワンライナーだった。3/7より少し小さい100万個の分数を小数になおして、一番小さいものを見つける。有理数にして比較しないのは時間がかかるから。公約数をみつけたりする時間だろうか。

require 'rational'
p Rational(*(2..1_000_000).inject([0,1]){|answer,d|
	answer[0]/answer[1].to_f < (d*3-1)/7/d.to_f ? [(d*3-1)/7,d] : answer
})

 @2011-03-14 なにこれ！

Project Euler Problem #71 « KeyZero Conversation

分数を初めてならった小学生が必ず間違える分数の足し算(通分せずに分母どうし分子どうしを加算する)にこんな意味があるとか！

 Problem 72

分単位のお時間がかかります。(訳:一時間はかからないけど……)

何倍も速くなるので「Integer#prime_division」を使う代わりに 100万要素の配列を使ってる。トレードオフで使用メモリは数MBから 100MB超になるが。 小手先のチューンよりアルゴリズムを改良しろってのはもっともだけど、かなしいかな、できることとできないことがあるのです。

LIMIT = 1_000_000
pfs = Array.new(LIMIT+1){ [] }
count = 0
2.upto(LIMIT){|d| # 0) LIMIT以下のすべての分母 d に対して、
	print d,"\r"
	count += d-1 # 1) とりあえず d-1 通りの分子を計上し、
	d.step(LIMIT, d){|_| pfs[_] << d } if pfs[d].empty?
	# 2) 8分の6など通分可能なものを差し引きする。
	(1..(pfs[d].size)).each{|r|
		cms = pfs[d].combination(r).map{|pf| pf.inject(&:*) }
		count -= (-1)**(r%2+1) * cms.map{|cm| (d-1)/cm }.inject(&:+)
	}
}
p count

Ruby 1.9からバックポートされてきた(のだと思われる見覚えのないメソッド) cycle, tap, combination, permutation, productといったメソッドが便利だ。あとは自然数を無限に生成し続ける無限リストのようなものをどれだけ簡単に書けるかだ。なにかショートカットがあるのだろうか。これでは長すぎる。

Enumerable::Enumerator.new(lambda{|&block| n=0; loop{ block.call n+=1 } }, :call).each{|x| p x }

それと、block.callの部分を yieldにできないのもわかりにくい。Procと blockと lambdaの微妙な違いによるものなのだろうか。

Rubyによる他所の Project Eulerの解答をみていてこういう書き方も知ってるけど、カウンタが Floatになっちゃうのが不満。

1.upto(1/0.0){|n| p n }

あれ？ Fixnumだ。Floatになるのは stepだった。

1.step(1/0.0){|n| p n } # 1.0, 2.0, 3.0,...

明示的に Fixnumの増分: 1を指定しても n は Float. この違いはなんだろう。

cycleの使い道として zipを想定していたが拒否されてしまった。

irb> [1,2,3,4,5].zip([0])
=> 1, 0], [2, nil], [3, nil], [4, nil], [5, nil
irb> [1,2,3,4,5].zip([0].cycle)
TypeError: can't convert Enumerable::Enumerator into Array
        from (irb):2:in `zip'
        from (irb):2
        from :0
irb> RUBY_DESCRIPTION
=> "ruby 1.8.7 (2010-01-10 patchlevel 249) [i386-mswin32]"

最終更新: 2011-03-12T02:24+0900

♪ [ProjectEuler] 65, 66, 67, 68, 69

 Problem 65

分母を一番深いところから順番に計算していく。

a = ([2] + (1..33).map{|k| [1,2*k,1] }.inject(&:+)).reverse
denom, numer = *a.inject([0,1]){|nd, x| [nd[1], x*nd[1]+nd[0]] }
require 'rational'
p Rational(numer, denom).numerator.to_s.chars.inject(0){|sum,c| sum - ?0 + c[0] }

 Problem 66

xを増やしながらの総当たりで、最後に見つかった Dが答え。と思ったんだけど Dが見つかるペースがどんどこ落ちていく。一日以上かけても 969個の Dのうち 270個が残ってる。

 Problem 67

Problem 18の延長で以前解いた。

 Problem 68

1. 16桁だから 10は一回しか使わない。10は external node.
2. 先頭の桁(外側のノードの最小値)を最大にするために、external nodeに 10,9,8,7,6を配置する。
3. そうすると一辺の和はとりうる値の中で最小の 14。
4. 二桁目を最大にする 6-5-3からスタートして考えよう。次は 10-3-1しかない。
5. ってなかんじで穴埋め。

 Problem 69

前問に引き続いて、数学でもプログラミングでもなく、オラクルで。

1. n/φ(n)を最大にするために……
2. 分母を小さくするために、nはたくさんの素数を因数に持っている方がいい。⇒すべて異なる素数の積からなる数。
3. 分子を大きくするために、nは上限の 100万に近い方がいい。
4. 分子と分母のバランスの取り方は勘で。

最終更新: 2011-03-12T03:46+0900

♪ [ProjectEuler] Q62, Q63, Q64

 Q62

"exactly five" って書いてあるから、同じ数字の並べ替えで作れる立方数が 6以上あってもダメだと思うんだ。

memo = Hash.new{|h,k| h[k] = [] }
n = 0
k_length = 1
loop{
	n += 1
	cube = n*n*n
	k = cube.to_s.split(//).sort.join('')
	if k.length != k_length
		answer = memo.values.select{|cubes| cubes.size == 5 }
		if not answer.empty?
			answer.each{|cubes| p cubes }
			exit
		end
		memo.clear
		k_length = k.length
	end
	memo[k] << cube
}

 Q63

# (x-1)/x <= log10(n) < 1 (n = ?,?,...)
count = 0
x = 0
loop{
	x += 1
	boundary = (x-1.0)/x
	lower_bound = (1..9).to_a.reverse.find{|n| Math.log10(n) < boundary } || 0
	count += 9 - lower_bound
	break if lower_bound == 9
}
p count

またまた「Dreamshire | Project Euler Problem 63 Solution」の解答を検討してみたい。二重のループなんてない。logの計算だって 9回だけ。どういうことだ？

1. 10の n乗は常に n+1桁になる。
2. 9の n乗は 10の n乗より小さいため n+1桁には絶対に届かず、ある程度までは n桁を保つが、そのうち n-1桁に落ちる。
3. 8の n乗は 9の n乗より早く n-1桁に落ちる。
4. 以下 1の n乗まで。
5. で、「ある程度」って具体的には？
6. 10を約0.954乗すると 9になる。9は 10より 0.046(=1-0.954)程度小さい数だ。
7. この 0.046がいくつ集まると 10一個分小さい(=桁が落ちる)ことになるだろう。21.7(=1÷0.046)だ。
8. 9の場合、21乗までは n乗が n桁を保っているが 22乗は違う。

というストーリーをひねり出した。「9は 10より 0.046(=1-0.954)程度小さい数だ」ってくだりがいかにも苦しい。小数だからごまかしがきいてるけど、ぴったり 10一個分小さくなる場合は n桁、n-1桁、どっち？ (たぶんまだ n桁だな。1^1がそう)

ともあれ、明かされてみればワンライナーの問題だったよ。

p (1..9).inject(0){|sum,n| sum + (1/(1-Math.log10(n))).floor }

常用対数を直接求めるメソッドが用意されてるあたりが Rubyだなとおもた。

 Q64

連分数っていうらしい。a_nの求め方、a += 1 while 0 <= r - (n - d*(a+1))**2 の条件部分が判然としない。スクリプト中のコメントにあるように、対象としてるルートの係数が必ず約分されて 1になることも理解できてない。

def next_frac(r, n, d) # (√r + n) / d = a + 1 / [(√r + n_) / d_]
	a = 0
	a += 1 while 0 <= r - (n - d*(a+1))**2
	d_ = (r - (n - d*a)**2) / d
	raise if (r - (n - d*a)**2) % d != 0 # why OK?
	n_ = -(n - d*a)
	return a, r, n_, d_
end

def period_of(r)
	rnd = [r, 0, 1]
	arr = []
	loop{
		a, *rnd = next_frac(*rnd)
		arr << rnd
		period = arr.size-1 - arr.index(rnd)
		return period if 0 != period
		return 0 if rnd[2] == 0 # √r is rational.
	}
end

count = 0
1.upto(10000){|n|
	count += 1 if period_of(n) % 2 == 1
}
p count

ところで、この問題を解くときに Math.sqrtを使うのってインチキくさくない？(だから使ってないんだけど)

最終更新: 2011-03-02T05:24+0900

♪ [ProjectEuler] Q61

 Q61

何も考えずにコーディングしただけ。一瞬 CPUが考え込みます。

generators = [
	lambda{ n = 0
		lambda{ n+=1; n*(n+1)/2 }
	}.call,
	lambda{ n = 0
		lambda{ n+=1; n*n }
	}.call,
	lambda{ n = 0
		lambda{ n+=1; n*(3*n-1)/2 }
	}.call,
	lambda{ n = 0
		lambda{ n+=1; n*(2*n-1) }
	}.call,
	lambda{ n = 0
		lambda{ n+=1; n*(5*n-3)/2 }
	}.call,
	lambda{ n = 0
		lambda{ n+=1; n*(3*n-2) }
	}.call,
]
# 数を準備
d4polynumbers = generators.map{|g|
	() while (p = g.call) < 1000
	a = [p]
	a.push(p) while (p = g.call) < 10000
	a
}
# 端緒(の集まり)
bunch_of_chain = d4polynumbers[d4polynumbers.size-1].map{|p|
	[[p, d4polynumbers.size-1]]
}
# 端緒を伸ばすもの
extender = lambda{|chain, pool|
	xx = chain.last.first.to_s[-2,2]
	( (0...(pool.size)).to_a - chain.map{|_| _.last } ).map{|i|
		[i, pool[i]]
	}.map{|i, nums|
		nums.find_all{|num|
			num.to_s[0,2] == xx
		}.map{|num|
			chain + [[num, i]]
		}
	}.inject(&:+)
}
# 伸ばしていく
(d4polynumbers.size-1).times{
	bunch_of_chain = bunch_of_chain.map{|chain|
		extender[chain, d4polynumbers]
	}.inject(&:+)
}
# 輪っか？
bunch_of_cyclic_chain = bunch_of_chain.reject{|chain|
	chain.first.first.to_s[0,2] != chain.last.first.to_s[-2,2]
}
# 出力
bunch_of_cyclic_chain.each{|chain|
	puts chain.map{|a,_| a }.join("\t")
	puts chain.map{|_,b| "P#{b+3}" }.join("\t")
	puts "sum: #{chain.map{|a,_| a }.inject(&:+)}"
}

先は長いのにもう失速してる。「良いもの。悪いもの。: Project Eulerを100問解いてみた」テトレーションとか聞いたこともない単語なんだけど……。

中学生の時に 3^{50} の一の位は何かという問題が出た。でも Problem 188は何乗したらいいかもわからない。下手の考え休むに似たりっていうけどどうしたもんかなあ。ない知恵を絞るのも悪くないと思うんだけど。

最終更新: 2011-02-20T21:45+0900

♪ [ProjectEuler] Q58, Q59, Q60

 Q58

10%未満っていうのは絶妙なポイントなのかな。全然 9%未満に落ちない。

def prime? x
	return false if x < 2
	return true if x == 2
	quo, rem = x.divmod(2)
	return false if rem == 0
	t = 1
	while t < quo
		t += 2
		quo, rem = x.divmod(t)
		return false if rem == 0
	end
	return true
end

x, t = 1, 0
primes_on_diagonals = 0
loop{
	t += 2
	3.times{
		x += t
		primes_on_diagonals += 1 if prime? x
	}
	x += t
	puts "#{primes_on_diagonals} primes out of #{2*t+1} (#{100*primes_on_diagonals/(2*t+1)}%, side length=#{t+1})"
	exit if 100 * primes_on_diagonals / (2*t+1) < 10
}

 Q59

1. 暗号化キーの長さが 3文字なのがわかってるので、数列を 3列に並べて各列を眺める。
2. 英単語には eが一番使われるとか、単語では theが再頻出だとかの統計データがあるらしい。(と書いてあるのを何度も目にしたが実際のデータは見たことがない)
3. 文章なら文字としてはスペースが一番多いはずだ。
4. 単語で theが一番多いなら、theの前後の空白は右下に向いて、(1行1列目)↘(2行2列目), (1行2列目)↘(2行3列目), (1行3列目)↘(2行1列目) の並びで現れるはずだ。
5. というかんじでマジックナンバーが出てきた。
6. あ、問題文のサマリに「Using a brute force attack, can you ...」って書いてある。しまった。

encrypted_text = [79,59,12,...,22,73,0,0] # last 2 elements are padding.
text = ""
0.step(encrypted_text.size-1, 3){|i|
	text += (encrypted_text[i+0] ^ (71 ^ " "[0])).chr
	text += (encrypted_text[i+1] ^ (79 ^ " "[0])).chr
	text += (encrypted_text[i+2] ^ (68 ^ " "[0])).chr
}
text.chop!.chop! # remove padding
puts text
puts "sum: #{text.bytes.inject(:+)}"

 Q60

1. 素数二つ(A,B)に分割できてローテートしても素数な素数(AB,BA)。⇒ a set of 2 primes (A,B)
2. (AX,XA),(BX,XB)が 1の条件を満たす 2 sets of 2 primes が見つかる。⇒ a set of 3 primes (A,B,X)
3. (AY,YA),(BY,YB),(XY,YX)が 1の条件を満たす 3 sets of 2 primesが見つかる。⇒ a set of 4 primes (A,B,X,Y)
4. (AZ,ZA),(BZ,ZB),(XZ,ZX),(YZ,ZY)が 1の条件を満たす 4 sets of 2 primesが見つかる。⇒ a set of 5 primes (A,B,X,Y,Z)

1を満たす素数を発見しながらそれを使って、1の集合から2へ、2の集合から3へ、3の集合から4へ、要素をプロモートしていけばよさそう。

# 寝る前にやる。

寝てしまった。答えが出ない。素数を分割するんでなく、素数のペアを組み合わせて素数かどうか判定した方がいいかもしれない。そろそろ身にしみて理解してきたけど、素数って印象よりありふれ過ぎてる。

 @2011-02-17

ちょっとくらい時間がかかってもいーやって考えてたけど、何日もかけても四つ組みが 7つと、五つ組が 0個しか見つからないことがわかったので、1分以内に答えを出すべくもうちょっと考える。

1. 小中学生の頃に読んだ数字遊びの本から得た知識によれば、10進表記の各桁の和が 3の倍数の数はそれ自体も 3の倍数だという。だというwww
2. 各桁の和を 3で割った余りが 0の素数は 3だけ。
3. 各桁の和を 3で割った余りが 1の素数と 2の素数を組にすると、それをつなげた数は 3の倍数になってしまい素数ではなくなるので組にできない。
4. 各桁の和を 3で割った余りで素数を分類すると 5つ組みとして考えられるのは (0,1,1,1,1), (0,2,2,2,2), (1,1,1,1,1), (2,2,2,2,2)の 4パターンだけ。
5. 3桁までの素数で総当たりしてみよう。
6. 見つからなかったので 4桁までで総当たり。

def prime? x
	return false if x < 2
	return true if x == 2
	quo, rem = x.divmod(2)
	return false if rem == 0
	t = 1
	while t < quo
		t += 2
		quo, rem = x.divmod(t)
		return false if rem == 0
	end
	return true
end

set012 = [[],[3],[]]
require 'mathn'
Prime.new.each{|prime|
	break if 10000 <= prime
	dmod3 = prime.to_s.bytes.inject(0){|sum,byte| sum+byte-?0 } % 3
	set012[dmod3] << prime
}
set1, set2 = set012[1], set012[2]
set2[0] = 3
# set1 = [3,7,13,...]
# set2 = [3,5,11,...]

make_group_of_two = lambda{|set|
	pair = {}
	0.upto(set.size-2){|i|
		(i+1).upto(set.size-1){|j|
			if prime?("#{set[i]}#{set[j]}".to_i) and prime?("#{set[j]}#{set[i]}".to_i) 
				(pair[[set[i]]]||=[]) << set[j]
			end
		}
	}
	return pair
}
group1, group2 = make_group_of_two.call(set1), make_group_of_two.call(set2)

extend_group = lambda{|g|
	group = {}
	g.each_pair{|rest, last1s| # rest + one of last1s = group
		last1s.each{|last1|
			next1s = last1s
			gg, out = rest.clone, last1
			gg.size.times{|i|
				gg[gg.size-1-i], out = out, gg[gg.size-1-i]
				next1s &= g[gg]||[]
			}
			if ! next1s.empty?
				group[rest+[last1]] = next1s
			end
		}
	}
	return group
}
group1, group2 = extend_group.call(group1), extend_group.call(group2) # sets of 3 primes
group1, group2 = extend_group.call(group1), extend_group.call(group2) # sets of 4 primes
group1, group2 = extend_group.call(group1), extend_group.call(group2) # sets of 5 primes

printer = lambda{|rest, last1s|
	last1s.each{|last1|
		puts %[#{rest.inject(&:+)+last1}:\t#{rest.join("\t")}\t#{last1}]
	}
}
group1.each(&printer)
group2.each(&printer)

分単位の時間で答えはでたけどもその五つ組の合計が意外に大きくて、10000以上の素数を組に加えても最小の組み合わせになりうる。計算量の増大の仕方がひどくて、これ以上桁数を増やして試行するのは無理だというのに。

 「だという わらわらわら」@昨日

じゃないよね。

\begin{array}{rcl} q & = & a_0 + 10a_1 + 10^2a_2 +……+ 10^na_n \quad\mbox{(}a_n\mbox{は 0以上 9以下の整数)}\\ & = & (a_0 + a_1 + a_2 +……+ a_n) + 9a_1 + 99a_2 +……+(10^n-1)a_n\\ \end{array}

a_0+a_1+a_2+……+a_n が 3の倍数の整数 qは 3の倍数です。

たしか 4の倍数についても同じような判定規則があった気がした。忘れたけど。

たしか 5の倍数についてもどこかの桁を見るだけで(略

 @2011-02-19

4は 100を作るから下2桁だけ。5は 10を作るから下1桁だけを見ればいい。

 「計算量の増大の仕方がひどくて」

一番時間を食ってるのは make_group_of_two. 異なる二要素の組み合わせということで n^2-n 回の素数判定を行ってる。素数判定自体も nの大きさに比例する(※1:1ではないけど)ループを持っている。大変なはずだ。

とりあえず、今の素数判定より賢い素数判定があるのはわかってるけどわからないので使ってない。(注:わかる => 知ってる, 理解できる) 丁寧にコードを読んだらわかるかもだけどそれはチートっぽい。大学入試の数論関係の問題だって、解答をチラ見したら誰だって理解できんだよ。

最終更新: 2011-02-12T22:42+0900

♪ [ProjectEuler] Q56, Q57

 Q56

Bignumはできれば使いたくない。aが 100未満なので 8桁ずつ。

answer = [0, 0, 0] # sum, a, b
1.upto(99){|a|
	digits = "1"
	1.upto(99){|b|
		sum = 0
		carry = 0
		0.step(digits.size-1, 8) {|i|
			l, r = [0, digits.size-i-8].max, digits.size-i
			carry, digits8 = (digits[l...r].to_i * a + carry).divmod(100000000)
			digits8 = "00000000#{digits8}"[-8,8]
			digits[l...r] = digits8
			digits8.each_byte{|byte|
				sum += byte - ?0
			}
		}
		if carry != 0
			digits8 = carry.to_s
			digits = digits8 + digits
			digits8.each_byte{|byte|
				sum += byte - ?0
			}
		end
		if answer[0] < sum
			*answer = sum, a, b
		end
	}
}
p answer

 Q57

とかいいながら Bignum。

count = 0
numer, denom = 1, 1
1000.times{
	numer, denom = numer + denom + denom, numer + denom
	count += 1 if numer.to_s.length != denom.to_s.length
}
p count

最終更新: 2011-02-10T04:56+0900

♪ [ProjectEuler] Q51, Q52, Q53, Q55

 Q51

ただただ、手と CPUを動かすだけで精一杯。(頭は役に立ってないよ)

primes = [2]
is_prime = lambda{|x|
	result = true
	primes.each{|prime|
		quo, rem = x.divmod(prime)
		if rem == 0
			result = false
			break
		end
		break if quo < prime
	}
	return result
}

# replace 2 digits or 3 digits. キ・メ・ウ・チ
def find_8_prime_family(a)
	return [] if a.size < 8
	a.map!{|x| x.to_s }
	h = Hash.new{|h,k| h[k] = [] }

	# 2 digits
	0.upto(a.first.size-3){|i|
		(i+1).upto(a.first.size-2){|j|
			h.clear
			a.each do |prime|
				if prime[i] == prime[j]
					h[prime[0...i]+prime[(i+1)...j]+prime[(j+1)...(prime.size)]].push prime
				end
			end
			h.each do |_,v|
				return v if v.size == 8
			end
		}
	}

	# 3 digits
	0.upto(a.first.size-4){|i|
		(i+1).upto(a.first.size-3){|j|
			(j+1).upto(a.first.size-2){|k|
				h.clear
				a.each do |prime|
					if prime[i] == prime[j] and prime[j] == prime[k]
						h[prime[0...i]+prime[(i+1)...j]+prime[(j+1)...k]+prime[(k+1)...(prime.size)]].push prime
					end
				end
				h.each do |_,v|
					return v if v.size == 8
				end
			}
		}
	}

	return []
end

x = 1
start = 0 # start of primes of a width.
loop {
	x += 2
	next unless is_prime.call x
	print "#{x}\r"
	if primes[start].to_s.length != x.to_s.length
		a = find_8_prime_family primes.last(primes.size-start)
		if ! a.empty?
			puts a.sort.join(" ")
			exit
		end
		start = primes.size
	end
	primes.push x
}

 Q52

桁数ごとに探索範囲を決めて、総当たり。

問題が xについても同じ数の組み合わせであることを求めてると思わなくてチェックしてないけど、結果的に xも 2x,3xなんかと同じ数字で構成されてた。

digits = 10
loop {
	digits *= 10
	(digits/2).upto((digits*10-1)/6){|x|
		print "#{x}\r"
		x2 = (x*2).to_s.split(//).sort
		if [3,4,5,6].all?{|n|
			x2 == (x*n).to_s.split(//).sort
		} then
			puts [1,2,3,4,5,6].map{|n| x*n }.join(" ")
			exit
		end
	}
}

 Q53

浮動小数点数なんてファジーなものを使っちゃったよ。Math.sqrtの使用をこれまで頑なに避けてたのも、結果が Floatになるからだったり。

count = 0
23.upto(100){|n|
	cmb = 1.0
	1.upto(n/2){|r|
		cmb /= r
		cmb *= (n-r+1)
		count += (n-r == r) ? 1 : 2 if 1_000_000 < cmb
	}
}
p count

 Q55

問題文が難しかった。3割ぐらいは推測。

あっけなく答えが出たので to_s.reverse.to_i みたいなのをなくすべく、Integer#reverse を自作してみたら、かえって遅くなったし。

class Integer
	# 負数については考えてない。
	def reverse
		x = 0
		this = self
		begin
			this, rem = this.divmod(10)
			x = 10*x + rem
		end while 0 < this
		x
	end
end

count = 0
10.upto(10_000-1){|x|
	is_lychrel = true
	50.times{
		x = x + x.reverse
		if x == x.reverse
			is_lychrel = false
			break
		end
	}
	count += 1 if is_lychrel
}
p count

最終更新: 2011-02-09T20:16+0900

♪ [ProjectEuler] Q50

 Q50

1. 100万未満の素数リストを作成。
2. リストをコピーして作業領域とする。
3. コピーしたリストの要素に、一つ右の位置の素数、二つ右の位置の素数、以下略を加えていき、都度、和が素数かどうか確かめる。

時間がかかるので逐一進捗を表示してる。この問題に魔法の一手なんてあるのかね。

primes = []
is_prime = lambda{|x|
	result = true
	primes.each{|prime|
		quotient, remainder = x.divmod(prime)
		if remainder == 0
			result = false
			break
		end
		break if quotient < prime
	}
	return result
}
2.upto(999_999){|x|
	primes.push x if is_prime.call x
}
puts "#{primes.size} primes under 1 million."

work = primes.dup
step = 0
primes_found = []
live_elements = work.size
while 0 < live_elements
	step += 1
	primes_found.clear
	live_elements = 0
	print "step #{step}\r"
	0.upto(work.size-1-step){|i|
		work[i] += primes[i+step]
		if work[i] < 1_000_000
			live_elements += 1
			primes_found.push work[i] if is_prime.call work[i]
		end
	}
	if primes_found.empty?
	elsif primes_found.size < 10
		puts "step #{step}: #{primes_found.join ' '}"
	else
		puts "step #{step}: #{primes_found.size} primes"
	end
end

魔法の一手はなくても……

• 素数列の作成と和の評価を同時進行にできるかも。
• 作業領域の範囲の狭め方が下手。
• 奇数ステップ(偶数個の素数の和)は不要。
• 100万に近い素数は不要。(それを一つ加えただけで 100万を超えるんでは項数を稼げない)
• 項数の上限は求められる。

答えを出した後でググるのが楽しい。フォーラムは読んでないけど、多分これ以上ないっていうような答えが書いてありそうで、面白くなさそうな気がしてる。(理解できない数学的知識が使われてたら、print XXXXXXX(answer); って書かれてるのと変わらないから)

最終更新: 2011-02-09T01:05+0900

♪ [ProjectEuler] Q47, Q48, Q49

 Q47

昨日よりちょっとはマシになったかと。アホすぎた素数判定を、素因数の数を数える処理と一体化した。でも 10秒以上かかります。

primes = [2]
have4primefactors = []

num_of_factors = lambda{|x|
	prime_factors = 0
	primes.each{|prime|
		quotient, remainder = x.divmod(prime)
		if quotient < prime
			prime_factors += 1
			break
		end
		if remainder == 0
			prime_factors += 1
			break if 4 < prime_factors
			x /= prime while x % prime == 0
			break if x == 1
		end
	}
	return prime_factors
}

x = 2
loop {
	x += 1
	print "#{x}\r"
	case num_of_factors.call(x)
	when 1
		primes.push x
		have4primefactors.clear
	when 4
		have4primefactors.push x
		p have4primefactors.first if have4primefactors.length == 4
	else
		have4primefactors.clear
	end
}

 Q48

恥ずかしいほどまっすぐで乱暴なスクリプトだけど、コンソールの表示も待てないくらいノーウェイトで答えが出るんだから仕方がない。

p (1..1000).inject(0){|sum,x| sum + x**x }

 Q49

1. 4桁の素数リストを作る
2. リストから条件を満たす二つの素数を選ぶ
3. 計算で求めた三番目の数(三つのうち一番大きい)が条件を満たすか調べる。

10秒くらいかかります。

primes = []

is_prime = lambda{|x|
	result = true
	primes.each{|prime|
		quotient, remainder = x.divmod(prime)
		if remainder == 0
			result = false
			break
		end
		break if quotient < prime
	}
	return result
}

2.upto(9999){|x|
	primes.push x if is_prime.call x
}

primes_4digit = primes.last(primes.length - primes.rindex{|x| x < 1000 } - 1)
0.upto(primes_4digit.size-1){|i|
	p = primes_4digit[i]
#	next if p == 1487
	(i+1).upto(primes_4digit.size-1){|j|
		q = primes_4digit[j]
		r = q + q - p
		next if p.to_s.split(//).sort != q.to_s.split(//).sort or
		        q.to_s.split(//).sort != r.to_s.split(//).sort
		k = nil
		(j+1).upto(primes_4digit.size-1){|_k|
			if r == primes_4digit[_k]
				k = _k
				break
			elsif r < primes_4digit[_k]
				break
			end
		}
		if k
			puts "#{p} #{q} #{r}"
#			exit
		end
	}
}

最終更新: 2011-02-07T05:28+0900

♪ [ProjectEuler] Q46, Q47

 Q46

squares[]はソート済みなのに .include?()でそれを活かせないのが不満。

primes = [] # Omit 2. Even prime is not needed.
squares = [1]
def prime?(n)
	return false if 0 == n%2
	3.step(n/2, 2) {|x|
		return false if 0 == n%x
	}
	return true
end

x = 1
loop {
	x += 2
	squares.push((squares.size+1)**2) if squares.last < x
	if prime?(x)
		primes.push x
		next
	end
	print "#{x}\r"
	next if primes.any?{|prime| prime < x and squares.include?((x-prime)/2) }
	p x # answer
	break
}

 Q47

何の工夫もないのですんごく時間がかかる。

def prime_gt2?(n)
	return false if 0 == n%2
	x, upper_bound = 3, n/2
	while x <= upper_bound
		upper_bound, remainder = n.divmod(x)
		return false if 0 == remainder
		x += 1
	end
	return true
end

primes = [2]
have4primefactors = []
have4primefactor = lambda{|x|
	num_of_factors = 0
	primes.each{|prime|
		if x % prime == 0
			num_of_factors += 1
			break if 4 < num_of_factors
			x /= prime while x % prime == 0
		end
	}
	return num_of_factors == 4
}

x = 2
loop {
	x += 1
	print "#{x}\r"
	if prime_gt2? x
		primes.push x
		have4primefactors.clear
	elsif have4primefactor.call(x)
		have4primefactors.push x
		p have4primefactors.first if have4primefactors.length == 4
	else
		have4primefactors.clear
	end
}

最終更新: 2011-02-05T10:26+0900

♪ [ProjectEuler] Q14, 『珠玉のプログラミング』, Collatz予想

ちょっと前の日記から……

Q14

この漸化式は『珠玉のプログラミング』で見た。どうして収束するのかわからなかった。

見たっていうのはコラム4の問題で。

4.6問題

5.入力xが正の整数であるとき、以下のループが終了することを示してください。

while x != 1 do
  if xが偶数なら
    x = x/2
  else
    x = 3*x+1

これがまんま「コラッツの問題 - Wikipedia)」と呼ばれる未解決の問題だということに、今日「d.y.d.」を読んでいて気がついた。本の巻末のヒントを読み直してみたらこんなことが書いてあるし。

もし、この問題が解けたら、近くの大学の数学科に急いで行って、博士号を申請しましょう。

ひどい(笑)。わからなくて当然だ。遅刻学生が黒板の問題を宿題だと思って解いて提出したら、それは未解決問題だったとかいうのは、お話の世界なんだから。