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

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

問８５「長方形の数え上げ」
注意深く数えると, 横が3, 縦が2の長方形の格子には, 18個の長方形が含まれている.

ぴったり2,000,000個の長方形を含むような長方形の格子は存在しない. 一番近い解を持つような格子の面積を求めよ.






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






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

def f(n):
	xmax, sk, xk, yk = int((-1+(n*8+1)**0.5)/2), 1, 1, 1
	for x in xrange(1, xmax):
		t = x*(x+1)
		ys = int(((-t)+(t*t+16*t*n)**0.5)/2/t)
		for y in [ys, ys+1]:
			s = x*(x+1)*y*(y+1)/4
			if s==n: return x*y, x, y, s
			if abs(n-s)<abs(n-sk): sk, xk, yk = s, x, y
	return xk*yk, xk, yk, sk

n = 2000000
a, x, y, m = f(n)
print "area:", a
print "grid:", x,"x",y
print "num of rectangular:", m

格子が縦1個横x個の場合、含まれる長方形の数nは、1からnまでの和であり、
n = x(x+1)/2
これをxの二次方程式として解くと、
x2 + x - 2n = 0
x = (-1 ± √1-4x1x(-2n)) / 2
x>0より
x = (-1 + √1+8n) / 2　・・・(A)

このxは、縦1個のときの横の個数で、横の個数の最大値xmaxとなります。
こうして、xのループを1からxmaxの値まで1つずつ変化させます。

縦x個横y個の格子に含まれる長方形の個数nは、
縦方向にx(x+1)/2、横方向にy(y+1)/2の積です。
n = x(x+1)/2 × y(y+1)/2
n = x(x+1)y(y+1)/4　・・・(B)
ここで、t = x(x+1)と置くと、
n = ty(y+1)/4
n = (ty2 + ty)/4
ty2 + ty - 4n = 0

yの二次方程式として解くと、
y = (-t ± √(t2-4t(-4n))) / 2t
y>0より
y = (-t + √(t2+16nt)) / 2t
yは整数値なので、上記yの前後の整数、つまり上記yをint関数で切り捨てた値と、これに+1した値のいずれかが横の個数です。
上記の縦横の個数の候補の組み合わせで含まれる長方形の個数を計算し、ぴったりn個ならば終了、それ以外はnとの差の絶対値が小さいものをとっておきます。


１．関数f(n)
・n個の長方形を含む格子に一番近い解を持つ格子の面積、及びそのときの縦の個数、横の個数、長方形の数を返します。
　ただし、縦の個数は最も小さい値とします。

・xmax, sk, xk, yk = int((-1+(n*8+1)**0.5)/2), 1, 1, 1
　初期値設定します。縦個数の最大値xmaxとして上記式(A)を実装します。
　平方根はmathモジュールにsqrt関数もありますが、
　累乗演算子**で0.5乗しても同じです。

・ys = int(((-t)+(t*t+16*t*n)**0.5)/2/t)
　横個数の最大値ysとして上記式(C)を実装します。

・s = x*(x+1)*y*(y+1)/4
　含まれる長方形の数sとして上記式(B)を実装します。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
area: 2772
grid: 36 x 77

num of rectangular: 1999998