REALHW89 规矩出方圆

题目链接

规矩出方圆

题目描述

在边长为 1 的正方形画布上，圆心放在 (0.5, 0.5)，半径为 0.5 的圆恰好内切于画布。

我们要在水平方向与竖直方向各放置 M 条切割线（含边界 0 和 1，总共 M 条），把正方形分成 (M-1)×(M-1) 个小矩形。每个小矩形只要与圆的交集面积大于 1e-10，就视为“被染色”。

当 M 为奇数时，期望通过巧妙地安排切割线的位置，使所有被染色的小矩形的总面积尽可能小，并输出这个最小面积，结果四舍_五入到小数点后 4 位。

为简化问题，可以利用关于中心对称的性质：最优方案可以令横纵切割线共享同一组坐标。

解题思路

本题并非简单的几何模拟，而是一个隐藏的组合优化问题，可以通过动态规划 (DP) 并结合分治优化 (D&C Optimization / Knuth-Yao Speedup) 来精确求解。之前尝试的梯度下降法属于数值近似，难以保证找到全局最优解且实现复杂易错。

1. 问题转化与建模

  1. 对称性：利用对称性，我们只需考虑第一象限，即 均在 的区域，圆心在 。为方便计算，我们将坐标系平移，使圆心在 ，半径 ，此时正方形范围为 。我们只需研究第一象限 的情况，最终结果乘以 4。

  2. 染色面积的理解：题目要求最小化的“染色面积”，可以理解为用 个高度不等的矩形（台阶）去覆盖四分之一圆的面积。我们的目标是调整这些台阶的分割点 ，使得所有台阶的总面积 最小。

  3. 目标函数：台阶的总面积 ，其中 是圆在高度 的半宽度。我们要最小化的其实是台阶面积与圆真实面积的差值，即“超出的面积” ，其中 是四分之一圆的面积。

2. 离散化与动态规划

  1. 离散化：连续空间上的优化是困难的。我们将 轴从 到 的区间离散化为 个（例如 ）等距的微小点 。问题转化为从这 个点中挑选 个点作为切割点。

  2. 动态规划 (DP)：

    • 状态定义: 表示用 个台阶（ 条水平切割线）覆盖到离散点 时，所能达到的最小“超出面积”。
    • 状态转移方程:
    • 是用一个台阶覆盖 到 区间所产生的“超出面积”： 其中 是圆面积从 到 的定积分。

3. 分治优化 (D&C Optimization)

朴素 DP 的复杂度为 ，会超时。但上述 cost 函数满足四边形不等式（或说最优决策点具有单调性），因此可以使用分治优化。

算法流程:

  1. 预处理：计算所有离散点 对应的 值和圆面积的定积分前缀和 。
  2. 主循环：进行 轮迭代，代表放置 条切割线。
  3. 分治函数 :
    • 此函数用于计算 在 属于区间 时的值。
    • 。我们暴力枚举 在 中为 找到最优决策点 。
    • 由于决策单调性， 的最优决策点在 中， 的最优决策点在 中。
    • 递归调用 和 。
    • 这样，每层 DP 的复杂度从 降为 。

最终结果：，即 4 倍的 (四分之一圆面积 + 最小超出面积)。

代码

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>

using namespace std;

const double R = 0.5;
const double INF = 1e100;

int M, P, N;
vector<double> y, v, s;
vector<double> dp0, dp1;

// 成本函数 cost(k, j)
double cost(int k, int j) {
    return (y[j] - y[k]) * v[k] - (s[j] - s[k]);
}

// 分治优化
void dnc(int l, int r, int k_min, int k_max) {
    if (l > r) return;

    int mid = l + (r - l) / 2;
    double best_val = INF;
    int k_best = -1;

    int k_upper = min(k_max, mid - 1);
    for (int k = k_min; k <= k_upper; ++k) {
        double current_val = dp0[k] + cost(k, mid);
        if (current_val < best_val) {
            best_val = current_val;
            k_best = k;
        }
    }
    dp1[mid] = best_val;

    dnc(l, mid - 1, k_min, k_best);
    dnc(mid + 1, r, k_best, k_max);
}

int main() {
    cin >> M;
    P = (M + 1) / 2;
    N = max(8000, min(20000, P * 200));

    y.resize(N + 1);
    v.resize(N + 1);
    s.resize(N + 1);
    dp0.assign(N + 1, INF);
    dp1.resize(N + 1);

    double d = R / N;
    for (int i = 0; i <= N; ++i) {
        y[i] = i * d;
        double z = R * R - y[i] * y[i];
        v[i] = (z > 0) ? sqrt(z) : 0;
    }

    s[0] = 0.0;
    for (int i = 1; i <= N; ++i) {
        double arg = y[i] / R;
        if (arg > 1.0) arg = 1.0;
        if (arg < -1.0) arg = -1.0;
        s[i] = 0.5 * (y[i] * v[i] + R * R * asin(arg));
    }
    
    dp0[0] = 0.0;

    for (int i = 1; i <= P; ++i) {
        fill(dp1.begin(), dp1.end(), 0.0);
        dnc(1, N, 0, N - 1);
        dp0 = dp1;
    }

    double ans = 4 * (s[N] + dp0[N]);
    cout << fixed << setprecision(4) << ans << endl;

    return 0;
}

