プロジェクトオイラー 問１２５

プロジェクトオイラー(Project Euler)の問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索です。

問１２５「連続する平方数の和で表せる回文数の和」
回文数 595 は連続する平方数の和で表すことができるという面白い性質を持つ:
6² + 7² + 8² + 9² + 10² + 11² + 12²

1000 未満で連続する平方数の和で表せる回文数はちょうど 11 あり, その合計は 4164 である.
正の整数の平方のみをこの問題では扱うため, 1 = 0² + 1² は含めないことに注意せよ.

回文数でありかつ連続する平方数の和で表せる, 10⁸ 未満のすべての数の合計を求めよ.



注意！！！
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。



私の回答例は以下の通りです。

def f(nmax):
	L = []
	M = [i*i for i in xrange(1, int(nmax ** 0.5)+1)]
	
	for i in xrange(len(M)):
		s = M[i]
		for j in M[i+1:]:
			s += j
			if s >= nmax: break
			t = str(s)
			if t == t[::-1]: L.append(s)
			
	L = set(L)
	return sum(L), len(L)
	
if __name__=="__main__":
	nmax = 10**8
	s, k = f(nmax)
	print "値n ＜", nmax, "の場合：合計 =", s, "、個数 =", k

平方数の和で表せる数を計算し、それが回文数ならためるようにします。

１．関数f(nmax)
・値nの最大値を受け取り、平方数の和で表せる回文数の合計値と個数を返します。

・M = [i*i for i in xrange(1, int(nmax ** 0.5)+1)]
　Mは、必要な平方数の二乗値のリストです。
　二乗値がnmax以下なので、二乗する値iはnmaxの平方根の整数部分までです。
　forループ上の終端値はこの値を +1 した値となります。

・for i in xrange(len(M)):
　ループ変数iは、二乗値リストMの要素のうち、和の計算で使用する最初の値のインデックス値です。

・s = M[i]
　sは二乗値の和です。初期値は二乗値リストMのi番目の値です。

・for j in M[i+1:]:
　　s += j
　ループ変数jは、二乗値リストMのインデックスi番目より後の値です。
　sにjを累積していきます。

・if s >= nmax: break
　二乗値の和sの最大値チェックです。nmax以上になったら処理終了です。
  
・t = str(s)
　if t == t[::-1]: L.append(s)
　二乗値の和sの回文数チェックです。
　二乗値の和sを文字型に変換してtとし、各桁を逆順に並べた文字列と一致するかを判定します。この回文数判定は問３６のときと同じです。
　回文数の場合、リストLにためます。

・L = set(L)
　リストLをset文で集合型にすることで、その要素の値をユニークにします。

・return sum(L), len(L)
　リストLの合計値と個数を返却します。


解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
値n ＜ 100000000 の場合：合計 = 2906969179 、個数 = 166

プロジェクトオイラー 問１２４

プロジェクトオイラー(Project Euler)の問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索です。

問１２４「順序付き根基」
n の"根基"(radical)は, rad(n) で表し, n の素因数の積を意味する.
例えば 504 = 23 × 32 × 7 であるから, rad(504) = 2 × 3 × 7 = 42 である.
1 ≦ n ≦ 10 に対して rad(n) を計算し, rad(n) を対象にソートし, rad(n) が同じ場合は n を対象にソートすると以下のようになる.

未ソート 　 ソート済み
 n rad(n) n rad(n) k
 1  1     1  1     1
 2  2     2  2     2
 3  3     4  2     3
 4  2     8  2     4
 5  5     3  3     5
 6  6     9  3     6
 7  7     5  5     7
 8  2     6  6     8
 9  3     7  7     9
10 10    10 10    10

E(k) をソートした表の n の列の k 番目の要素とする.
例えば, E(4) = 8, E(6) = 9 である.

rad(n) を 1 ≦ n ≦ 100000 でソートした場合, E(10000) を求めよ.



注意！！！
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。



私の回答例は以下の通りです。

def f(nmin, nmax, k):
	P = p(nmax)
	E = []
	for n in xrange(nmin, nmax + 1):
		if n < 1: continue
		L = [1]
		m = n
		for i in P:
			if m % i == 0:
				L.append(i)
				while m % i == 0: m //= i
			if m == 1: break
			
		r = mul(L)
		E.append((r, n, L[1:]))
		
	Es = sorted(E)
	return Es[k-1]

def p(n):
	L = [0,0]+[1]*(n-1)
	i = 2
	while i*i<=n:
		while not L[i]: i += 1
		for j in xrange(i+i, n+1, i): L[j] = 0
		i += 1
	return [i for i in xrange(n+1) if L[i]]

