プロジェクトオイラー 問８０別解

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

問８０「平方根の10進展開」
よく知られているように、自然数の平方根が整数ではないならば、それは無理数である。
そのような平方根の10進展開(decimal expansion)は繰り返しを持たない無限小数になる。
2の平方根は1.41421356237309504880..., であり、
その頭から100桁の数字を合計すると475になる。

初めの100個の自然数の平方根のうち、無理数について、
それぞれの頭から100桁の数字を足した数の総和を求めよ。


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






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

def f(n, m):
	L, c = [i*i for i in xrange(n+1)], 0
	for i in xrange(n+1):
		if i in L: continue
		x, a = i*10**m, i*10**(m*2)
		while x*x-a>1: x = (x+a/x)/2
		c += sum([int(s) for s in str(x)[:m]])
	return c
n, m = 100, 100
print f(n, m)

平方根をニュートン法で求めます。
ニュートン法は方程式の解を近似計算する数値計算法です。
ニュートン法ではf(x)=0となるxを求める際に、xの近傍のxnでの接線とそのX軸との交点xn+1を順に求めていくことを繰り返してxの近似値を求めます。

（参考）
・ニュートン法（wikipedia）

y=f(x)において、x=xnでの接線の傾きは微分係数f'(xn)なので以下の式が成り立ちます。
f(xn) = f'(xn) × (xn - xn+1)
式変形して、
xn+1 = Xn - f(xn)/f'(x)

当問題では平方根を求める元の数をaとして
f(x) = x2 - a としてf(x)=0となるxの近似値を求めます。
なお、xもaもそのままの値では小数点以下の桁が丸められてしまうので10のm乗倍などして大きな整数値にしてから処理して、x2 - a が1を超える場合に近似値を求め続けていきます。

開平法と比べ、ループ変数iの内側が簡単な式になっています。
処理時間も私の作った開平法ロジックと比べて30%程度です。


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

40886