规矩出方圆

思路

在单位正方形上有一个内切圆（圆心 $(0.5, 0.5)$ ，半径 $0.5$ ）。我们放 $M$ 条水平和竖直切割线（含边界），将正方形切成若干小矩形。只要一个矩形和圆的交集面积 $> 10^{-10}$ ，就算"被染色"。目标是让被染色的小矩形总面积最小。

先想几个关键问题：

什么样的矩形会被"浪费"？ 被染色的矩形分两类——完全在圆内的（不可避免，面积贡献 = 圆面积），和横跨圆边界的"边界格子"。边界格子的整块面积被计入，但实际只有一部分落在圆内。我们的目标就是让这些边界格子尽可能薄。

怎么利用对称性？ 圆关于中心 $(0.5, 0.5)$ 有四重对称性。最优方案可以让横、纵切割线用同一组坐标，并且关于 $0.5$ 对称。

核心转化——角度参数化。 设 $k = (M-1)/2$ ，对称切割点位于 $c_i = 0.5 - 0.5\cos(\alpha_i)$ ，其中 $\alpha_0 < \alpha_1 < \cdots < \alpha_{k-1}$ 是 $(0, \pi/2)$ 中的角度。这等价于把圆弧上的点投影到坐标轴。每条切割线正好穿过圆上某个角度对应的点。

最优的自洽条件。 当切割线位于 $c_i$ 时，圆在 $y = c_i$ 处的两个交点的 $x$ 坐标分别是 $0.5 \pm \sqrt{c_i(1 - c_i)} = 0.5 \pm 0.5\sin(\alpha_i)$ 。如果令 $\alpha_i + \alpha_{k-1-i} = \pi/2$ （角度关于 $\pi/4$ 对称），那么 $\sin(\alpha_i) = \cos(\alpha_{k-1-i})$ ，圆的交点恰好落在另一条切割线上！这意味着边界格子在角落处恰好"贴着"圆弧。

在这个对称约束下，一个格子 $(i, j)$ 位于圆外（不被染色）当且仅当 $i + j \le k - 1$ 。这使得"未着色面积"的表达式变得非常干净：

$ $\text{未着色} = 4 \sum_{i=0}^{k-1} w_i \cdot \text{prefix}[k - i]$ $

其中 $w_i = 0.5(\cos\alpha_{i-1} - \cos\alpha_i)$ 是第 $i$ 条带的宽度， $\text{prefix}[m]$ 是前 $m$ 条带宽度之和。

解析梯度 + L-BFGS。 由于角度对称，自由变量只有 $\lfloor k/2 \rfloor$ 个（ $k$ 为奇数时中间角度固定为 $\pi/4$ ）。对每个自由角度 $\beta_j$ ，梯度可以解析推导：

$ $\frac{\partial \text{val}}{\partial \beta_j} = 4(w_j \cos\beta_j - w_{k-1-j} \sin\beta_j)$ $

用 L-BFGS 算法几十次迭代就能收敛到最优解，对 $M = 199$ 也只需不到 0.1 秒。

代码

import math

def solve():
    M = int(input())
    k = (M - 1) // 2
    n_free = k // 2
    k_odd = k % 2
    pi4 = math.pi / 4
    cos_pi4 = 0.7071067811865476

    if n_free == 0:
        if k == 0:
            print("1.0000")
        elif k == 1:
            w = 0.5 * (1.0 - cos_pi4)
            print(f"{1.0 - 4.0 * w * w:.4f}")
        elif k == 2:
            pass
        else:
            print("1.0000")
        if k <= 1:
            return

    n = n_free if n_free > 0 else 1
    if k == 2:
        n = 1

    def eval_grad(params):
        sb = sorted(params)
        cos_b = [math.cos(b) for b in sb]
        sin_b = [math.sin(b) for b in sb]

        cos_a = list(cos_b)
        if k_odd:
            cos_a.append(cos_pi4)
        for i in range(n - 1, -1, -1):
            cos_a.append(sin_b[i])

        prev = 1.0
        prefix = [0.0] * (k + 1)
        widths = [0.0] * k
        for i in range(k):
            widths[i] = 0.5 * (prev - cos_a[i])
            prefix[i + 1] = prefix[i] + widths[i]
            prev = cos_a[i]

        uc = 0.0
        prev = 1.0
        for i in range(k):
            uc += widths[i] * prefix[k - i]
        val = 1.0 - 4.0 * uc

        sort_order = sorted(range(n), key=lambda i: params[i])
        grad_sorted = [0.0] * n
        for j in range(n):
            grad_sorted[j] = 4.0 * (widths[j] * cos_b[j] - widths[k-1-j] * sin_b[j])

        grad = [0.0] * n
        for i, si in enumerate(sort_order):
            grad[si] = grad_sorted[i]

        return val, grad

    params = [(i + 1) * pi4 / (n + 1) for i in range(n)]

    lo_bound = 1e-10
    hi_bound = pi4 - 1e-10

    m_lbfgs = min(10, n)
    s_list = []
    y_list = []
    rho_list = []

    val, grad = eval_grad(params)

    for iteration in range(5000):
        if math.sqrt(sum(g*g for g in grad)) < 1e-13:
            break

        q = list(grad)
        alphas_bfgs = []
        for i in range(len(s_list) - 1, -1, -1):
            alpha = rho_list[i] * sum(s_list[i][j]*q[j] for j in range(n))
            alphas_bfgs.append(alpha)
            for j in range(n):
                q[j] -= alpha * y_list[i][j]
        alphas_bfgs.reverse()

        if s_list:
            ys = sum(y_list[-1][j]*s_list[-1][j] for j in range(n))
            yy = sum(y_list[-1][j]**2 for j in range(n))
            gamma = ys / yy if yy > 0 else 1.0
        else:
            gamma = 0.01

        r = [gamma * qi for qi in q]

        for i in range(len(s_list)):
            beta = rho_list[i] * sum(y_list[i][j]*r[j] for j in range(n))
            for j in range(n):
                r[j] += s_list[i][j] * (alphas_bfgs[i] - beta)

        direction = [-ri for ri in r]

        step = 1.0
        old_val = val
        for ls in range(30):
            new_params = [max(lo_bound, min(hi_bound, params[i] + step * direction[i])) for i in range(n)]
            new_val, new_grad = eval_grad(new_params)
            if new_val < old_val - 1e-15 * abs(old_val):
                break
            step *= 0.5
        else:
            break

        s = [new_params[j] - params[j] for j in range(n)]
        y = [new_grad[j] - grad[j] for j in range(n)]
        sy = sum(s[j]*y[j] for j in range(n))

        if sy > 1e-15:
            if len(s_list) >= m_lbfgs:
                s_list.pop(0)
                y_list.pop(0)
                rho_list.pop(0)
            s_list.append(s)
            y_list.append(y)
            rho_list.append(1.0 / sy)

        params = new_params
        val = new_val
        grad = new_grad

    print(f"{val:.4f}")

solve()

复杂度分析

时间复杂度：每次 L-BFGS 迭代需要 $O(k)$ 计算函数值和梯度，总迭代次数约为 $O(k)$ 级别，故总复杂度约 $O(k^2)$ ，其中 $k = (M-1)/2$ 。
空间复杂度： $O(k)$ ，存储角度参数、宽度数组和 L-BFGS 历史向量。