def mul(L): return reduce(lambda x,y:x*y, L)

if __name__=="__main__":
	nmin, nmax, k = 1, 100000, 10000
	r, n, L = f(nmin, nmax, k)
	print nmin, "≦ n ≦", nmax, "のとき、E(", k, ") =", n
	print "rad(", n, ") =", " x ".join([str(i) for i in L]), "=", r 

１．関数f(nmin, nmax, k)
・nの最小値と最大値、問題で求めるk番目として指定するkを受け取り、
　E(k)をソートした表の、k番目の(根基r, n, nの素数リスト)のタプルを返します。

・P = p(nmax)
　Pは必要な素数のリストです。nmax以下の素数が小さい順に入っています。

・for n in xrange(nmin, nmax + 1):
　ループ変数nは、nmimからnmaxまで１つずつ大きくしながら処理します。

・if n < 1: continue
　nが1未満は処理対象外なので処理せず次のnへ進みます。
 
・L = [1]
　Lはnを構成する素因数リストです。
　根基の計算の掛け算のため、1を設定して初期化しておきます

・m = n
　mはnを素因数で割った商です。素因数が判明する都度、値が変化します。

・for i in P:
　ループ変数iは素因数リストPから小さい順に取り出した値です。

・if m % i == 0:
　　L.append(i)
　　while m % i == 0: m //= i
　mが素数iで割り切れる場合、素因数リストLに追加し、mを素因数iで割れるだけ割り込んでおきます。

・if m == 1: break
　mが1になれば素因数分解の終了で、素数iのループを終了します。

・r = mul(L)
　rはnの根基です。関数mulで素数リストLの全要素の積を求めます。
 
・E.append((r, n, L[1:]))
　Eは、(nの根基r, n, nの素因数リストL)のタプルをためるリストです。
 
・Es = sorted(E)
　Esは上記Eに求める範囲のタプルをため終わったところで、ソートします。
　ソートキーは、指定しない場合は先頭からになるので、
　この場合、根基r、nの順にソートキーとして採用されます

・return Es[k-1]
　求めるk番目の値は、ソート済リストEsのインデックスk-1番目の値です。

２．関数P(n)
・n以下の素数を小さい順にもつリストを返します。問３７で作成した関数を流用しました。

３．関数mul(L)
・リストLの全要素の積を返します。問１１０で作成した関数を流用しました。


解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
1 ≦ n ≦ 100000 のとき、E( 10000 ) = 21417
rad( 21417 ) = 3 x 11 x 59 = 1947

プロジェクトオイラー 問１２３

プロジェクトオイラー(Project Euler)の問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索です。

問１２３「素数の自乗で割った余り」
p_n を n 番目の素数とする. (p1 = 2, p2 = 3, ...)
r を (p_n - 1)ⁿ + (p_n + 1)ⁿ を p_n² で割った余りとする.

例えば, n = 3 のとき, p₃ = 5 であり, 43 + 63 = 280 ≡ 5 mod 25.
余り r が 10⁹ より大きくなる n の最小値は 7037 である.
余り r が 10¹⁰ より大きくなる最初の n を求めよ.



注意！！！
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。



私の回答例は以下の通りです。

def f(rmin):
	for i, p in enumerate(P()):
		n = i + 1
		if n % 2:
			r = (2*n*p) % (p*p)
			if r > rmin: return n, p, r
		
def P():
	L, i = [2], 1
	yield 2
	while True:
		i += 2
		for j in L:
			if not i%j: break
		else:
			L.append(i)
			yield i
			
if __name__=="__main__":
	rmin = 10**10
	n, p, r = f(rmin)
	print "n =", n, "(素数p =", p, ", 余りr =", r, ")"

１．関数f(rmin)
・最小の余りrminを受け取り、問題に合うn番目の素数pとそのときの余りrについて、n,p,rを返します。

・この問題は問１２０の続きの問題です。
　問１２０では、
(a - 1)ⁿ + (a + 1)ⁿ を n² で割った余りをrとしていましたが、
　本問では、
(p_n - 1)ⁿ + (P_n + 1)ⁿ を p_n²で割った余りをrとしています。

・問１２０の式に、a=p_nを代入すると以下のことがわかります。
　余りrは、nが偶数か奇数かでパターンが異なります。
　nが偶数の場合、余りrはaに関わらず 2 です。
　nが奇数の場合、余りrは 2np_n をp_n²で割った余りです。
　
・pは素数なので、2以外はすべて奇数です。
　題意から答えは明らかに2よりも大きく、必ず奇数です。
　このため、nが偶数のときは処理せず飛ばします。