import java.util.Scanner;
import java.util.Arrays;

public class Main {
    static final double R = 0.5;
    static final double INF = 1e100;
    static int M, P, N;
    static double[] y, v, s;
    static double[] dp0, dp1;

    public static void main(String[] args) {
        Scanner sc = new Scanner(System.in);
        M = sc.nextInt();
        P = (M + 1) / 2;
        N = Math.max(8000, Math.min(20000, P * 200));

        y = new double[N + 1];
        v = new double[N + 1];
        s = new double[N + 1];
        dp0 = new double[N + 1];
        dp1 = new double[N + 1];

        double d = R / N;
        for (int i = 0; i <= N; i++) {
            y[i] = i * d;
            double z = R * R - y[i] * y[i];
            v[i] = (z > 0) ? Math.sqrt(z) : 0;
        }

        s[0] = 0.0;
        for (int i = 1; i <= N; i++) {
            double arg = y[i] / R;
            if (arg > 1.0) arg = 1.0;
            if (arg < -1.0) arg = -1.0;
            s[i] = 0.5 * (y[i] * v[i] + R * R * Math.asin(arg));
        }

        Arrays.fill(dp0, INF);
        dp0[0] = 0.0;

        for (int i = 1; i <= P; i++) {
            Arrays.fill(dp1, 0.0);
            dnc(1, N, 0, N - 1);
            dp0 = dp1.clone();
        }

        double ans = 4 * (s[N] + dp0[N]);
        System.out.printf("%.4f\n", ans);
    }
    
    private static double cost(int k, int j) {
        return (y[j] - y[k]) * v[k] - (s[j] - s[k]);
    }

    private static void dnc(int l, int r, int k_min, int k_max) {
        if (l > r) return;

        int mid = l + (r - l) / 2;
        double best_val = INF;
        int k_best = -1;

        int k_upper = Math.min(k_max, mid - 1);
        for (int k = k_min; k <= k_upper; k++) {
            double current_val = dp0[k] + cost(k, mid);
            if (current_val < best_val) {
                best_val = current_val;
                k_best = k;
            }
        }
        dp1[mid] = best_val;

        dnc(l, mid - 1, k_min, k_best);
        dnc(mid + 1, r, k_best, k_max);
    }
}

import sys
import math

# 增大递归深度限制
sys.setrecursionlimit(20000)

def solve():
    M = int(input())

    R = 0.5
    P = (M + 1) // 2
    N = max(8000, min(20000, P * 200))
    INF = 1e100

    y = [0.0] * (N + 1)
    v = [0.0] * (N + 1)
    s = [0.0] * (N + 1)
    
    d = R / N
    for i in range(N + 1):
        y[i] = i * d
        z = R * R - y[i] * y[i]
        v[i] = math.sqrt(z) if z > 0 else 0

    s[0] = 0.0
    for i in range(1, N + 1):
        # 边界情况处理
        arg = y[i] / R
        if arg > 1.0: arg = 1.0
        if arg < -1.0: arg = -1.0
        s[i] = 0.5 * (y[i] * v[i] + R * R * math.asin(arg))

    dp0 = [INF] * (N + 1)
    dp0[0] = 0.0
    dp1 = [0.0] * (N + 1)

    def cost(k, j):
        return (y[j] - y[k]) * v[k] - (s[j] - s[k])

    def dnc(l, r, k_min, k_max):
        if l > r:
            return

        mid = (l + r) // 2
        best_val = INF
        k_best = -1

        k_upper = min(k_max, mid - 1)
        for k in range(k_min, k_upper + 1):
            current_val = dp0[k] + cost(k, mid)
            if current_val < best_val:
                best_val = current_val
                k_best = k
        
        dp1[mid] = best_val

        dnc(l, mid - 1, k_min, k_best)
        dnc(mid + 1, r, k_best, k_max)

    for _ in range(P):
        dnc(1, N, 0, N - 1)
        dp0 = list(dp1)

    ans = 4 * (s[N] + dp0[N])
    print(f"{ans:.4f}")

solve()

算法及复杂度

  • 算法: 动态规划 (DP) + 分治优化 (D&C Optimization)
  • 时间复杂度: , 其中 是切割线数量 (), 是离散化点的数量 (一个较大的常数，如 8000-20000)。 次迭代，每次分治优化的 DP 计算为 。
  • 空间复杂度: , 主要用于存储 DP 数组和预计算的 数组。