Project Euler 226

题目

A Scoop of Blancmange

The blancmange curve is the set of points \((x, y)\) such that \(0 \le x \le 1\) and \(y = \sum \limits_{n = 0}^{\infty} {\dfrac{s(2^n x)}{2^n}}\), where \(s(x)\) is the distance from \(x\) to the nearest integer.

The area under the blancmange curve is equal to \(\dfrac{1}{2}\), shown in pink in the diagram below.

Let \(C\) be the circle with centre \((\dfrac{1}{4},\dfrac{1}{2})\) and radius \(\dfrac{1}{4}\), shown in black in the diagram.

What area under the blancmange curve is enclosed by \(C\)?

Give your answer rounded to eight decimal places in the form \(0.abcdefgh\)

解决方案

该页面给出了求解牛奶冻曲线积分的递推公式。如果\(f(x)=\sum \limits_{n = 0}^{\infty} {\dfrac{s(2^n x)}{2^n}}\)，那么有\(I(x)\)：

\[I(x)=\int_{0}^x f(x)dx= \left \{\begin{aligned} &\dfrac{I(2x)}{4}+\dfrac{x^2}{2} & & \text{if}\quad 0\le x\le \dfrac{1}{2} \\ &\dfrac{1}{2}-I(1-x) & & \text{if}\quad \dfrac{1}{2}\le x\le 1 \\ &\dfrac{n}{2}+I(x-n) & & \text{if}\quad n\le x\le(n+1) \\ \end{aligned}\right. \]

为方便计算，我们假设圆\(C\)和曲线中只有两个交点，其中右交点在\(B(0.5,f(0.5))\)，左交点通过二分法来求出，假设为为\(A(a,f(a))\)。

接下来使用该页面的一些公式计算了弓形\(\mathop{AB}\limits^{\frown}\)的面积，从而计算出弧\(\mathop{AB}\limits^{\frown}\)竖直向下的那一块面积\(S\)，也就是说，

\[S=\int_a^{\frac{1}{2}} \frac{1}{2} - \sqrt{\frac{1}{4}^2 - \left(x - \frac{1}{2}\right)^2}dx\]

最终答案为\(I(\dfrac{1}{2})-I(a)-S\)。

代码

from math import asin, sin

f = lambda x: sum(abs(x * (1 << n) - round(x * (1 << n))) / (1 << n) for n in range(60))
g = lambda x: 0.5 - (x / 2 - x * x) ** 0.5


def I(x):
    if x < 1e-13:
        return 0
    if x <= 0.5:
        return I(2 * x) / 4 + x * x / 2
    else:
        return 0.5 - I(1 - x)


l, r = 0, 0.5
for _ in range(100):
    mid = (l + r) * 0.5
    if f(mid) > g(mid):
        r = mid
    else:
        l = mid
xl, xr = l, 0.5
yl, yr = g(xl), g(xr)
cr = 1 / 4
d = ((xl - xr) ** 2 + (yl - yr) ** 2) ** 0.5
theta = 2 * asin(d / (2 * cr))
arc_area = cr * cr / 2 * (theta - sin(theta))
circle_integral = (yl + yr) * (xr - xl) / 2 - arc_area
ans = I(xr) - I(xl) - circle_integral
print("{:.8f}".format(ans))