２．関数P()
・素数を小さい方から昇順に返すジェネレーターです。呼び出しの都度、次の素数を返却します。
・問１０８で作成した関数を改良しました




解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
n = 21035 (素数p = 237737 , 余りr = 10001595590 )

プロジェクトオイラー 問１２２

プロジェクトオイラー(Project Euler)の問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索です。

問１２２「効率的なべき乗計算」
n15を求めるのに最も単純な方法では 14 回の掛け算が必要である.

n × n × ... × n = n¹⁵
しかし2進法を用いれば, 6 回の掛け算で計算できる.

n × n = n²
n² × n² = n⁴
n⁴ × n⁴ = n⁸
n⁸ × n⁴ = n¹²
n¹² × n² = n¹⁴
n¹⁴ × n = n¹⁵
ところがたった 5 回の掛け算のみでも計算できる.

n × n = n²
n² × n = n³
n³ × n³ = n⁶
n⁶ × n⁶ = n¹²
n¹² × n³ = n¹⁵
m(k)を n^k を求めるのに必要最低限な掛け算の回数と定義する.
たとえば m(15)=5 である.

1 ≦ k ≦ 200 に対し, Σ m(k) を求めよ.



注意！！！
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。



私の回答例は以下の通りです。

def f(n1, n2):
	L = [[] for i in xrange(n2 + 1)]
	M = [0 for i in xrange(n2 + 1)]
	L[1].append([])
	
	for i in xrange(1, n2):
		for m in L[i]:
			s = m + [i]
			t = len(s)
			for j in s[::-1]:
				for k in s[::-1]:
					u = j + k
					if not (i < u  <= n2): break
					if 0 < M[u] < t: continue
					if s in L[u]: continue
					L[u].append(s)
					M[u] = t
	return sum(M[n1:])
	
if __name__=="__main__":
	print f(1, 200)

掛け算するときの掛ける数の指数を足していき、指定数に達するまでの最小回数を求めます。
１度掛け算して作成した累乗値の指数は再利用できます。

掛け算で使用した指数値のユニークな配列を作成し、
この要素と掛け算で到達したその指数値自身から、
値を２つ取り出し、その和の累乗値の配列として順に計算していきます。
このとき、指数値のユニークな値の配列は複数パターンできることがあり、いずれも必要です。
例えば、n=5の場合、到達するパターンは[1,2,3]と[1,2,4]の２つあり、
これを元にさらに大きな指数値について計算していくとき、
どちらの系統を使った方が最小回数になるか、この段階では判断つかないためです。

また、掛け算を実施した回数を別のリストに保存しておき、最後にその合計を求めます。

１．関数f(n1, n2)
・題意にあるm(k)について、n1≦k≦n2のときのm(k)の和を返します。

・L = [[] for i in xrange(n2 + 1)]
　M = [0 for i in xrange(n2 + 1)]
　Lは、そのインデックス位置の値の累積値を求めるときの「掛け合わせた指数値のユニークリスト」のリストです。
　Mは、そのインデックス位置の値の累積値を求めるときの最小掛け算回数のリストです。
　例えば、n=5の場合、L[5]=[[1, 2, 3], [1, 2, 4]]、M[5]=3　です。

・L[1].append([])
　初期値です。これでn**1の場合に、この指数値までに使用した指数値のリストは無いが、この指数値そのもの（つまり、1）を使用して累乗していくことが可能になります。

・for i in xrange(1, n2):
　ループ変数iは、1からn2に至るインデックスです。累乗値計算の範囲です。

・for m in L[i]:
　ループ変数mは「掛け合わせた指数値のユニークリスト」です。

・s = m + [i]
　sは「掛け合わせた指数値のユニークリスト」に、それより到達した指数値iそのものを追加したリストです。
　次の累乗値計算の元になる指数値のリストです。

・t = len(s)
　tは、リストsの次の累乗で生成できる累乗値までの掛け算回数の候補です。

・for j in s[::-1]:
　for k in s[::-1]:
　ループ変数jとループ変数kは、リストsから取り出した２つの要素です。
　大きい値から順に、次の値を計算していきます。

・u = j + k
　uは、次の掛け算で到達できる累乗値の指数です。

・if not (i < u  <= n2): break
　i < u < n2 なら書き換え対象です。
　そうでない場合は、kはそれ以下の値しかないので、条件は満たせずループは終了です。

