2012年7月26日木曜日

プロジェクトオイラー Problem27

「プロジェクト オイラー」の問題をpythonでやってみます。

出典: Project Euler (日本語翻訳サイト)
参考:
サイエンスもの置き場 プロジェクトオイラーを遊び倒すガイド(初級編)

Problem 27
オイラーは以下の二次式を考案している:
n2 + n + 41
この式は, nを0から39までの連続する整数としたときに40個の素数を生成する.
しかし, n = 40のとき402 + 40 + 41 = 40(40 + 1) + 41となり41で割り切れる.
また, n = 41のときは412 + 41 + 41であり明らかに41で割り切れる.


計算機を用いて, 二次式 n2 - 79n + 1601という式が発見できた.
これはn = 0 から 79 の連続する整数で素数を生成する.
係数の積は, -79 × 1601 で -126479である.

さて, |a| < 1000, |b| < 1000 として以下の二次式を考える.

n2 + an + b (ここで|a|は絶対値)
n=0から始めて連続する整数で素数を生成したときに最長の長さとなる上の二次式の,係数a, bの積を答えよ.



私の解答例は以下です。
-----
def p(n):
	if n<=2: return []
	L = [2]
	for t in xrange(3, n+1, 2):
		for s in L:
			if not t%s: break
			elif t<s*s: L.append(t); break
	return L
	
def g(a, b, n, L):
	for i in xrange(n):
		s = i*i+a*i+b
		if s not in L: break
	return i
	
def f(n):
	L, M = p(n), [0, 0, 0]
	for b in L+[-s for s in L]:
		for a in xrange(-999, 1000, 2):
			if abs(1+a+b) not in L: continue
			s = g(a, b, n, L)
			if M[2]<s: M = [a, b, s]
	return M
	
a, b, s = f(1000)
print "ab =",a*b,", n2",a,"n +",b,", the number of prime =", s
-----

総当りするにはa,bの範囲が広すぎるので、範囲を絞り込むために、
f(n) = n^2 + a*n + b
に小さい値を代入して様子を見ます。

・f(0)=b
 bは素数 ・・・(A)
 素数は「2かそれ以上の奇数」の集合に含まれるので、
 bは「2かそれ以上の奇数」に含まれる ・・・(B)

・f(1)=1+a+b
 「1+a+b」は素数 ・・・(C)
 (C)に(B)を当てはめると、
 aは奇数 ・・・(D)

これらをふまえて、

1.関数p(n)
・n未満の素数リストを返します。
problem 010を改良しました。

2.関数g(a, b, n, L)
・aとbは「n^2 + a*n + b」の係数、nはこの式のnの上限値です。
 Lは範囲内の素数のプラスマイナス両符号分のリストです。

・a,b,nによる計算値が、素数リストLに無い値になったときにループを抜け、
 そのときのインデックス値を返します。
 このインデックス値が、計算値のうち、素数が続いた数になります。

3.関数f(n)
・for b in L+[-s for s in L]:
 bは、上記(A)より、n未満の素数のプラスマイナス両符号分を範囲とします。
 リストは+演算子で、両方の要素を合わせたリストになります。
 
・for a in xrange(-999, 1000, 2):
 aは、上記(D)より、絶対値1000未満の奇数を範囲とします。
 
・if abs(1+a+b) not in L: continue
 範囲を絞ったa,bを総当りして、上記(C)より、「1+a+b」が対象素数リストに無ければ次の候補へ処理と進めます。

・s = g(a, b, n, L)
 対象素数リストにあれば、そのときのa,bは条件を満たすので関数gで、計算値に
 何個素数が続いて現れるかチェックします。、
 
・if M[2]<s: M = [a, b, s]
 連続して出現する素数の個数が最大の場合、リストMにa,bの値とともに素数個数も保存します。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
ab = -59231 , n2 -61 n + 971 , the number of prime = 62

0 件のコメント:

コメントを投稿