UVA 12457：橢圓周長計算與數值積分優化

問題描述

UVA 12457 要求計算橢圓的周長。給定長軸 a 和短軸 b，需輸出周長，精度至少達到小數點後 6 位。

解題思路

方法一：積分 + 羅姆貝格加速

橢圓周長可用積分公式表示：

L = 4a ∫₀^{π/2} √(1 - e² sin²θ) dθ，其中離心率 e² = (a² - b²)/a²。

直接數值積分（如矩形法）需極小的步長才能達到 1e-6 精度，易導致 TLE。因此引入羅姆貝格積分（Romberg Integration），利用理查森外推加速收斂，大幅減少函數評估次數。

核心步驟：

1. 以梯形法則計算初步積分值 T₁。
2. 反覆將步長減半，計算 T₂, T₃ …
3. 應用外推公式生成 S、C、R 列，直到相邻 R 值差值小於給定容差（如 5e-7）。

代碼採用遞推陣列 T, S, C, R 實現，最終返回 R 列的最新值作為積分結果，再乘以 4a 即得周長。

方法二：級數展開（更簡潔高效）

利用橢圓周長的無窮級數展開：

L = 2πa [1 - Σ_{n=1}^{∞} ( (2n-1)!! / (2n)!! )² * e^{2n} / (2n-1) ]

其中 e² = 1 - b²/a²。只需累加足夠多項（如 100,000 項）即可滿足精度要求，代碼極短且運行極快。

關鍵細節

• π 的取法：acos(-1.0) 或 4.0*atan(1.0)，均可獲得雙精度 π。
• 輸入順序：題目給定 b a，注意讀入時對應長短軸。
• 輸出格式：Case %d: %.6lf，取絕對值以防負零。

原始碼參考

羅姆貝格積分版

[cpp]
#include <iostream>
#include <stdio.h>
#include <math.h>

using namespace std;
#define Pi 3.1415926535898
#define ds 0.0000005
#define M 40

double a = 0, b = 0;

double f(double si) {
    return sqrt(1 - ((a*a - b*b) / (a*a)) * sin(si) * sin(si));
}

long Sq(int k) {
    long m = 1;
    for (; k > 0; k--) m *= 2;
    return m;
}

double rom(double a, double b, double e) {
    double T[M], S[M], C[M], R[M], y;
    int k = 0, i;
    double h = (b - a) * 1.0;
    T[1] = h / 2.0 * (f(a) + f(b));
    h /= 2.0;

    for (k = 1; k <= M; k++, h /= 2) {
        y = 0;
        for (i = 1; i <= Sq(k - 1); i++) {
            y += f(a + (2 * i - 1) * h);
        }
        T[Sq(k)] = T[Sq(k - 1)] / 2.0 + h * y;
        S[Sq(k - 1)] = 4.0 / 3 * T[Sq(k)] - 1.0 / 3 * T[Sq(k - 1)];
        if (k - 2 >= 0)
            C[Sq(k - 2)] = 16.0 / 15 * S[Sq(k - 1)] - 1.0 / 15 * S[Sq(k - 2)];
        if (k - 3 >= 0)
            R[Sq(k - 3)] = 64.0 / 63 * C[Sq(k - 2)] - 1.0 / 63 * C[Sq(k - 3)];
        if (k - 4 >= 0) {
            if (fabs(R[Sq(k - 3)] - R[Sq(k - 4)]) < e)
                break;
        }
    }
    return R[Sq(k - 3)];
}

int main() {
    int t = 0, k = 0;
    double ans = 0;
    scanf("%d", &t);
    for (k = 1; k <= t; k++) {
        scanf("%lf %lf", &b, &a);
        ans = 4 * a * rom(Pi / 2, 0, ds);
        printf("Case %d: %.6lf\n", k, fabs(ans));
    }
    return 0;
}

級數展開版

[cpp]
#include <stdio.h>
#include <string.h>
#include <math.h>
const double pi = (4.0 * atan(1.0));
double a, b;

int main() {
    int i, j, T, cas = 0;
    scanf("%d", &T);
    while (T--) {
        scanf("%lf%lf", &b, &a);
        double e = 1.0 - (b * b) / (a * a);
        double sum = 1.0, tmp = 1.0;
        for (i = 1; i < 100000; i++) {
            double r = (2.0 * i - 1) / (2.0 * i);
            sum -= (tmp *= r * r * e) / (2.0 * i - 1);
        }
        printf("Case %d: %.8f\n", ++cas, 2.0 * pi * a * sum);
    }
    return 0;
}

總結

本題展示了選擇合適演算法的重要性：同樣求解橢圓周長，羅姆貝格積分需較多代碼且易調試；級數展開則代碼極簡、運行極快。實戰中應優先考慮數學閉式解或快速收斂級數。