・if 0 < M[u] < t: continue
　0 < M[u] < tの場合、今までの結果よりも掛け算回数が少ないパターンなので、次のパターンを考えます。
　そうでない場合、ここまでの処理結果よりも同じか長い掛け算回数のパターンなので処理続行です。

・if s in L[u]: continue
　リストsは指数値uのときの「掛け合わせた指数値のユニークリスト」に既にあれば、次のパターンを考えます。

・L[u].append(s)
　M[u] = t
　ここまでくれば題意の満たす初めての「掛け合わせた指数値のユニークリスト」になので、掛け算回数tとともに保存します。

・return sum(M[n1:])
　iのループが終了したところで、n1以上の掛け算回数の和を返却します。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
1582

プロジェクトオイラー 問１２１

プロジェクトオイラー(Project Euler)の問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索です。

問１２１「円盤ゲームの賞金額」
袋の中に赤い円盤 1 枚と青い円盤 1 枚が入っている.
ある賭けゲームにおいて, プレイヤーは円盤を 1 枚ランダムに取り出しその色を記録する.
各ターンの終わりに円盤は袋に戻され, 追加の赤い円盤 1 枚が袋に足され, そして次の円盤がランダムに取り出される.
プレイヤーはゲームをプレイするのに 1 ポンドを支払い, ゲーム終了時に青い円盤を赤い円盤より多く取り出していれば勝利する.
ゲームが 4 ターンプレイされたとすると, プレイヤーが勝利する確率はちょうど11/120 となる.
したがって, 胴元が赤字の見込みになるまでにこのゲームで勝ちに対して割り当てるべき賞金の最大は 10 ポンドとなるであろう.
支払いはすべてポンドの整数倍であり,
またゲームをプレイするのに支払われたもともとの 1 ポンドを含んでいるため,
与えられた例では実際にはプレイヤーは 9 ポンドを獲得することに注意しよう.
15 ターンがプレイされるゲーム 1 回に割り当てられるべき賞金の最大を求めよ.


注意！！！
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。



私の回答例は以下の通りです。

def mul(L): return reduce(lambda x,y: x*y, L)

def f(n):
	import itertools
	m = n//2 + 1
	c = 0
	L = [i+2 for i in xrange(n)]
	Ls = set(L)
	
	for i in xrange(m, n+1):
		for t in itertools.combinations(L, i):
			u = Ls - set(t)
			if u: c += mul([s-1 for s in u])
			else: c += 1
			
	s = mul(L)//c
	return s
	
if __name__=="__main__":
	print f(15)

１．関数f(n)
・n回のターンで得られる賞金の最大値を求めます。

・m = n//2 + 1
　mは勝利に必要な勝ち数です。n回のターンのうち半分以上で勝ちです。

・L = [i+2 for i in xrange(n)]
　Lは各ターンごとの全円盤数のリストです。
　リストLのi番目の要素は、iターン目に円盤を取り出す確率の分母でもあります。

・Ls = set(L)
　Lsは全円盤数リストLを集合型にしたものです。後で集合演算で使います。

・for i in xrange(m, n+1):
　ループ変数iは勝ち数です。
　必要な最低勝ち数mから全勝の場合のnまでの値を取ります。

・for t in itertools.combinations(L, i):
　ループ変数tは勝ち数iの場合に青円盤を取り出すターンでの全円盤数、の組み合わせです。
　つまり、青円盤を取り出すターンごとのそれぞれの確率の分母の組み合わせでもあり、
　各ターンごとの全円盤数リストLから、勝ち数i個の要素を取り出した組み合わせです。

・u = Ls - set(t)
　uは全円盤数リストLsからループ変数tの要素を集合演算で差し引いたリストです。
　つまり、赤円盤を取り出すターンでの全円盤数、の組み合わせです。
　これは、青円盤を取り出す１つのパターンの中で、異なる赤円盤を取り出す。場合の数です。
　なので、uの要素を全部掛け合わせると、青円盤を取り出す１つのパターンの。場合の数になります。

・if u: c += mul([s-1 for s in u])
　else: c += 1
　cは青円盤を取り出す１つのパターンの、場合の数の累積値です。
　青円盤を取り出す1つのパターンの場合の数は、uの要素を全部掛け合わせた値なので
　これを累積します。
　ただし、全部青円盤、つまり全ターンで勝つ場合、赤円盤が出ないためリストuが空で
積として計算できないので固定で１パターン分を累積します。

・s = mul(L)//c
　return s
　求める値sは１ゲーム当たりの最大賞金で、勝率の逆数です。
　全部のパターン数を、勝ちパターン数で割った値の小数点以下切り捨て値です。
　問題中の例の４ターンの場合、勝率11/120なので、最大賞金は120/11 = 10.9 ...で10です。
　全部の場合の数は、各ターンごとの全円盤数の積です。
　これを勝つパターン数cで割った値の整数値が、求める最大賞金です。

