Project Euler 226 - inamori’s diary

プロジェクトオイラー
 http://projecteuler.net/

Q226.
$y = \sum_{n=0}^{\infty}{\frac{s(2^nx)}{2^n}}$ （s(x)はxから最も近い整数までの距離）が描く曲線の下の部分と、中心(1/4, 1/2）半径1/4の円の重なる面積を求めよ。

交点を求めて、あとは台形法で面積を求めた。
フラクタルについて利用できないかと思ったが、精度の問題があって意味がないように思えた。



from math import sqrt, asin, pi
def s(x):

    return 0.5 - abs(0.5 - x)
def f(x):

    y = 0

    n = 0

    while x > 0:

        y += s(x) / 2 ** n

        x *= 2

        x = x - int(x)

        n += 1

    return y
def C(x):

    return 0.5 - sqrt(x * (0.5 - x))
def cross(f, g, x0, x1):

    for n in xrange(N - 3):

        x2 = (x0 + x1) / 2

        if g(x2) < f(x2):

            x1 = x2

        else:

            x0 = x2

    return (x0 + x1) / 2
def integral(f, x0, x1, n):

    s = 0

    x = x0

    h = (x1 - x0) / n

    for k in xrange(n + 1):

        s += f(x)

        x += h

    s = (s + s - f(x0) - f(x1)) * h / 2

    return s
def g(x):

    return x
N = 2 ** 16

x0 = cross(f, C, 0, 0.125)

S1 = integral(f, x0, 0.5, N)

r = 0.25

S2 = (0.5 - x0) * 0.5 - r ** 2 / 2 * (pi / 2 + asin( (0.25 - x0) / r

                )) - (0.25 - x0) / 2 * sqrt(r ** 2 - (0.25 - x0) ** 2)

print S1 - S2