２．関数mul(L)
・リストLの全要素の積です。問１１０で作成した関数を流用しました。



解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
2269

プロジェクトオイラー 問１２０

プロジェクトオイラー(Project Euler)の問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索です。

問１２０「自乗で割った余り」
(a-1)ⁿ + (a+1)ⁿ を a² で割った余りを r と定義する.
例えば, a=7, n=3 の時, r=42 である: 6³ + 8³ = 728 ≡ 42 mod 49.
n が変われば r も変わるが, a=7 の時 r の最大値 rmax は 42 であることがわかる.
3 ≦ a ≦ 1000 において, Σ rmax を求めよ.


注意！！！
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。



私の回答例は以下の通りです。

def f(s, e): return sum([a*(a - (2- a%2)) for a in xrange(s, e+1)])
start, end = 3, 1000
print f(start, end)

まず、(a+1)ⁿ + (a-1)ⁿ　を展開してみます。
(a+1)ⁿ + (a-1)ⁿ ÷ a² の余りについて、手計算して様子を見てみます。

n=2のとき
  (a+1)² + (a-1)²
= (a² + 2*a + 1) + (a² - 2*a + 1)
= 2(a² + 1)
なので、a²で割った余りは、2です。

n=3のとき
  (a+1)³ + (a-1)³
= (a³ + 3*a² + 3*a + 1) + (a³ - 3*a² + 3*a - 1)
= 2(a³ + 3*a)
= 2(aⁿ * n*a)
なので、a²で割った余りは、2*n*aのa²での余り、2*n*a÷a²の余りです。

n=4のとき
  (a+1)⁴ + (a-1)⁴
= (a⁴ + 4*a³ + 6*a² + 4*a + 1) + (a⁴ - 4*a³ + 6*a² - 4*a + 1)
= 2(a⁴ + 6*a²) + 2
なので、a²での余りは、2です。

上記より、展開式をaの乗数の多い項から順に並べたとき、
・aの係数は二項係数でパスカルの三角形のそれぞれの値です。
・(a+1)と(a-1)のn乗同士をたすことで、偶数番目の項は互いに打ち消し合います。
・題意からn≧3なので、aの2乗以上の項は全部割り切れ、余りはでません。
なので、
nが奇数の場合、余りはa¹の項だけみればよく、2*n*a÷(a²)の余りです。
nが偶数の場合、余りはa⁰の項だけみればよく、2です。
なので、nが奇数の場合だけ考えれば十分です。
割り算の余りは「割られる数」＝「割る数の倍数」のときに0になり、
割られる数がこの値よりも小さいができるだけ大きな値というのが求める値です。
また、「割られる数」が「割る数」を超えると、
それまでに出現した余りが再度循環して出現するので、
「割られる数」が「割る数」以下の場合だけ考えれば十分です。
特に、「割られる数」＝「割る数」の場合、商が1で余り0。
「割られる数」がこの値よりも小さいができるだけ大きな値の場合、
商が0で残りが余りになり、「割られる数」＝「余り」になります。
そこでまず、「割られる数」＝「余り」の場合のnを計算します。
2an = a²
n = a/2
nがこの値よりも小さくできるだけ大きな値のとき、余りが最大になります。
ところで、題意からnもaも整数なので、aが奇数か偶数かで場合分けが必要です。
(1)aが奇数の場合
n = a/2よりも小さくできるだけ大きな値は、n = (a-1)/2 です。
このときの求める値は、2an = 2a(a-1)/2 = a(a-1) です。
(2)aが偶数の場合
n = a/2よりも小さくできるだけ大きな値は、n = (a-2)/2 です。
このときの求める値は、2an = 2a(a-2)/2 = a(a-2) です。
aを奇数か偶数かで場合分けして上記の値を累計すればいいことになります。

１．関数f(s, e)
s以上e以下の値について、問題文にある最大値を返します。
・sum([a*(a - (2- a%2)) for a in xrange(s, e+1)])
　%演算子は余りなので、
　「a%2」は、aが奇数なら1、偶数なら0です。
　「2 - a%2」は、aが奇数なら1、偶数なら2です。
　「a*(a - (2- a%2))」は、aが奇数ならa(a-1)、aが偶数ならa(a-2)　です。
　sからeまで順にaのforループで上記の計算結果を配列に入れて、sum関数で合計します。
　ループは内包表記をしています。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
333082500